SRC/zcposv.f

   1 *> \brief <b> ZCPOSV computes the solution to system of linear equations A * X = B for PO matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZCPOSV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zcposv.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zcposv.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zcposv.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZCPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
  22 *                          SWORK, RWORK, ITER, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO
  26 *       INTEGER            INFO, ITER, LDA, LDB, LDX, N, NRHS
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   RWORK( * )
  30 *       COMPLEX            SWORK( * )
  31 *       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( N, * ),
  32 *      $                   X( LDX, * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> ZCPOSV computes the solution to a complex system of linear equations
  42 *>    A * X = B,
  43 *> where A is an N-by-N Hermitian positive definite matrix and X and B
  44 *> are N-by-NRHS matrices.
  45 *>
  46 *> ZCPOSV first attempts to factorize the matrix in COMPLEX and use this
  47 *> factorization within an iterative refinement procedure to produce a
  48 *> solution with COMPLEX*16 normwise backward error quality (see below).
  49 *> If the approach fails the method switches to a COMPLEX*16
  50 *> factorization and solve.
  51 *>
  52 *> The iterative refinement is not going to be a winning strategy if
  53 *> the ratio COMPLEX performance over COMPLEX*16 performance is too
  54 *> small. A reasonable strategy should take the number of right-hand
  55 *> sides and the size of the matrix into account. This might be done
  56 *> with a call to ILAENV in the future. Up to now, we always try
  57 *> iterative refinement.
  58 *>
  59 *> The iterative refinement process is stopped if
  60 *>     ITER > ITERMAX
  61 *> or for all the RHS we have:
  62 *>     RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
  63 *> where
  64 *>     o ITER is the number of the current iteration in the iterative
  65 *>       refinement process
  66 *>     o RNRM is the infinity-norm of the residual
  67 *>     o XNRM is the infinity-norm of the solution
  68 *>     o ANRM is the infinity-operator-norm of the matrix A
  69 *>     o EPS is the machine epsilon returned by DLAMCH('Epsilon')
  70 *> The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
  71 *> respectively.
  72 *> \endverbatim
  73 *
  74 *  Arguments:
  75 *  ==========
  76 *
  77 *> \param[in] UPLO
  78 *> \verbatim
  79 *>          UPLO is CHARACTER*1
  80 *>          = 'U':  Upper triangle of A is stored;
  81 *>          = 'L':  Lower triangle of A is stored.
  82 *> \endverbatim
  83 *>
  84 *> \param[in] N
  85 *> \verbatim
  86 *>          N is INTEGER
  87 *>          The number of linear equations, i.e., the order of the
  88 *>          matrix A.  N >= 0.
  89 *> \endverbatim
  90 *>
  91 *> \param[in] NRHS
  92 *> \verbatim
  93 *>          NRHS is INTEGER
  94 *>          The number of right hand sides, i.e., the number of columns
  95 *>          of the matrix B.  NRHS >= 0.
  96 *> \endverbatim
  97 *>
  98 *> \param[in,out] A
  99 *> \verbatim
 100 *>          A is COMPLEX*16 array,
 101 *>          dimension (LDA,N)
 102 *>          On entry, the Hermitian matrix A. If UPLO = 'U', the leading
 103 *>          N-by-N upper triangular part of A contains the upper
 104 *>          triangular part of the matrix A, and the strictly lower
 105 *>          triangular part of A is not referenced.  If UPLO = 'L', the
 106 *>          leading N-by-N lower triangular part of A contains the lower
 107 *>          triangular part of the matrix A, and the strictly upper
 108 *>          triangular part of A is not referenced.
 109 *>
 110 *>          Note that the imaginary parts of the diagonal
 111 *>          elements need not be set and are assumed to be zero.
 112 *>
 113 *>          On exit, if iterative refinement has been successfully used
 114 *>          (INFO.EQ.0 and ITER.GE.0, see description below), then A is
 115 *>          unchanged, if double precision factorization has been used
 116 *>          (INFO.EQ.0 and ITER.LT.0, see description below), then the
 117 *>          array A contains the factor U or L from the Cholesky
 118 *>          factorization A = U**H*U or A = L*L**H.
 119 *> \endverbatim
 120 *>
 121 *> \param[in] LDA
 122 *> \verbatim
 123 *>          LDA is INTEGER
 124 *>          The leading dimension of the array A.  LDA >= max(1,N).
 125 *> \endverbatim
 126 *>
 127 *> \param[in] B
 128 *> \verbatim
 129 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
 130 *>          The N-by-NRHS right hand side matrix B.
 131 *> \endverbatim
 132 *>
 133 *> \param[in] LDB
 134 *> \verbatim
 135 *>          LDB is INTEGER
 136 *>          The leading dimension of the array B.  LDB >= max(1,N).
 137 *> \endverbatim
 138 *>
 139 *> \param[out] X
 140 *> \verbatim
 141 *>          X is COMPLEX*16 array, dimension (LDX,NRHS)
 142 *>          If INFO = 0, the N-by-NRHS solution matrix X.
 143 *> \endverbatim
 144 *>
 145 *> \param[in] LDX
 146 *> \verbatim
 147 *>          LDX is INTEGER
 148 *>          The leading dimension of the array X.  LDX >= max(1,N).
 149 *> \endverbatim
 150 *>
 151 *> \param[out] WORK
 152 *> \verbatim
 153 *>          WORK is COMPLEX*16 array, dimension (N*NRHS)
 154 *>          This array is used to hold the residual vectors.
 155 *> \endverbatim
 156 *>
 157 *> \param[out] SWORK
 158 *> \verbatim
 159 *>          SWORK is COMPLEX array, dimension (N*(N+NRHS))
 160 *>          This array is used to use the single precision matrix and the
 161 *>          right-hand sides or solutions in single precision.
 162 *> \endverbatim
 163 *>
 164 *> \param[out] RWORK
 165 *> \verbatim
 166 *>          RWORK is DOUBLE PRECISION array, dimension (N)
 167 *> \endverbatim
 168 *>
 169 *> \param[out] ITER
 170 *> \verbatim
 171 *>          ITER is INTEGER
 172 *>          < 0: iterative refinement has failed, COMPLEX*16
 173 *>               factorization has been performed
 174 *>               -1 : the routine fell back to full precision for
 175 *>                    implementation- or machine-specific reasons
 176 *>               -2 : narrowing the precision induced an overflow,
 177 *>                    the routine fell back to full precision
 178 *>               -3 : failure of CPOTRF
 179 *>               -31: stop the iterative refinement after the 30th
 180 *>                    iterations
 181 *>          > 0: iterative refinement has been successfully used.
 182 *>               Returns the number of iterations
 183 *> \endverbatim
 184 *>
 185 *> \param[out] INFO
 186 *> \verbatim
 187 *>          INFO is INTEGER
 188 *>          = 0:  successful exit
 189 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 190 *>          > 0:  if INFO = i, the leading minor of order i of
 191 *>                (COMPLEX*16) A is not positive definite, so the
 192 *>                factorization could not be completed, and the solution
 193 *>                has not been computed.
 194 *> \endverbatim
 195 *
 196 *  Authors:
 197 *  ========
 198 *
 199 *> \author Univ. of Tennessee
 200 *> \author Univ. of California Berkeley
 201 *> \author Univ. of Colorado Denver
 202 *> \author NAG Ltd.
 203 *
 204 *> \date June 2016
 205 *
 206 *> \ingroup complex16POsolve
 207 *
 208 *  =====================================================================
 209       SUBROUTINE ZCPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
 210      $                   SWORK, RWORK, ITER, INFO )
 211 *
 212 *  -- LAPACK driver routine (version 3.6.1) --
 213 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 214 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 215 *     June 2016
 216 *
 217 *     .. Scalar Arguments ..
 218       CHARACTER          UPLO
 219       INTEGER            INFO, ITER, LDA, LDB, LDX, N, NRHS
 220 *     ..
 221 *     .. Array Arguments ..
 222       DOUBLE PRECISION   RWORK( * )
 223       COMPLEX            SWORK( * )
 224       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( N, * ),
 225      $                   X( LDX, * )
 226 *     ..
 227 *
 228 *  =====================================================================
 229 *
 230 *     .. Parameters ..
 231       LOGICAL            DOITREF
 232       PARAMETER          ( DOITREF = .TRUE. )
 233 *
 234       INTEGER            ITERMAX
 235       PARAMETER          ( ITERMAX = 30 )
 236 *
 237       DOUBLE PRECISION   BWDMAX
 238       PARAMETER          ( BWDMAX = 1.0E+00 )
 239 *
 240       COMPLEX*16         NEGONE, ONE
 241       PARAMETER          ( NEGONE = ( -1.0D+00, 0.0D+00 ),
 242      $                   ONE = ( 1.0D+00, 0.0D+00 ) )
 243 *
 244 *     .. Local Scalars ..
 245       INTEGER            I, IITER, PTSA, PTSX
 246       DOUBLE PRECISION   ANRM, CTE, EPS, RNRM, XNRM
 247       COMPLEX*16         ZDUM
 248 *
 249 *     .. External Subroutines ..
 250       EXTERNAL           ZAXPY, ZHEMM, ZLACPY, ZLAT2C, ZLAG2C, CLAG2Z,
 251      $                   CPOTRF, CPOTRS, XERBLA
 252 *     ..
 253 *     .. External Functions ..
 254       INTEGER            IZAMAX
 255       DOUBLE PRECISION   DLAMCH, ZLANHE
 256       LOGICAL            LSAME
 257       EXTERNAL           IZAMAX, DLAMCH, ZLANHE, LSAME
 258 *     ..
 259 *     .. Intrinsic Functions ..
 260       INTRINSIC          ABS, DBLE, MAX, SQRT
 261 *     .. Statement Functions ..
 262       DOUBLE PRECISION   CABS1
 263 *     ..
 264 *     .. Statement Function definitions ..
 265       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
 266 *     ..
 267 *     .. Executable Statements ..
 268 *
 269       INFO = 0
 270       ITER = 0
 271 *
 272 *     Test the input parameters.
 273 *
 274       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 275          INFO = -1
 276       ELSE IF( N.LT.0 ) THEN
 277          INFO = -2
 278       ELSE IF( NRHS.LT.0 ) THEN
 279          INFO = -3
 280       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 281          INFO = -5
 282       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 283          INFO = -7
 284       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
 285          INFO = -9
 286       END IF
 287       IF( INFO.NE.0 ) THEN
 288          CALL XERBLA( 'ZCPOSV', -INFO )
 289          RETURN
 290       END IF
 291 *
 292 *     Quick return if (N.EQ.0).
 293 *
 294       IF( N.EQ.0 )
 295      $   RETURN
 296 *
 297 *     Skip single precision iterative refinement if a priori slower
 298 *     than double precision factorization.
 299 *
 300       IF( .NOT.DOITREF ) THEN
 301          ITER = -1
 302          GO TO 40
 303       END IF
 304 *
 305 *     Compute some constants.
 306 *
 307       ANRM = ZLANHE( 'I', UPLO, N, A, LDA, RWORK )
 308       EPS = DLAMCH( 'Epsilon' )
 309       CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX
 310 *
 311 *     Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
 312 *
 313       PTSA = 1
 314       PTSX = PTSA + N*N
 315 *
 316 *     Convert B from double precision to single precision and store the
 317 *     result in SX.
 318 *
 319       CALL ZLAG2C( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO )
 320 *
 321       IF( INFO.NE.0 ) THEN
 322          ITER = -2
 323          GO TO 40
 324       END IF
 325 *
 326 *     Convert A from double precision to single precision and store the
 327 *     result in SA.
 328 *
 329       CALL ZLAT2C( UPLO, N, A, LDA, SWORK( PTSA ), N, INFO )
 330 *
 331       IF( INFO.NE.0 ) THEN
 332          ITER = -2
 333          GO TO 40
 334       END IF
 335 *
 336 *     Compute the Cholesky factorization of SA.
 337 *
 338       CALL CPOTRF( UPLO, N, SWORK( PTSA ), N, INFO )
 339 *
 340       IF( INFO.NE.0 ) THEN
 341          ITER = -3
 342          GO TO 40
 343       END IF
 344 *
 345 *     Solve the system SA*SX = SB.
 346 *
 347       CALL CPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
 348      $             INFO )
 349 *
 350 *     Convert SX back to COMPLEX*16
 351 *
 352       CALL CLAG2Z( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
 353 *
 354 *     Compute R = B - AX (R is WORK).
 355 *
 356       CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N )
 357 *
 358       CALL ZHEMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
 359      $            WORK, N )
 360 *
 361 *     Check whether the NRHS normwise backward errors satisfy the
 362 *     stopping criterion. If yes, set ITER=0 and return.
 363 *
 364       DO I = 1, NRHS
 365          XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) )
 366          RNRM = CABS1( WORK( IZAMAX( N, WORK( 1, I ), 1 ), I ) )
 367          IF( RNRM.GT.XNRM*CTE )
 368      $      GO TO 10
 369       END DO
 370 *
 371 *     If we are here, the NRHS normwise backward errors satisfy the
 372 *     stopping criterion. We are good to exit.
 373 *
 374       ITER = 0
 375       RETURN
 376 *
 377    10 CONTINUE
 378 *
 379       DO 30 IITER = 1, ITERMAX
 380 *
 381 *        Convert R (in WORK) from double precision to single precision
 382 *        and store the result in SX.
 383 *
 384          CALL ZLAG2C( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO )
 385 *
 386          IF( INFO.NE.0 ) THEN
 387             ITER = -2
 388             GO TO 40
 389          END IF
 390 *
 391 *        Solve the system SA*SX = SR.
 392 *
 393          CALL CPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
 394      $                INFO )
 395 *
 396 *        Convert SX back to double precision and update the current
 397 *        iterate.
 398 *
 399          CALL CLAG2Z( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
 400 *
 401          DO I = 1, NRHS
 402             CALL ZAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 )
 403          END DO
 404 *
 405 *        Compute R = B - AX (R is WORK).
 406 *
 407          CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N )
 408 *
 409          CALL ZHEMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
 410      $               WORK, N )
 411 *
 412 *        Check whether the NRHS normwise backward errors satisfy the
 413 *        stopping criterion. If yes, set ITER=IITER>0 and return.
 414 *
 415          DO I = 1, NRHS
 416             XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) )
 417             RNRM = CABS1( WORK( IZAMAX( N, WORK( 1, I ), 1 ), I ) )
 418             IF( RNRM.GT.XNRM*CTE )
 419      $         GO TO 20
 420          END DO
 421 *
 422 *        If we are here, the NRHS normwise backward errors satisfy the
 423 *        stopping criterion, we are good to exit.
 424 *
 425          ITER = IITER
 426 *
 427          RETURN
 428 *
 429    20    CONTINUE
 430 *
 431    30 CONTINUE
 432 *
 433 *     If we are at this place of the code, this is because we have
 434 *     performed ITER=ITERMAX iterations and never satisified the
 435 *     stopping criterion, set up the ITER flag accordingly and follow
 436 *     up on double precision routine.
 437 *
 438       ITER = -ITERMAX - 1
 439 *
 440    40 CONTINUE
 441 *
 442 *     Single-precision iterative refinement failed to converge to a
 443 *     satisfactory solution, so we resort to double precision.
 444 *
 445       CALL ZPOTRF( UPLO, N, A, LDA, INFO )
 446 *
 447       IF( INFO.NE.0 )
 448      $   RETURN
 449 *
 450       CALL ZLACPY( 'All', N, NRHS, B, LDB, X, LDX )
 451       CALL ZPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO )
 452 *
 453       RETURN
 454 *
 455 *     End of ZCPOSV.
 456 *
 457       END