SRC/dsposv.f

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