00001 SUBROUTINE ZSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, 00002 $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, 00003 $ RWORK, INFO ) 00004 * 00005 * -- LAPACK driver routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER FACT, UPLO 00012 INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS 00013 DOUBLE PRECISION RCOND 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER IPIV( * ) 00017 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ) 00018 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), 00019 $ WORK( * ), X( LDX, * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * ZSYSVX uses the diagonal pivoting factorization to compute the 00026 * solution to a complex system of linear equations A * X = B, 00027 * where A is an N-by-N symmetric matrix and X and B are N-by-NRHS 00028 * matrices. 00029 * 00030 * Error bounds on the solution and a condition estimate are also 00031 * provided. 00032 * 00033 * Description 00034 * =========== 00035 * 00036 * The following steps are performed: 00037 * 00038 * 1. If FACT = 'N', the diagonal pivoting method is used to factor A. 00039 * The form of the factorization is 00040 * A = U * D * U**T, if UPLO = 'U', or 00041 * A = L * D * L**T, if UPLO = 'L', 00042 * where U (or L) is a product of permutation and unit upper (lower) 00043 * triangular matrices, and D is symmetric and block diagonal with 00044 * 1-by-1 and 2-by-2 diagonal blocks. 00045 * 00046 * 2. If some D(i,i)=0, so that D is exactly singular, then the routine 00047 * returns with INFO = i. Otherwise, the factored form of A is used 00048 * to estimate the condition number of the matrix A. If the 00049 * reciprocal of the condition number is less than machine precision, 00050 * INFO = N+1 is returned as a warning, but the routine still goes on 00051 * to solve for X and compute error bounds as described below. 00052 * 00053 * 3. The system of equations is solved for X using the factored form 00054 * of A. 00055 * 00056 * 4. Iterative refinement is applied to improve the computed solution 00057 * matrix and calculate error bounds and backward error estimates 00058 * for it. 00059 * 00060 * Arguments 00061 * ========= 00062 * 00063 * FACT (input) CHARACTER*1 00064 * Specifies whether or not the factored form of A has been 00065 * supplied on entry. 00066 * = 'F': On entry, AF and IPIV contain the factored form 00067 * of A. A, AF and IPIV will not be modified. 00068 * = 'N': The matrix A will be copied to AF and factored. 00069 * 00070 * UPLO (input) CHARACTER*1 00071 * = 'U': Upper triangle of A is stored; 00072 * = 'L': Lower triangle of A is stored. 00073 * 00074 * N (input) INTEGER 00075 * The number of linear equations, i.e., the order of the 00076 * matrix A. N >= 0. 00077 * 00078 * NRHS (input) INTEGER 00079 * The number of right hand sides, i.e., the number of columns 00080 * of the matrices B and X. NRHS >= 0. 00081 * 00082 * A (input) COMPLEX*16 array, dimension (LDA,N) 00083 * The symmetric matrix A. If UPLO = 'U', the leading N-by-N 00084 * upper triangular part of A contains the upper triangular part 00085 * of the matrix A, and the strictly lower triangular part of A 00086 * is not referenced. If UPLO = 'L', the leading N-by-N lower 00087 * triangular part of A contains the lower triangular part of 00088 * the matrix A, and the strictly upper triangular part of A is 00089 * not referenced. 00090 * 00091 * LDA (input) INTEGER 00092 * The leading dimension of the array A. LDA >= max(1,N). 00093 * 00094 * AF (input or output) COMPLEX*16 array, dimension (LDAF,N) 00095 * If FACT = 'F', then AF is an input argument and on entry 00096 * contains the block diagonal matrix D and the multipliers used 00097 * to obtain the factor U or L from the factorization 00098 * A = U*D*U**T or A = L*D*L**T as computed by ZSYTRF. 00099 * 00100 * If FACT = 'N', then AF is an output argument and on exit 00101 * returns the block diagonal matrix D and the multipliers used 00102 * to obtain the factor U or L from the factorization 00103 * A = U*D*U**T or A = L*D*L**T. 00104 * 00105 * LDAF (input) INTEGER 00106 * The leading dimension of the array AF. LDAF >= max(1,N). 00107 * 00108 * IPIV (input or output) INTEGER array, dimension (N) 00109 * If FACT = 'F', then IPIV is an input argument and on entry 00110 * contains details of the interchanges and the block structure 00111 * of D, as determined by ZSYTRF. 00112 * If IPIV(k) > 0, then rows and columns k and IPIV(k) were 00113 * interchanged and D(k,k) is a 1-by-1 diagonal block. 00114 * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and 00115 * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) 00116 * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = 00117 * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were 00118 * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 00119 * 00120 * If FACT = 'N', then IPIV is an output argument and on exit 00121 * contains details of the interchanges and the block structure 00122 * of D, as determined by ZSYTRF. 00123 * 00124 * B (input) COMPLEX*16 array, dimension (LDB,NRHS) 00125 * The N-by-NRHS right hand side matrix B. 00126 * 00127 * LDB (input) INTEGER 00128 * The leading dimension of the array B. LDB >= max(1,N). 00129 * 00130 * X (output) COMPLEX*16 array, dimension (LDX,NRHS) 00131 * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X. 00132 * 00133 * LDX (input) INTEGER 00134 * The leading dimension of the array X. LDX >= max(1,N). 00135 * 00136 * RCOND (output) DOUBLE PRECISION 00137 * The estimate of the reciprocal condition number of the matrix 00138 * A. If RCOND is less than the machine precision (in 00139 * particular, if RCOND = 0), the matrix is singular to working 00140 * precision. This condition is indicated by a return code of 00141 * INFO > 0. 00142 * 00143 * FERR (output) DOUBLE PRECISION array, dimension (NRHS) 00144 * The estimated forward error bound for each solution vector 00145 * X(j) (the j-th column of the solution matrix X). 00146 * If XTRUE is the true solution corresponding to X(j), FERR(j) 00147 * is an estimated upper bound for the magnitude of the largest 00148 * element in (X(j) - XTRUE) divided by the magnitude of the 00149 * largest element in X(j). The estimate is as reliable as 00150 * the estimate for RCOND, and is almost always a slight 00151 * overestimate of the true error. 00152 * 00153 * BERR (output) DOUBLE PRECISION array, dimension (NRHS) 00154 * The componentwise relative backward error of each solution 00155 * vector X(j) (i.e., the smallest relative change in 00156 * any element of A or B that makes X(j) an exact solution). 00157 * 00158 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 00159 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00160 * 00161 * LWORK (input) INTEGER 00162 * The length of WORK. LWORK >= max(1,2*N), and for best 00163 * performance, when FACT = 'N', LWORK >= max(1,2*N,N*NB), where 00164 * NB is the optimal blocksize for ZSYTRF. 00165 * 00166 * If LWORK = -1, then a workspace query is assumed; the routine 00167 * only calculates the optimal size of the WORK array, returns 00168 * this value as the first entry of the WORK array, and no error 00169 * message related to LWORK is issued by XERBLA. 00170 * 00171 * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 00172 * 00173 * INFO (output) INTEGER 00174 * = 0: successful exit 00175 * < 0: if INFO = -i, the i-th argument had an illegal value 00176 * > 0: if INFO = i, and i is 00177 * <= N: D(i,i) is exactly zero. The factorization 00178 * has been completed but the factor D is exactly 00179 * singular, so the solution and error bounds could 00180 * not be computed. RCOND = 0 is returned. 00181 * = N+1: D is nonsingular, but RCOND is less than machine 00182 * precision, meaning that the matrix is singular 00183 * to working precision. Nevertheless, the 00184 * solution and error bounds are computed because 00185 * there are a number of situations where the 00186 * computed solution can be more accurate than the 00187 * value of RCOND would suggest. 00188 * 00189 * ===================================================================== 00190 * 00191 * .. Parameters .. 00192 DOUBLE PRECISION ZERO 00193 PARAMETER ( ZERO = 0.0D+0 ) 00194 * .. 00195 * .. Local Scalars .. 00196 LOGICAL LQUERY, NOFACT 00197 INTEGER LWKOPT, NB 00198 DOUBLE PRECISION ANORM 00199 * .. 00200 * .. External Functions .. 00201 LOGICAL LSAME 00202 INTEGER ILAENV 00203 DOUBLE PRECISION DLAMCH, ZLANSY 00204 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANSY 00205 * .. 00206 * .. External Subroutines .. 00207 EXTERNAL XERBLA, ZLACPY, ZSYCON, ZSYRFS, ZSYTRF, ZSYTRS 00208 * .. 00209 * .. Intrinsic Functions .. 00210 INTRINSIC MAX 00211 * .. 00212 * .. Executable Statements .. 00213 * 00214 * Test the input parameters. 00215 * 00216 INFO = 0 00217 NOFACT = LSAME( FACT, 'N' ) 00218 LQUERY = ( LWORK.EQ.-1 ) 00219 IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN 00220 INFO = -1 00221 ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) 00222 $ THEN 00223 INFO = -2 00224 ELSE IF( N.LT.0 ) THEN 00225 INFO = -3 00226 ELSE IF( NRHS.LT.0 ) THEN 00227 INFO = -4 00228 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00229 INFO = -6 00230 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN 00231 INFO = -8 00232 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00233 INFO = -11 00234 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00235 INFO = -13 00236 ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN 00237 INFO = -18 00238 END IF 00239 * 00240 IF( INFO.EQ.0 ) THEN 00241 LWKOPT = MAX( 1, 2*N ) 00242 IF( NOFACT ) THEN 00243 NB = ILAENV( 1, 'ZSYTRF', UPLO, N, -1, -1, -1 ) 00244 LWKOPT = MAX( LWKOPT, N*NB ) 00245 END IF 00246 WORK( 1 ) = LWKOPT 00247 END IF 00248 * 00249 IF( INFO.NE.0 ) THEN 00250 CALL XERBLA( 'ZSYSVX', -INFO ) 00251 RETURN 00252 ELSE IF( LQUERY ) THEN 00253 RETURN 00254 END IF 00255 * 00256 IF( NOFACT ) THEN 00257 * 00258 * Compute the factorization A = U*D*U' or A = L*D*L'. 00259 * 00260 CALL ZLACPY( UPLO, N, N, A, LDA, AF, LDAF ) 00261 CALL ZSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO ) 00262 * 00263 * Return if INFO is non-zero. 00264 * 00265 IF( INFO.GT.0 )THEN 00266 RCOND = ZERO 00267 RETURN 00268 END IF 00269 END IF 00270 * 00271 * Compute the norm of the matrix A. 00272 * 00273 ANORM = ZLANSY( 'I', UPLO, N, A, LDA, RWORK ) 00274 * 00275 * Compute the reciprocal of the condition number of A. 00276 * 00277 CALL ZSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, INFO ) 00278 * 00279 * Compute the solution vectors X. 00280 * 00281 CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX ) 00282 CALL ZSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO ) 00283 * 00284 * Use iterative refinement to improve the computed solutions and 00285 * compute error bounds and backward error estimates for them. 00286 * 00287 CALL ZSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, 00288 $ LDX, FERR, BERR, WORK, RWORK, INFO ) 00289 * 00290 * Set INFO = N+1 if the matrix is singular to working precision. 00291 * 00292 IF( RCOND.LT.DLAMCH( 'Epsilon' ) ) 00293 $ INFO = N + 1 00294 * 00295 WORK( 1 ) = LWKOPT 00296 * 00297 RETURN 00298 * 00299 * End of ZSYSVX 00300 * 00301 END
1.7.2