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