LAPACK 3.3.1
Linear Algebra PACKage

zcposv.f

Go to the documentation of this file.
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.1)                 --
00005 *
00006 *  -- April 2011                                                      --
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
 All Files Functions