TESTING/LIN/ctpt01.f

   1 *> \brief \b CTPT01
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE CTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       CHARACTER          DIAG, UPLO
  15 *       INTEGER            N
  16 *       REAL               RCOND, RESID
  17 *       ..
  18 *       .. Array Arguments ..
  19 *       REAL               RWORK( * )
  20 *       COMPLEX            AINVP( * ), AP( * )
  21 *       ..
  22 *
  23 *
  24 *> \par Purpose:
  25 *  =============
  26 *>
  27 *> \verbatim
  28 *>
  29 *> CTPT01 computes the residual for a triangular matrix A times its
  30 *> inverse when A is stored in packed format:
  31 *>    RESID = norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ),
  32 *> where EPS is the machine epsilon.
  33 *> \endverbatim
  34 *
  35 *  Arguments:
  36 *  ==========
  37 *
  38 *> \param[in] UPLO
  39 *> \verbatim
  40 *>          UPLO is CHARACTER*1
  41 *>          Specifies whether the matrix A is upper or lower triangular.
  42 *>          = 'U':  Upper triangular
  43 *>          = 'L':  Lower triangular
  44 *> \endverbatim
  45 *>
  46 *> \param[in] DIAG
  47 *> \verbatim
  48 *>          DIAG is CHARACTER*1
  49 *>          Specifies whether or not the matrix A is unit triangular.
  50 *>          = 'N':  Non-unit triangular
  51 *>          = 'U':  Unit triangular
  52 *> \endverbatim
  53 *>
  54 *> \param[in] N
  55 *> \verbatim
  56 *>          N is INTEGER
  57 *>          The order of the matrix A.  N >= 0.
  58 *> \endverbatim
  59 *>
  60 *> \param[in] AP
  61 *> \verbatim
  62 *>          AP is COMPLEX array, dimension (N*(N+1)/2)
  63 *>          The original upper or lower triangular matrix A, packed
  64 *>          columnwise in a linear array.  The j-th column of A is stored
  65 *>          in the array AP as follows:
  66 *>          if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
  67 *>          if UPLO = 'L',
  68 *>             AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] AINVP
  72 *> \verbatim
  73 *>          AINVP is COMPLEX array, dimension (N*(N+1)/2)
  74 *>          On entry, the (triangular) inverse of the matrix A, packed
  75 *>          columnwise in a linear array as in AP.
  76 *>          On exit, the contents of AINVP are destroyed.
  77 *> \endverbatim
  78 *>
  79 *> \param[out] RCOND
  80 *> \verbatim
  81 *>          RCOND is REAL
  82 *>          The reciprocal condition number of A, computed as
  83 *>          1/(norm(A) * norm(AINV)).
  84 *> \endverbatim
  85 *>
  86 *> \param[out] RWORK
  87 *> \verbatim
  88 *>          RWORK is REAL array, dimension (N)
  89 *> \endverbatim
  90 *>
  91 *> \param[out] RESID
  92 *> \verbatim
  93 *>          RESID is REAL
  94 *>          norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS )
  95 *> \endverbatim
  96 *
  97 *  Authors:
  98 *  ========
  99 *
 100 *> \author Univ. of Tennessee
 101 *> \author Univ. of California Berkeley
 102 *> \author Univ. of Colorado Denver
 103 *> \author NAG Ltd.
 104 *
 105 *> \date November 2011
 106 *
 107 *> \ingroup complex_lin
 108 *
 109 *  =====================================================================
 110       SUBROUTINE CTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, RWORK, RESID )
 111 *
 112 *  -- LAPACK test routine (version 3.4.0) --
 113 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 114 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 115 *     November 2011
 116 *
 117 *     .. Scalar Arguments ..
 118       CHARACTER          DIAG, UPLO
 119       INTEGER            N
 120       REAL               RCOND, RESID
 121 *     ..
 122 *     .. Array Arguments ..
 123       REAL               RWORK( * )
 124       COMPLEX            AINVP( * ), AP( * )
 125 *     ..
 126 *
 127 *  =====================================================================
 128 *
 129 *     .. Parameters ..
 130       REAL               ZERO, ONE
 131       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 132 *     ..
 133 *     .. Local Scalars ..
 134       LOGICAL            UNITD
 135       INTEGER            J, JC
 136       REAL               AINVNM, ANORM, EPS
 137 *     ..
 138 *     .. External Functions ..
 139       LOGICAL            LSAME
 140       REAL               CLANTP, SLAMCH
 141       EXTERNAL           LSAME, CLANTP, SLAMCH
 142 *     ..
 143 *     .. External Subroutines ..
 144       EXTERNAL           CTPMV
 145 *     ..
 146 *     .. Intrinsic Functions ..
 147       INTRINSIC          REAL
 148 *     ..
 149 *     .. Executable Statements ..
 150 *
 151 *     Quick exit if N = 0.
 152 *
 153       IF( N.LE.0 ) THEN
 154          RCOND = ONE
 155          RESID = ZERO
 156          RETURN
 157       END IF
 158 *
 159 *     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
 160 *
 161       EPS = SLAMCH( 'Epsilon' )
 162       ANORM = CLANTP( '1', UPLO, DIAG, N, AP, RWORK )
 163       AINVNM = CLANTP( '1', UPLO, DIAG, N, AINVP, RWORK )
 164       IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
 165          RCOND = ZERO
 166          RESID = ONE / EPS
 167          RETURN
 168       END IF
 169       RCOND = ( ONE / ANORM ) / AINVNM
 170 *
 171 *     Compute A * AINV, overwriting AINV.
 172 *
 173       UNITD = LSAME( DIAG, 'U' )
 174       IF( LSAME( UPLO, 'U' ) ) THEN
 175          JC = 1
 176          DO 10 J = 1, N
 177             IF( UNITD )
 178      $         AINVP( JC+J-1 ) = ONE
 179 *
 180 *           Form the j-th column of A*AINV.
 181 *
 182             CALL CTPMV( 'Upper', 'No transpose', DIAG, J, AP,
 183      $                  AINVP( JC ), 1 )
 184 *
 185 *           Subtract 1 from the diagonal to form A*AINV - I.
 186 *
 187             AINVP( JC+J-1 ) = AINVP( JC+J-1 ) - ONE
 188             JC = JC + J
 189    10    CONTINUE
 190       ELSE
 191          JC = 1
 192          DO 20 J = 1, N
 193             IF( UNITD )
 194      $         AINVP( JC ) = ONE
 195 *
 196 *           Form the j-th column of A*AINV.
 197 *
 198             CALL CTPMV( 'Lower', 'No transpose', DIAG, N-J+1, AP( JC ),
 199      $                  AINVP( JC ), 1 )
 200 *
 201 *           Subtract 1 from the diagonal to form A*AINV - I.
 202 *
 203             AINVP( JC ) = AINVP( JC ) - ONE
 204             JC = JC + N - J + 1
 205    20    CONTINUE
 206       END IF
 207 *
 208 *     Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS)
 209 *
 210       RESID = CLANTP( '1', UPLO, 'Non-unit', N, AINVP, RWORK )
 211 *
 212       RESID = ( ( RESID*RCOND ) / REAL( N ) ) / EPS
 213 *
 214       RETURN
 215 *
 216 *     End of CTPT01
 217 *
 218       END