1 SUBROUTINE ZSPSVX( 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.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1999 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* ZSPSVX uses the diagonal pivoting factorization A = U*D*U**T or 25* A = L*D*L**T to compute the solution to a complex system of linear 26* equations A * X = B, where A is an N-by-N symmetric 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**T, if UPLO = 'U', or 39* A = L * D * L**T, if UPLO = 'L', 40* where U (or L) is a product of permutation and unit upper (lower) 41* triangular matrices and D is symmetric 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 65* of A. AP, 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 symmetric 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**T or A = L*D*L**T as computed by ZSPTRF, 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**T or A = L*D*L**T as computed by ZSPTRF, 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 ZSPTRF. 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 ZSPTRF. 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 symmetric matrix A: 178* 179* a11 a12 a13 a14 180* a22 a23 a24 181* a33 a34 (aij = 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, ZLANSP 201 EXTERNAL LSAME, DLAMCH, ZLANSP 202* .. 203* .. External Subroutines .. 204 EXTERNAL XERBLA, ZCOPY, ZLACPY, ZSPCON, ZSPRFS, ZSPTRF, 205 $ ZSPTRS 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( 'ZSPSVX', -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 ZSPTRF( UPLO, N, AFP, IPIV, INFO ) 241* 242* Return if INFO is non-zero. 243* 244 IF( INFO.NE.0 ) THEN 245 IF( INFO.GT.0 ) 246 $ RCOND = ZERO 247 RETURN 248 END IF 249 END IF 250* 251* Compute the norm of the matrix A. 252* 253 ANORM = ZLANSP( 'I', UPLO, N, AP, RWORK ) 254* 255* Compute the reciprocal of the condition number of A. 256* 257 CALL ZSPCON( UPLO, N, AFP, IPIV, ANORM, RCOND, WORK, INFO ) 258* 259* Set INFO = N+1 if the matrix is singular to working precision. 260* 261 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 262 $ INFO = N + 1 263* 264* Compute the solution vectors X. 265* 266 CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 267 CALL ZSPTRS( UPLO, N, NRHS, AFP, IPIV, X, LDX, INFO ) 268* 269* Use iterative refinement to improve the computed solutions and 270* compute error bounds and backward error estimates for them. 271* 272 CALL ZSPRFS( UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, FERR, 273 $ BERR, WORK, RWORK, INFO ) 274* 275 RETURN 276* 277* End of ZSPSVX 278* 279 END 280