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