Annotation of rpl/lapack/lapack/zhpsvx.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZHPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X,
! 2: $ LDX, RCOND, FERR, BERR, WORK, RWORK, INFO )
! 3: *
! 4: * -- LAPACK driver routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER FACT, UPLO
! 11: INTEGER INFO, LDB, LDX, N, NRHS
! 12: DOUBLE PRECISION RCOND
! 13: * ..
! 14: * .. Array Arguments ..
! 15: INTEGER IPIV( * )
! 16: DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
! 17: COMPLEX*16 AFP( * ), AP( * ), B( LDB, * ), WORK( * ),
! 18: $ X( LDX, * )
! 19: * ..
! 20: *
! 21: * Purpose
! 22: * =======
! 23: *
! 24: * ZHPSVX uses the diagonal pivoting factorization A = U*D*U**H or
! 25: * A = L*D*L**H to compute the solution to a complex system of linear
! 26: * equations A * X = B, where A is an N-by-N Hermitian matrix stored
! 27: * in packed format and X and B are N-by-NRHS matrices.
! 28: *
! 29: * Error bounds on the solution and a condition estimate are also
! 30: * provided.
! 31: *
! 32: * Description
! 33: * ===========
! 34: *
! 35: * The following steps are performed:
! 36: *
! 37: * 1. If FACT = 'N', the diagonal pivoting method is used to factor A as
! 38: * A = U * D * U**H, if UPLO = 'U', or
! 39: * A = L * D * L**H, if UPLO = 'L',
! 40: * where U (or L) is a product of permutation and unit upper (lower)
! 41: * triangular matrices and D is Hermitian and block diagonal with
! 42: * 1-by-1 and 2-by-2 diagonal blocks.
! 43: *
! 44: * 2. If some D(i,i)=0, so that D is exactly singular, then the routine
! 45: * returns with INFO = i. Otherwise, the factored form of A is used
! 46: * to estimate the condition number of the matrix A. If the
! 47: * reciprocal of the condition number is less than machine precision,
! 48: * INFO = N+1 is returned as a warning, but the routine still goes on
! 49: * to solve for X and compute error bounds as described below.
! 50: *
! 51: * 3. The system of equations is solved for X using the factored form
! 52: * of A.
! 53: *
! 54: * 4. Iterative refinement is applied to improve the computed solution
! 55: * matrix and calculate error bounds and backward error estimates
! 56: * for it.
! 57: *
! 58: * Arguments
! 59: * =========
! 60: *
! 61: * FACT (input) CHARACTER*1
! 62: * Specifies whether or not the factored form of A has been
! 63: * supplied on entry.
! 64: * = 'F': On entry, AFP and IPIV contain the factored form of
! 65: * A. AFP and IPIV will not be modified.
! 66: * = 'N': The matrix A will be copied to AFP and factored.
! 67: *
! 68: * UPLO (input) CHARACTER*1
! 69: * = 'U': Upper triangle of A is stored;
! 70: * = 'L': Lower triangle of A is stored.
! 71: *
! 72: * N (input) INTEGER
! 73: * The number of linear equations, i.e., the order of the
! 74: * matrix A. N >= 0.
! 75: *
! 76: * NRHS (input) INTEGER
! 77: * The number of right hand sides, i.e., the number of columns
! 78: * of the matrices B and X. NRHS >= 0.
! 79: *
! 80: * AP (input) COMPLEX*16 array, dimension (N*(N+1)/2)
! 81: * The upper or lower triangle of the Hermitian matrix A, packed
! 82: * columnwise in a linear array. The j-th column of A is stored
! 83: * in the array AP as follows:
! 84: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
! 85: * if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
! 86: * See below for further details.
! 87: *
! 88: * AFP (input or output) COMPLEX*16 array, dimension (N*(N+1)/2)
! 89: * If FACT = 'F', then AFP is an input argument and on entry
! 90: * contains the block diagonal matrix D and the multipliers used
! 91: * to obtain the factor U or L from the factorization
! 92: * A = U*D*U**H or A = L*D*L**H as computed by ZHPTRF, stored as
! 93: * a packed triangular matrix in the same storage format as A.
! 94: *
! 95: * If FACT = 'N', then AFP is an output argument and on exit
! 96: * contains the block diagonal matrix D and the multipliers used
! 97: * to obtain the factor U or L from the factorization
! 98: * A = U*D*U**H or A = L*D*L**H as computed by ZHPTRF, stored as
! 99: * a packed triangular matrix in the same storage format as A.
! 100: *
! 101: * IPIV (input or output) INTEGER array, dimension (N)
! 102: * If FACT = 'F', then IPIV is an input argument and on entry
! 103: * contains details of the interchanges and the block structure
! 104: * of D, as determined by ZHPTRF.
! 105: * If IPIV(k) > 0, then rows and columns k and IPIV(k) were
! 106: * interchanged and D(k,k) is a 1-by-1 diagonal block.
! 107: * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
! 108: * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
! 109: * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
! 110: * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
! 111: * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
! 112: *
! 113: * If FACT = 'N', then IPIV is an output argument and on exit
! 114: * contains details of the interchanges and the block structure
! 115: * of D, as determined by ZHPTRF.
! 116: *
! 117: * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
! 118: * The N-by-NRHS right hand side matrix B.
! 119: *
! 120: * LDB (input) INTEGER
! 121: * The leading dimension of the array B. LDB >= max(1,N).
! 122: *
! 123: * X (output) COMPLEX*16 array, dimension (LDX,NRHS)
! 124: * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
! 125: *
! 126: * LDX (input) INTEGER
! 127: * The leading dimension of the array X. LDX >= max(1,N).
! 128: *
! 129: * RCOND (output) DOUBLE PRECISION
! 130: * The estimate of the reciprocal condition number of the matrix
! 131: * A. If RCOND is less than the machine precision (in
! 132: * particular, if RCOND = 0), the matrix is singular to working
! 133: * precision. This condition is indicated by a return code of
! 134: * INFO > 0.
! 135: *
! 136: * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
! 137: * The estimated forward error bound for each solution vector
! 138: * X(j) (the j-th column of the solution matrix X).
! 139: * If XTRUE is the true solution corresponding to X(j), FERR(j)
! 140: * is an estimated upper bound for the magnitude of the largest
! 141: * element in (X(j) - XTRUE) divided by the magnitude of the
! 142: * largest element in X(j). The estimate is as reliable as
! 143: * the estimate for RCOND, and is almost always a slight
! 144: * overestimate of the true error.
! 145: *
! 146: * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
! 147: * The componentwise relative backward error of each solution
! 148: * vector X(j) (i.e., the smallest relative change in
! 149: * any element of A or B that makes X(j) an exact solution).
! 150: *
! 151: * WORK (workspace) COMPLEX*16 array, dimension (2*N)
! 152: *
! 153: * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
! 154: *
! 155: * INFO (output) INTEGER
! 156: * = 0: successful exit
! 157: * < 0: if INFO = -i, the i-th argument had an illegal value
! 158: * > 0: if INFO = i, and i is
! 159: * <= N: D(i,i) is exactly zero. The factorization
! 160: * has been completed but the factor D is exactly
! 161: * singular, so the solution and error bounds could
! 162: * not be computed. RCOND = 0 is returned.
! 163: * = N+1: D is nonsingular, but RCOND is less than machine
! 164: * precision, meaning that the matrix is singular
! 165: * to working precision. Nevertheless, the
! 166: * solution and error bounds are computed because
! 167: * there are a number of situations where the
! 168: * computed solution can be more accurate than the
! 169: * value of RCOND would suggest.
! 170: *
! 171: * Further Details
! 172: * ===============
! 173: *
! 174: * The packed storage scheme is illustrated by the following example
! 175: * when N = 4, UPLO = 'U':
! 176: *
! 177: * Two-dimensional storage of the Hermitian matrix A:
! 178: *
! 179: * a11 a12 a13 a14
! 180: * a22 a23 a24
! 181: * a33 a34 (aij = conjg(aji))
! 182: * a44
! 183: *
! 184: * Packed storage of the upper triangle of A:
! 185: *
! 186: * AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
! 187: *
! 188: * =====================================================================
! 189: *
! 190: * .. Parameters ..
! 191: DOUBLE PRECISION ZERO
! 192: PARAMETER ( ZERO = 0.0D+0 )
! 193: * ..
! 194: * .. Local Scalars ..
! 195: LOGICAL NOFACT
! 196: DOUBLE PRECISION ANORM
! 197: * ..
! 198: * .. External Functions ..
! 199: LOGICAL LSAME
! 200: DOUBLE PRECISION DLAMCH, ZLANHP
! 201: EXTERNAL LSAME, DLAMCH, ZLANHP
! 202: * ..
! 203: * .. External Subroutines ..
! 204: EXTERNAL XERBLA, ZCOPY, ZHPCON, ZHPRFS, ZHPTRF, ZHPTRS,
! 205: $ ZLACPY
! 206: * ..
! 207: * .. Intrinsic Functions ..
! 208: INTRINSIC MAX
! 209: * ..
! 210: * .. Executable Statements ..
! 211: *
! 212: * Test the input parameters.
! 213: *
! 214: INFO = 0
! 215: NOFACT = LSAME( FACT, 'N' )
! 216: IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
! 217: INFO = -1
! 218: ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
! 219: $ THEN
! 220: INFO = -2
! 221: ELSE IF( N.LT.0 ) THEN
! 222: INFO = -3
! 223: ELSE IF( NRHS.LT.0 ) THEN
! 224: INFO = -4
! 225: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 226: INFO = -9
! 227: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
! 228: INFO = -11
! 229: END IF
! 230: IF( INFO.NE.0 ) THEN
! 231: CALL XERBLA( 'ZHPSVX', -INFO )
! 232: RETURN
! 233: END IF
! 234: *
! 235: IF( NOFACT ) THEN
! 236: *
! 237: * Compute the factorization A = U*D*U' or A = L*D*L'.
! 238: *
! 239: CALL ZCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 )
! 240: CALL ZHPTRF( UPLO, N, AFP, IPIV, INFO )
! 241: *
! 242: * Return if INFO is non-zero.
! 243: *
! 244: IF( INFO.GT.0 )THEN
! 245: RCOND = ZERO
! 246: RETURN
! 247: END IF
! 248: END IF
! 249: *
! 250: * Compute the norm of the matrix A.
! 251: *
! 252: ANORM = ZLANHP( 'I', UPLO, N, AP, RWORK )
! 253: *
! 254: * Compute the reciprocal of the condition number of A.
! 255: *
! 256: CALL ZHPCON( UPLO, N, AFP, IPIV, ANORM, RCOND, WORK, INFO )
! 257: *
! 258: * Compute the solution vectors X.
! 259: *
! 260: CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
! 261: CALL ZHPTRS( UPLO, N, NRHS, AFP, IPIV, X, LDX, INFO )
! 262: *
! 263: * Use iterative refinement to improve the computed solutions and
! 264: * compute error bounds and backward error estimates for them.
! 265: *
! 266: CALL ZHPRFS( UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, FERR,
! 267: $ BERR, WORK, RWORK, INFO )
! 268: *
! 269: * Set INFO = N+1 if the matrix is singular to working precision.
! 270: *
! 271: IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
! 272: $ INFO = N + 1
! 273: *
! 274: RETURN
! 275: *
! 276: * End of ZHPSVX
! 277: *
! 278: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhpsvx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack