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