SRC/dpbrfs.f

   1       SUBROUTINE DPBRFS( UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, B,
   2      $                   LDB, X, LDX, FERR, BERR, WORK, IWORK, INFO )
   3 *
   4 *  -- LAPACK routine (version 3.2) --
   5 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
   6 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
   7 *     November 2006
   8 *
   9 *     Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
  10 *
  11 *     .. Scalar Arguments ..
  12       CHARACTER          UPLO
  13       INTEGER            INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS
  14 *     ..
  15 *     .. Array Arguments ..
  16       INTEGER            IWORK( * )
  17       DOUBLE PRECISION   AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
  18      $                   BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
  19 *     ..
  20 *
  21 *  Purpose
  22 *  =======
  23 *
  24 *  DPBRFS improves the computed solution to a system of linear
  25 *  equations when the coefficient matrix is symmetric positive definite
  26 *  and banded, and provides error bounds and backward error estimates
  27 *  for the solution.
  28 *
  29 *  Arguments
  30 *  =========
  31 *
  32 *  UPLO    (input) CHARACTER*1
  33 *          = 'U':  Upper triangle of A is stored;
  34 *          = 'L':  Lower triangle of A is stored.
  35 *
  36 *  N       (input) INTEGER
  37 *          The order of the matrix A.  N >= 0.
  38 *
  39 *  KD      (input) INTEGER
  40 *          The number of superdiagonals of the matrix A if UPLO = 'U',
  41 *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
  42 *
  43 *  NRHS    (input) INTEGER
  44 *          The number of right hand sides, i.e., the number of columns
  45 *          of the matrices B and X.  NRHS >= 0.
  46 *
  47 *  AB      (input) DOUBLE PRECISION array, dimension (LDAB,N)
  48 *          The upper or lower triangle of the symmetric band matrix A,
  49 *          stored in the first KD+1 rows of the array.  The j-th column
  50 *          of A is stored in the j-th column of the array AB as follows:
  51 *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
  52 *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
  53 *
  54 *  LDAB    (input) INTEGER
  55 *          The leading dimension of the array AB.  LDAB >= KD+1.
  56 *
  57 *  AFB     (input) DOUBLE PRECISION array, dimension (LDAFB,N)
  58 *          The triangular factor U or L from the Cholesky factorization
  59 *          A = U**T*U or A = L*L**T of the band matrix A as computed by
  60 *          DPBTRF, in the same storage format as A (see AB).
  61 *
  62 *  LDAFB   (input) INTEGER
  63 *          The leading dimension of the array AFB.  LDAFB >= KD+1.
  64 *
  65 *  B       (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
  66 *          The right hand side matrix B.
  67 *
  68 *  LDB     (input) INTEGER
  69 *          The leading dimension of the array B.  LDB >= max(1,N).
  70 *
  71 *  X       (input/output) DOUBLE PRECISION array, dimension (LDX,NRHS)
  72 *          On entry, the solution matrix X, as computed by DPBTRS.
  73 *          On exit, the improved solution matrix X.
  74 *
  75 *  LDX     (input) INTEGER
  76 *          The leading dimension of the array X.  LDX >= max(1,N).
  77 *
  78 *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
  79 *          The estimated forward error bound for each solution vector
  80 *          X(j) (the j-th column of the solution matrix X).
  81 *          If XTRUE is the true solution corresponding to X(j), FERR(j)
  82 *          is an estimated upper bound for the magnitude of the largest
  83 *          element in (X(j) - XTRUE) divided by the magnitude of the
  84 *          largest element in X(j).  The estimate is as reliable as
  85 *          the estimate for RCOND, and is almost always a slight
  86 *          overestimate of the true error.
  87 *
  88 *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
  89 *          The componentwise relative backward error of each solution
  90 *          vector X(j) (i.e., the smallest relative change in
  91 *          any element of A or B that makes X(j) an exact solution).
  92 *
  93 *  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
  94 *
  95 *  IWORK   (workspace) INTEGER array, dimension (N)
  96 *
  97 *  INFO    (output) INTEGER
  98 *          = 0:  successful exit
  99 *          < 0:  if INFO = -i, the i-th argument had an illegal value
 100 *
 101 *  Internal Parameters
 102 *  ===================
 103 *
 104 *  ITMAX is the maximum number of steps of iterative refinement.
 105 *
 106 *  =====================================================================
 107 *
 108 *     .. Parameters ..
 109       INTEGER            ITMAX
 110       PARAMETER          ( ITMAX = 5 )
 111       DOUBLE PRECISION   ZERO
 112       PARAMETER          ( ZERO = 0.0D+0 )
 113       DOUBLE PRECISION   ONE
 114       PARAMETER          ( ONE = 1.0D+0 )
 115       DOUBLE PRECISION   TWO
 116       PARAMETER          ( TWO = 2.0D+0 )
 117       DOUBLE PRECISION   THREE
 118       PARAMETER          ( THREE = 3.0D+0 )
 119 *     ..
 120 *     .. Local Scalars ..
 121       LOGICAL            UPPER
 122       INTEGER            COUNT, I, J, K, KASE, L, NZ
 123       DOUBLE PRECISION   EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
 124 *     ..
 125 *     .. Local Arrays ..
 126       INTEGER            ISAVE( 3 )
 127 *     ..
 128 *     .. External Subroutines ..
 129       EXTERNAL           DAXPY, DCOPY, DLACN2, DPBTRS, DSBMV, XERBLA
 130 *     ..
 131 *     .. Intrinsic Functions ..
 132       INTRINSIC          ABS, MAX, MIN
 133 *     ..
 134 *     .. External Functions ..
 135       LOGICAL            LSAME
 136       DOUBLE PRECISION   DLAMCH
 137       EXTERNAL           LSAME, DLAMCH
 138 *     ..
 139 *     .. Executable Statements ..
 140 *
 141 *     Test the input parameters.
 142 *
 143       INFO = 0
 144       UPPER = LSAME( UPLO, 'U' )
 145       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 146          INFO = -1
 147       ELSE IF( N.LT.0 ) THEN
 148          INFO = -2
 149       ELSE IF( KD.LT.0 ) THEN
 150          INFO = -3
 151       ELSE IF( NRHS.LT.0 ) THEN
 152          INFO = -4
 153       ELSE IF( LDAB.LT.KD+1 ) THEN
 154          INFO = -6
 155       ELSE IF( LDAFB.LT.KD+1 ) THEN
 156          INFO = -8
 157       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 158          INFO = -10
 159       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
 160          INFO = -12
 161       END IF
 162       IF( INFO.NE.0 ) THEN
 163          CALL XERBLA( 'DPBRFS', -INFO )
 164          RETURN
 165       END IF
 166 *
 167 *     Quick return if possible
 168 *
 169       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
 170          DO 10 J = 1, NRHS
 171             FERR( J ) = ZERO
 172             BERR( J ) = ZERO
 173    10    CONTINUE
 174          RETURN
 175       END IF
 176 *
 177 *     NZ = maximum number of nonzero elements in each row of A, plus 1
 178 *
 179       NZ = MIN( N+1, 2*KD+2 )
 180       EPS = DLAMCH( 'Epsilon' )
 181       SAFMIN = DLAMCH( 'Safe minimum' )
 182       SAFE1 = NZ*SAFMIN
 183       SAFE2 = SAFE1 / EPS
 184 *
 185 *     Do for each right hand side
 186 *
 187       DO 140 J = 1, NRHS
 188 *
 189          COUNT = 1
 190          LSTRES = THREE
 191    20    CONTINUE
 192 *
 193 *        Loop until stopping criterion is satisfied.
 194 *
 195 *        Compute residual R = B - A * X
 196 *
 197          CALL DCOPY( N, B( 1, J ), 1, WORK( N+1 ), 1 )
 198          CALL DSBMV( UPLO, N, KD, -ONE, AB, LDAB, X( 1, J ), 1, ONE,
 199      $               WORK( N+1 ), 1 )
 200 *
 201 *        Compute componentwise relative backward error from formula
 202 *
 203 *        max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
 204 *
 205 *        where abs(Z) is the componentwise absolute value of the matrix
 206 *        or vector Z.  If the i-th component of the denominator is less
 207 *        than SAFE2, then SAFE1 is added to the i-th components of the
 208 *        numerator and denominator before dividing.
 209 *
 210          DO 30 I = 1, N
 211             WORK( I ) = ABS( B( I, J ) )
 212    30    CONTINUE
 213 *
 214 *        Compute abs(A)*abs(X) + abs(B).
 215 *
 216          IF( UPPER ) THEN
 217             DO 50 K = 1, N
 218                S = ZERO
 219                XK = ABS( X( K, J ) )
 220                L = KD + 1 - K
 221                DO 40 I = MAX( 1, K-KD ), K - 1
 222                   WORK( I ) = WORK( I ) + ABS( AB( L+I, K ) )*XK
 223                   S = S + ABS( AB( L+I, K ) )*ABS( X( I, J ) )
 224    40          CONTINUE
 225                WORK( K ) = WORK( K ) + ABS( AB( KD+1, K ) )*XK + S
 226    50       CONTINUE
 227          ELSE
 228             DO 70 K = 1, N
 229                S = ZERO
 230                XK = ABS( X( K, J ) )
 231                WORK( K ) = WORK( K ) + ABS( AB( 1, K ) )*XK
 232                L = 1 - K
 233                DO 60 I = K + 1, MIN( N, K+KD )
 234                   WORK( I ) = WORK( I ) + ABS( AB( L+I, K ) )*XK
 235                   S = S + ABS( AB( L+I, K ) )*ABS( X( I, J ) )
 236    60          CONTINUE
 237                WORK( K ) = WORK( K ) + S
 238    70       CONTINUE
 239          END IF
 240          S = ZERO
 241          DO 80 I = 1, N
 242             IF( WORK( I ).GT.SAFE2 ) THEN
 243                S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
 244             ELSE
 245                S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
 246      $             ( WORK( I )+SAFE1 ) )
 247             END IF
 248    80    CONTINUE
 249          BERR( J ) = S
 250 *
 251 *        Test stopping criterion. Continue iterating if
 252 *           1) The residual BERR(J) is larger than machine epsilon, and
 253 *           2) BERR(J) decreased by at least a factor of 2 during the
 254 *              last iteration, and
 255 *           3) At most ITMAX iterations tried.
 256 *
 257          IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
 258      $       COUNT.LE.ITMAX ) THEN
 259 *
 260 *           Update solution and try again.
 261 *
 262             CALL DPBTRS( UPLO, N, KD, 1, AFB, LDAFB, WORK( N+1 ), N,
 263      $                   INFO )
 264             CALL DAXPY( N, ONE, WORK( N+1 ), 1, X( 1, J ), 1 )
 265             LSTRES = BERR( J )
 266             COUNT = COUNT + 1
 267             GO TO 20
 268          END IF
 269 *
 270 *        Bound error from formula
 271 *
 272 *        norm(X - XTRUE) / norm(X) .le. FERR =
 273 *        norm( abs(inv(A))*
 274 *           ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
 275 *
 276 *        where
 277 *          norm(Z) is the magnitude of the largest component of Z
 278 *          inv(A) is the inverse of A
 279 *          abs(Z) is the componentwise absolute value of the matrix or
 280 *             vector Z
 281 *          NZ is the maximum number of nonzeros in any row of A, plus 1
 282 *          EPS is machine epsilon
 283 *
 284 *        The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
 285 *        is incremented by SAFE1 if the i-th component of
 286 *        abs(A)*abs(X) + abs(B) is less than SAFE2.
 287 *
 288 *        Use DLACN2 to estimate the infinity-norm of the matrix
 289 *           inv(A) * diag(W),
 290 *        where W = abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) )))
 291 *
 292          DO 90 I = 1, N
 293             IF( WORK( I ).GT.SAFE2 ) THEN
 294                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
 295             ELSE
 296                WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
 297             END IF
 298    90    CONTINUE
 299 *
 300          KASE = 0
 301   100    CONTINUE
 302          CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
 303      $                KASE, ISAVE )
 304          IF( KASE.NE.0 ) THEN
 305             IF( KASE.EQ.1 ) THEN
 306 *
 307 *              Multiply by diag(W)*inv(A**T).
 308 *
 309                CALL DPBTRS( UPLO, N, KD, 1, AFB, LDAFB, WORK( N+1 ), N,
 310      $                      INFO )
 311                DO 110 I = 1, N
 312                   WORK( N+I ) = WORK( N+I )*WORK( I )
 313   110          CONTINUE
 314             ELSE IF( KASE.EQ.2 ) THEN
 315 *
 316 *              Multiply by inv(A)*diag(W).
 317 *
 318                DO 120 I = 1, N
 319                   WORK( N+I ) = WORK( N+I )*WORK( I )
 320   120          CONTINUE
 321                CALL DPBTRS( UPLO, N, KD, 1, AFB, LDAFB, WORK( N+1 ), N,
 322      $                      INFO )
 323             END IF
 324             GO TO 100
 325          END IF
 326 *
 327 *        Normalize error.
 328 *
 329          LSTRES = ZERO
 330          DO 130 I = 1, N
 331             LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
 332   130    CONTINUE
 333          IF( LSTRES.NE.ZERO )
 334      $      FERR( J ) = FERR( J ) / LSTRES
 335 *
 336   140 CONTINUE
 337 *
 338       RETURN
 339 *
 340 *     End of DPBRFS
 341 *
 342       END