Return to zsprfs.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

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