LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zcposv.f
Go to the documentation of this file.
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
 All Files Functions