SRC/dptrfs.f

   1 *> \brief \b DPTRFS
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DPTRFS + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dptrfs.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dptrfs.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dptrfs.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DPTRFS( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR,
  22 *                          BERR, WORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, LDB, LDX, N, NRHS
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   B( LDB, * ), BERR( * ), D( * ), DF( * ),
  29 *      $                   E( * ), EF( * ), FERR( * ), WORK( * ),
  30 *      $                   X( LDX, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> DPTRFS improves the computed solution to a system of linear
  40 *> equations when the coefficient matrix is symmetric positive definite
  41 *> and tridiagonal, and provides error bounds and backward error
  42 *> estimates for the solution.
  43 *> \endverbatim
  44 *
  45 *  Arguments:
  46 *  ==========
  47 *
  48 *> \param[in] N
  49 *> \verbatim
  50 *>          N is INTEGER
  51 *>          The order of the matrix A.  N >= 0.
  52 *> \endverbatim
  53 *>
  54 *> \param[in] NRHS
  55 *> \verbatim
  56 *>          NRHS is INTEGER
  57 *>          The number of right hand sides, i.e., the number of columns
  58 *>          of the matrix B.  NRHS >= 0.
  59 *> \endverbatim
  60 *>
  61 *> \param[in] D
  62 *> \verbatim
  63 *>          D is DOUBLE PRECISION array, dimension (N)
  64 *>          The n diagonal elements of the tridiagonal matrix A.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] E
  68 *> \verbatim
  69 *>          E is DOUBLE PRECISION array, dimension (N-1)
  70 *>          The (n-1) subdiagonal elements of the tridiagonal matrix A.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] DF
  74 *> \verbatim
  75 *>          DF is DOUBLE PRECISION array, dimension (N)
  76 *>          The n diagonal elements of the diagonal matrix D from the
  77 *>          factorization computed by DPTTRF.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] EF
  81 *> \verbatim
  82 *>          EF is DOUBLE PRECISION array, dimension (N-1)
  83 *>          The (n-1) subdiagonal elements of the unit bidiagonal factor
  84 *>          L from the factorization computed by DPTTRF.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] B
  88 *> \verbatim
  89 *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
  90 *>          The right hand side matrix B.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] LDB
  94 *> \verbatim
  95 *>          LDB is INTEGER
  96 *>          The leading dimension of the array B.  LDB >= max(1,N).
  97 *> \endverbatim
  98 *>
  99 *> \param[in,out] X
 100 *> \verbatim
 101 *>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
 102 *>          On entry, the solution matrix X, as computed by DPTTRS.
 103 *>          On exit, the improved solution matrix X.
 104 *> \endverbatim
 105 *>
 106 *> \param[in] LDX
 107 *> \verbatim
 108 *>          LDX is INTEGER
 109 *>          The leading dimension of the array X.  LDX >= max(1,N).
 110 *> \endverbatim
 111 *>
 112 *> \param[out] FERR
 113 *> \verbatim
 114 *>          FERR is DOUBLE PRECISION array, dimension (NRHS)
 115 *>          The forward error bound for each solution vector
 116 *>          X(j) (the j-th column of the solution matrix X).
 117 *>          If XTRUE is the true solution corresponding to X(j), FERR(j)
 118 *>          is an estimated upper bound for the magnitude of the largest
 119 *>          element in (X(j) - XTRUE) divided by the magnitude of the
 120 *>          largest element in X(j).
 121 *> \endverbatim
 122 *>
 123 *> \param[out] BERR
 124 *> \verbatim
 125 *>          BERR is DOUBLE PRECISION array, dimension (NRHS)
 126 *>          The componentwise relative backward error of each solution
 127 *>          vector X(j) (i.e., the smallest relative change in
 128 *>          any element of A or B that makes X(j) an exact solution).
 129 *> \endverbatim
 130 *>
 131 *> \param[out] WORK
 132 *> \verbatim
 133 *>          WORK is DOUBLE PRECISION array, dimension (2*N)
 134 *> \endverbatim
 135 *>
 136 *> \param[out] INFO
 137 *> \verbatim
 138 *>          INFO is INTEGER
 139 *>          = 0:  successful exit
 140 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 141 *> \endverbatim
 142 *
 143 *> \par Internal Parameters:
 144 *  =========================
 145 *>
 146 *> \verbatim
 147 *>  ITMAX is the maximum number of steps of iterative refinement.
 148 *> \endverbatim
 149 *
 150 *  Authors:
 151 *  ========
 152 *
 153 *> \author Univ. of Tennessee
 154 *> \author Univ. of California Berkeley
 155 *> \author Univ. of Colorado Denver
 156 *> \author NAG Ltd.
 157 *
 158 *> \date September 2012
 159 *
 160 *> \ingroup doublePTcomputational
 161 *
 162 *  =====================================================================
 163       SUBROUTINE DPTRFS( N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR,
 164      $                   BERR, WORK, INFO )
 165 *
 166 *  -- LAPACK computational routine (version 3.4.2) --
 167 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 168 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 169 *     September 2012
 170 *
 171 *     .. Scalar Arguments ..
 172       INTEGER            INFO, LDB, LDX, N, NRHS
 173 *     ..
 174 *     .. Array Arguments ..
 175       DOUBLE PRECISION   B( LDB, * ), BERR( * ), D( * ), DF( * ),
 176      $                   E( * ), EF( * ), FERR( * ), WORK( * ),
 177      $                   X( LDX, * )
 178 *     ..
 179 *
 180 *  =====================================================================
 181 *
 182 *     .. Parameters ..
 183       INTEGER            ITMAX
 184       PARAMETER          ( ITMAX = 5 )
 185       DOUBLE PRECISION   ZERO
 186       PARAMETER          ( ZERO = 0.0D+0 )
 187       DOUBLE PRECISION   ONE
 188       PARAMETER          ( ONE = 1.0D+0 )
 189       DOUBLE PRECISION   TWO
 190       PARAMETER          ( TWO = 2.0D+0 )
 191       DOUBLE PRECISION   THREE
 192       PARAMETER          ( THREE = 3.0D+0 )
 193 *     ..
 194 *     .. Local Scalars ..
 195       INTEGER            COUNT, I, IX, J, NZ
 196       DOUBLE PRECISION   BI, CX, DX, EPS, EX, LSTRES, S, SAFE1, SAFE2,
 197      $                   SAFMIN
 198 *     ..
 199 *     .. External Subroutines ..
 200       EXTERNAL           DAXPY, DPTTRS, XERBLA
 201 *     ..
 202 *     .. Intrinsic Functions ..
 203       INTRINSIC          ABS, MAX
 204 *     ..
 205 *     .. External Functions ..
 206       INTEGER            IDAMAX
 207       DOUBLE PRECISION   DLAMCH
 208       EXTERNAL           IDAMAX, DLAMCH
 209 *     ..
 210 *     .. Executable Statements ..
 211 *
 212 *     Test the input parameters.
 213 *
 214       INFO = 0
 215       IF( N.LT.0 ) THEN
 216          INFO = -1
 217       ELSE IF( NRHS.LT.0 ) THEN
 218          INFO = -2
 219       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 220          INFO = -8
 221       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
 222          INFO = -10
 223       END IF
 224       IF( INFO.NE.0 ) THEN
 225          CALL XERBLA( 'DPTRFS', -INFO )
 226          RETURN
 227       END IF
 228 *
 229 *     Quick return if possible
 230 *
 231       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
 232          DO 10 J = 1, NRHS
 233             FERR( J ) = ZERO
 234             BERR( J ) = ZERO
 235    10    CONTINUE
 236          RETURN
 237       END IF
 238 *
 239 *     NZ = maximum number of nonzero elements in each row of A, plus 1
 240 *
 241       NZ = 4
 242       EPS = DLAMCH( 'Epsilon' )
 243       SAFMIN = DLAMCH( 'Safe minimum' )
 244       SAFE1 = NZ*SAFMIN
 245       SAFE2 = SAFE1 / EPS
 246 *
 247 *     Do for each right hand side
 248 *
 249       DO 90 J = 1, NRHS
 250 *
 251          COUNT = 1
 252          LSTRES = THREE
 253    20    CONTINUE
 254 *
 255 *        Loop until stopping criterion is satisfied.
 256 *
 257 *        Compute residual R = B - A * X.  Also compute
 258 *        abs(A)*abs(x) + abs(b) for use in the backward error bound.
 259 *
 260          IF( N.EQ.1 ) THEN
 261             BI = B( 1, J )
 262             DX = D( 1 )*X( 1, J )
 263             WORK( N+1 ) = BI - DX
 264             WORK( 1 ) = ABS( BI ) + ABS( DX )
 265          ELSE
 266             BI = B( 1, J )
 267             DX = D( 1 )*X( 1, J )
 268             EX = E( 1 )*X( 2, J )
 269             WORK( N+1 ) = BI - DX - EX
 270             WORK( 1 ) = ABS( BI ) + ABS( DX ) + ABS( EX )
 271             DO 30 I = 2, N - 1
 272                BI = B( I, J )
 273                CX = E( I-1 )*X( I-1, J )
 274                DX = D( I )*X( I, J )
 275                EX = E( I )*X( I+1, J )
 276                WORK( N+I ) = BI - CX - DX - EX
 277                WORK( I ) = ABS( BI ) + ABS( CX ) + ABS( DX ) + ABS( EX )
 278    30       CONTINUE
 279             BI = B( N, J )
 280             CX = E( N-1 )*X( N-1, J )
 281             DX = D( N )*X( N, J )
 282             WORK( N+N ) = BI - CX - DX
 283             WORK( N ) = ABS( BI ) + ABS( CX ) + ABS( DX )
 284          END IF
 285 *
 286 *        Compute componentwise relative backward error from formula
 287 *
 288 *        max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
 289 *
 290 *        where abs(Z) is the componentwise absolute value of the matrix
 291 *        or vector Z.  If the i-th component of the denominator is less
 292 *        than SAFE2, then SAFE1 is added to the i-th components of the
 293 *        numerator and denominator before dividing.
 294 *
 295          S = ZERO
 296          DO 40 I = 1, N
 297             IF( WORK( I ).GT.SAFE2 ) THEN
 298                S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
 299             ELSE
 300                S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
 301      $             ( WORK( I )+SAFE1 ) )
 302             END IF
 303    40    CONTINUE
 304          BERR( J ) = S
 305 *
 306 *        Test stopping criterion. Continue iterating if
 307 *           1) The residual BERR(J) is larger than machine epsilon, and
 308 *           2) BERR(J) decreased by at least a factor of 2 during the
 309 *              last iteration, and
 310 *           3) At most ITMAX iterations tried.
 311 *
 312          IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
 313      $       COUNT.LE.ITMAX ) THEN
 314 *
 315 *           Update solution and try again.
 316 *
 317             CALL DPTTRS( N, 1, DF, EF, WORK( N+1 ), N, INFO )
 318             CALL DAXPY( N, ONE, WORK( N+1 ), 1, X( 1, J ), 1 )
 319             LSTRES = BERR( J )
 320             COUNT = COUNT + 1
 321             GO TO 20
 322          END IF
 323 *
 324 *        Bound error from formula
 325 *
 326 *        norm(X - XTRUE) / norm(X) .le. FERR =
 327 *        norm( abs(inv(A))*
 328 *           ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
 329 *
 330 *        where
 331 *          norm(Z) is the magnitude of the largest component of Z
 332 *          inv(A) is the inverse of A
 333 *          abs(Z) is the componentwise absolute value of the matrix or
 334 *             vector Z
 335 *          NZ is the maximum number of nonzeros in any row of A, plus 1
 336 *          EPS is machine epsilon
 337 *
 338 *        The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
 339 *        is incremented by SAFE1 if the i-th component of
 340 *        abs(A)*abs(X) + abs(B) is less than SAFE2.
 341 *
 342          DO 50 I = 1, N
 343             IF( WORK( I ).GT.SAFE2 ) THEN
 344                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
 345             ELSE
 346                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
 347             END IF
 348    50    CONTINUE
 349          IX = IDAMAX( N, WORK, 1 )
 350          FERR( J ) = WORK( IX )
 351 *
 352 *        Estimate the norm of inv(A).
 353 *
 354 *        Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
 355 *
 356 *           m(i,j) =  abs(A(i,j)), i = j,
 357 *           m(i,j) = -abs(A(i,j)), i .ne. j,
 358 *
 359 *        and e = [ 1, 1, ..., 1 ]**T.  Note M(A) = M(L)*D*M(L)**T.
 360 *
 361 *        Solve M(L) * x = e.
 362 *
 363          WORK( 1 ) = ONE
 364          DO 60 I = 2, N
 365             WORK( I ) = ONE + WORK( I-1 )*ABS( EF( I-1 ) )
 366    60    CONTINUE
 367 *
 368 *        Solve D * M(L)**T * x = b.
 369 *
 370          WORK( N ) = WORK( N ) / DF( N )
 371          DO 70 I = N - 1, 1, -1
 372             WORK( I ) = WORK( I ) / DF( I ) + WORK( I+1 )*ABS( EF( I ) )
 373    70    CONTINUE
 374 *
 375 *        Compute norm(inv(A)) = max(x(i)), 1<=i<=n.
 376 *
 377          IX = IDAMAX( N, WORK, 1 )
 378          FERR( J ) = FERR( J )*ABS( WORK( IX ) )
 379 *
 380 *        Normalize error.
 381 *
 382          LSTRES = ZERO
 383          DO 80 I = 1, N
 384             LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
 385    80    CONTINUE
 386          IF( LSTRES.NE.ZERO )
 387      $      FERR( J ) = FERR( J ) / LSTRES
 388 *
 389    90 CONTINUE
 390 *
 391       RETURN
 392 *
 393 *     End of DPTRFS
 394 *
 395       END