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