00001 SUBROUTINE ZCPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK, 00002 + SWORK, RWORK, ITER, INFO ) 00003 * 00004 * -- LAPACK PROTOTYPE driver routine (version 3.3.0) -- 00005 * 00006 * November 2010 00007 * 00008 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00009 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00010 * .. 00011 * .. Scalar Arguments .. 00012 CHARACTER UPLO 00013 INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS 00014 * .. 00015 * .. Array Arguments .. 00016 DOUBLE PRECISION RWORK( * ) 00017 COMPLEX SWORK( * ) 00018 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ), 00019 + X( LDX, * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * ZCPOSV computes the solution to a complex system of linear equations 00026 * A * X = B, 00027 * where A is an N-by-N Hermitian positive definite matrix and X and B 00028 * are N-by-NRHS matrices. 00029 * 00030 * ZCPOSV first attempts to factorize the matrix in COMPLEX and use this 00031 * factorization within an iterative refinement procedure to produce a 00032 * solution with COMPLEX*16 normwise backward error quality (see below). 00033 * If the approach fails the method switches to a COMPLEX*16 00034 * factorization and solve. 00035 * 00036 * The iterative refinement is not going to be a winning strategy if 00037 * the ratio COMPLEX performance over COMPLEX*16 performance is too 00038 * small. A reasonable strategy should take the number of right-hand 00039 * sides and the size of the matrix into account. This might be done 00040 * with a call to ILAENV in the future. Up to now, we always try 00041 * iterative refinement. 00042 * 00043 * The iterative refinement process is stopped if 00044 * ITER > ITERMAX 00045 * or for all the RHS we have: 00046 * RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX 00047 * where 00048 * o ITER is the number of the current iteration in the iterative 00049 * refinement process 00050 * o RNRM is the infinity-norm of the residual 00051 * o XNRM is the infinity-norm of the solution 00052 * o ANRM is the infinity-operator-norm of the matrix A 00053 * o EPS is the machine epsilon returned by DLAMCH('Epsilon') 00054 * The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 00055 * respectively. 00056 * 00057 * Arguments 00058 * ========= 00059 * 00060 * UPLO (input) CHARACTER*1 00061 * = 'U': Upper triangle of A is stored; 00062 * = 'L': Lower triangle of A is stored. 00063 * 00064 * N (input) INTEGER 00065 * The number of linear equations, i.e., the order of the 00066 * matrix A. N >= 0. 00067 * 00068 * NRHS (input) INTEGER 00069 * The number of right hand sides, i.e., the number of columns 00070 * of the matrix B. NRHS >= 0. 00071 * 00072 * A (input/output) COMPLEX*16 array, 00073 * dimension (LDA,N) 00074 * On entry, the Hermitian matrix A. If UPLO = 'U', the leading 00075 * N-by-N upper triangular part of A contains the upper 00076 * triangular part of the matrix A, and the strictly lower 00077 * triangular part of A is not referenced. If UPLO = 'L', the 00078 * leading N-by-N lower triangular part of A contains the lower 00079 * triangular part of the matrix A, and the strictly upper 00080 * triangular part of A is not referenced. 00081 * 00082 * Note that the imaginary parts of the diagonal 00083 * elements need not be set and are assumed to be zero. 00084 * 00085 * On exit, if iterative refinement has been successfully used 00086 * (INFO.EQ.0 and ITER.GE.0, see description below), then A is 00087 * unchanged, if double precision factorization has been used 00088 * (INFO.EQ.0 and ITER.LT.0, see description below), then the 00089 * array A contains the factor U or L from the Cholesky 00090 * factorization A = U**H*U or A = L*L**H. 00091 * 00092 * LDA (input) INTEGER 00093 * The leading dimension of the array A. LDA >= max(1,N). 00094 * 00095 * B (input) COMPLEX*16 array, dimension (LDB,NRHS) 00096 * The N-by-NRHS right hand side matrix B. 00097 * 00098 * LDB (input) INTEGER 00099 * The leading dimension of the array B. LDB >= max(1,N). 00100 * 00101 * X (output) COMPLEX*16 array, dimension (LDX,NRHS) 00102 * If INFO = 0, the N-by-NRHS solution matrix X. 00103 * 00104 * LDX (input) INTEGER 00105 * The leading dimension of the array X. LDX >= max(1,N). 00106 * 00107 * WORK (workspace) COMPLEX*16 array, dimension (N*NRHS) 00108 * This array is used to hold the residual vectors. 00109 * 00110 * SWORK (workspace) COMPLEX array, dimension (N*(N+NRHS)) 00111 * This array is used to use the single precision matrix and the 00112 * right-hand sides or solutions in single precision. 00113 * 00114 * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 00115 * 00116 * ITER (output) INTEGER 00117 * < 0: iterative refinement has failed, COMPLEX*16 00118 * factorization has been performed 00119 * -1 : the routine fell back to full precision for 00120 * implementation- or machine-specific reasons 00121 * -2 : narrowing the precision induced an overflow, 00122 * the routine fell back to full precision 00123 * -3 : failure of CPOTRF 00124 * -31: stop the iterative refinement after the 30th 00125 * iterations 00126 * > 0: iterative refinement has been sucessfully used. 00127 * Returns the number of iterations 00128 * 00129 * INFO (output) INTEGER 00130 * = 0: successful exit 00131 * < 0: if INFO = -i, the i-th argument had an illegal value 00132 * > 0: if INFO = i, the leading minor of order i of 00133 * (COMPLEX*16) A is not positive definite, so the 00134 * factorization could not be completed, and the solution 00135 * has not been computed. 00136 * 00137 * ========= 00138 * 00139 * .. Parameters .. 00140 LOGICAL DOITREF 00141 PARAMETER ( DOITREF = .TRUE. ) 00142 * 00143 INTEGER ITERMAX 00144 PARAMETER ( ITERMAX = 30 ) 00145 * 00146 DOUBLE PRECISION BWDMAX 00147 PARAMETER ( BWDMAX = 1.0E+00 ) 00148 * 00149 COMPLEX*16 NEGONE, ONE 00150 PARAMETER ( NEGONE = ( -1.0D+00, 0.0D+00 ), 00151 + ONE = ( 1.0D+00, 0.0D+00 ) ) 00152 * 00153 * .. Local Scalars .. 00154 INTEGER I, IITER, PTSA, PTSX 00155 DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM 00156 COMPLEX*16 ZDUM 00157 * 00158 * .. External Subroutines .. 00159 EXTERNAL ZAXPY, ZHEMM, ZLACPY, ZLAT2C, ZLAG2C, CLAG2Z, 00160 + CPOTRF, CPOTRS, XERBLA 00161 * .. 00162 * .. External Functions .. 00163 INTEGER IZAMAX 00164 DOUBLE PRECISION DLAMCH, ZLANHE 00165 LOGICAL LSAME 00166 EXTERNAL IZAMAX, DLAMCH, ZLANHE, LSAME 00167 * .. 00168 * .. Intrinsic Functions .. 00169 INTRINSIC ABS, DBLE, MAX, SQRT 00170 * .. Statement Functions .. 00171 DOUBLE PRECISION CABS1 00172 * .. 00173 * .. Statement Function definitions .. 00174 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 00175 * .. 00176 * .. Executable Statements .. 00177 * 00178 INFO = 0 00179 ITER = 0 00180 * 00181 * Test the input parameters. 00182 * 00183 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00184 INFO = -1 00185 ELSE IF( N.LT.0 ) THEN 00186 INFO = -2 00187 ELSE IF( NRHS.LT.0 ) THEN 00188 INFO = -3 00189 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00190 INFO = -5 00191 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00192 INFO = -7 00193 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00194 INFO = -9 00195 END IF 00196 IF( INFO.NE.0 ) THEN 00197 CALL XERBLA( 'ZCPOSV', -INFO ) 00198 RETURN 00199 END IF 00200 * 00201 * Quick return if (N.EQ.0). 00202 * 00203 IF( N.EQ.0 ) 00204 + RETURN 00205 * 00206 * Skip single precision iterative refinement if a priori slower 00207 * than double precision factorization. 00208 * 00209 IF( .NOT.DOITREF ) THEN 00210 ITER = -1 00211 GO TO 40 00212 END IF 00213 * 00214 * Compute some constants. 00215 * 00216 ANRM = ZLANHE( 'I', UPLO, N, A, LDA, RWORK ) 00217 EPS = DLAMCH( 'Epsilon' ) 00218 CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX 00219 * 00220 * Set the indices PTSA, PTSX for referencing SA and SX in SWORK. 00221 * 00222 PTSA = 1 00223 PTSX = PTSA + N*N 00224 * 00225 * Convert B from double precision to single precision and store the 00226 * result in SX. 00227 * 00228 CALL ZLAG2C( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO ) 00229 * 00230 IF( INFO.NE.0 ) THEN 00231 ITER = -2 00232 GO TO 40 00233 END IF 00234 * 00235 * Convert A from double precision to single precision and store the 00236 * result in SA. 00237 * 00238 CALL ZLAT2C( UPLO, N, A, LDA, SWORK( PTSA ), N, INFO ) 00239 * 00240 IF( INFO.NE.0 ) THEN 00241 ITER = -2 00242 GO TO 40 00243 END IF 00244 * 00245 * Compute the Cholesky factorization of SA. 00246 * 00247 CALL CPOTRF( UPLO, N, SWORK( PTSA ), N, INFO ) 00248 * 00249 IF( INFO.NE.0 ) THEN 00250 ITER = -3 00251 GO TO 40 00252 END IF 00253 * 00254 * Solve the system SA*SX = SB. 00255 * 00256 CALL CPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N, 00257 + INFO ) 00258 * 00259 * Convert SX back to COMPLEX*16 00260 * 00261 CALL CLAG2Z( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO ) 00262 * 00263 * Compute R = B - AX (R is WORK). 00264 * 00265 CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N ) 00266 * 00267 CALL ZHEMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE, 00268 + WORK, N ) 00269 * 00270 * Check whether the NRHS normwise backward errors satisfy the 00271 * stopping criterion. If yes, set ITER=0 and return. 00272 * 00273 DO I = 1, NRHS 00274 XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) ) 00275 RNRM = CABS1( WORK( IZAMAX( N, WORK( 1, I ), 1 ), I ) ) 00276 IF( RNRM.GT.XNRM*CTE ) 00277 + GO TO 10 00278 END DO 00279 * 00280 * If we are here, the NRHS normwise backward errors satisfy the 00281 * stopping criterion. We are good to exit. 00282 * 00283 ITER = 0 00284 RETURN 00285 * 00286 10 CONTINUE 00287 * 00288 DO 30 IITER = 1, ITERMAX 00289 * 00290 * Convert R (in WORK) from double precision to single precision 00291 * and store the result in SX. 00292 * 00293 CALL ZLAG2C( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO ) 00294 * 00295 IF( INFO.NE.0 ) THEN 00296 ITER = -2 00297 GO TO 40 00298 END IF 00299 * 00300 * Solve the system SA*SX = SR. 00301 * 00302 CALL CPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N, 00303 + INFO ) 00304 * 00305 * Convert SX back to double precision and update the current 00306 * iterate. 00307 * 00308 CALL CLAG2Z( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO ) 00309 * 00310 DO I = 1, NRHS 00311 CALL ZAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 ) 00312 END DO 00313 * 00314 * Compute R = B - AX (R is WORK). 00315 * 00316 CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N ) 00317 * 00318 CALL ZHEMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE, 00319 + WORK, N ) 00320 * 00321 * Check whether the NRHS normwise backward errors satisfy the 00322 * stopping criterion. If yes, set ITER=IITER>0 and return. 00323 * 00324 DO I = 1, NRHS 00325 XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) ) 00326 RNRM = CABS1( WORK( IZAMAX( N, WORK( 1, I ), 1 ), I ) ) 00327 IF( RNRM.GT.XNRM*CTE ) 00328 + GO TO 20 00329 END DO 00330 * 00331 * If we are here, the NRHS normwise backward errors satisfy the 00332 * stopping criterion, we are good to exit. 00333 * 00334 ITER = IITER 00335 * 00336 RETURN 00337 * 00338 20 CONTINUE 00339 * 00340 30 CONTINUE 00341 * 00342 * If we are at this place of the code, this is because we have 00343 * performed ITER=ITERMAX iterations and never satisified the 00344 * stopping criterion, set up the ITER flag accordingly and follow 00345 * up on double precision routine. 00346 * 00347 ITER = -ITERMAX - 1 00348 * 00349 40 CONTINUE 00350 * 00351 * Single-precision iterative refinement failed to converge to a 00352 * satisfactory solution, so we resort to double precision. 00353 * 00354 CALL ZPOTRF( UPLO, N, A, LDA, INFO ) 00355 * 00356 IF( INFO.NE.0 ) 00357 + RETURN 00358 * 00359 CALL ZLACPY( 'All', N, NRHS, B, LDB, X, LDX ) 00360 CALL ZPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO ) 00361 * 00362 RETURN 00363 * 00364 * End of ZCPOSV. 00365 * 00366 END
1.7.2