LAPACK: SRC/cptrfs.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CPTRFS( UPLO, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
00002      $                   FERR, BERR, WORK, RWORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          UPLO
00011       INTEGER            INFO, LDB, LDX, N, NRHS
00012 *     ..
00013 *     .. Array Arguments ..
00014       REAL               BERR( * ), D( * ), DF( * ), FERR( * ),
00015      $                   RWORK( * )
00016       COMPLEX            B( LDB, * ), E( * ), EF( * ), WORK( * ),
00017      $                   X( LDX, * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  CPTRFS improves the computed solution to a system of linear
00024 *  equations when the coefficient matrix is Hermitian positive definite
00025 *  and tridiagonal, and provides error bounds and backward error
00026 *  estimates for the solution.
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  UPLO    (input) CHARACTER*1
00032 *          Specifies whether the superdiagonal or the subdiagonal of the
00033 *          tridiagonal matrix A is stored and the form of the
00034 *          factorization:
00035 *          = 'U':  E is the superdiagonal of A, and A = U**H*D*U;
00036 *          = 'L':  E is the subdiagonal of A, and A = L*D*L**H.
00037 *          (The two forms are equivalent if A is real.)
00038 *
00039 *  N       (input) INTEGER
00040 *          The order of the matrix A.  N >= 0.
00041 *
00042 *  NRHS    (input) INTEGER
00043 *          The number of right hand sides, i.e., the number of columns
00044 *          of the matrix B.  NRHS >= 0.
00045 *
00046 *  D       (input) REAL array, dimension (N)
00047 *          The n real diagonal elements of the tridiagonal matrix A.
00048 *
00049 *  E       (input) COMPLEX array, dimension (N-1)
00050 *          The (n-1) off-diagonal elements of the tridiagonal matrix A
00051 *          (see UPLO).
00052 *
00053 *  DF      (input) REAL array, dimension (N)
00054 *          The n diagonal elements of the diagonal matrix D from
00055 *          the factorization computed by CPTTRF.
00056 *
00057 *  EF      (input) COMPLEX array, dimension (N-1)
00058 *          The (n-1) off-diagonal elements of the unit bidiagonal
00059 *          factor U or L from the factorization computed by CPTTRF
00060 *          (see UPLO).
00061 *
00062 *  B       (input) COMPLEX array, dimension (LDB,NRHS)
00063 *          The right hand side matrix B.
00064 *
00065 *  LDB     (input) INTEGER
00066 *          The leading dimension of the array B.  LDB >= max(1,N).
00067 *
00068 *  X       (input/output) COMPLEX array, dimension (LDX,NRHS)
00069 *          On entry, the solution matrix X, as computed by CPTTRS.
00070 *          On exit, the improved solution matrix X.
00071 *
00072 *  LDX     (input) INTEGER
00073 *          The leading dimension of the array X.  LDX >= max(1,N).
00074 *
00075 *  FERR    (output) REAL array, dimension (NRHS)
00076 *          The forward error bound for each solution vector
00077 *          X(j) (the j-th column of the solution matrix X).
00078 *          If XTRUE is the true solution corresponding to X(j), FERR(j)
00079 *          is an estimated upper bound for the magnitude of the largest
00080 *          element in (X(j) - XTRUE) divided by the magnitude of the
00081 *          largest element in X(j).
00082 *
00083 *  BERR    (output) REAL array, dimension (NRHS)
00084 *          The componentwise relative backward error of each solution
00085 *          vector X(j) (i.e., the smallest relative change in
00086 *          any element of A or B that makes X(j) an exact solution).
00087 *
00088 *  WORK    (workspace) COMPLEX array, dimension (N)
00089 *
00090 *  RWORK   (workspace) REAL array, dimension (N)
00091 *
00092 *  INFO    (output) INTEGER
00093 *          = 0:  successful exit
00094 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00095 *
00096 *  Internal Parameters
00097 *  ===================
00098 *
00099 *  ITMAX is the maximum number of steps of iterative refinement.
00100 *
00101 *  =====================================================================
00102 *
00103 *     .. Parameters ..
00104       INTEGER            ITMAX
00105       PARAMETER          ( ITMAX = 5 )
00106       REAL               ZERO
00107       PARAMETER          ( ZERO = 0.0E+0 )
00108       REAL               ONE
00109       PARAMETER          ( ONE = 1.0E+0 )
00110       REAL               TWO
00111       PARAMETER          ( TWO = 2.0E+0 )
00112       REAL               THREE
00113       PARAMETER          ( THREE = 3.0E+0 )
00114 *     ..
00115 *     .. Local Scalars ..
00116       LOGICAL            UPPER
00117       INTEGER            COUNT, I, IX, J, NZ
00118       REAL               EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
00119       COMPLEX            BI, CX, DX, EX, ZDUM
00120 *     ..
00121 *     .. External Functions ..
00122       LOGICAL            LSAME
00123       INTEGER            ISAMAX
00124       REAL               SLAMCH
00125       EXTERNAL           LSAME, ISAMAX, SLAMCH
00126 *     ..
00127 *     .. External Subroutines ..
00128       EXTERNAL           CAXPY, CPTTRS, XERBLA
00129 *     ..
00130 *     .. Intrinsic Functions ..
00131       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, REAL
00132 *     ..
00133 *     .. Statement Functions ..
00134       REAL               CABS1
00135 *     ..
00136 *     .. Statement Function definitions ..
00137       CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
00138 *     ..
00139 *     .. Executable Statements ..
00140 *
00141 *     Test the input parameters.
00142 *
00143       INFO = 0
00144       UPPER = LSAME( UPLO, 'U' )
00145       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00146          INFO = -1
00147       ELSE IF( N.LT.0 ) THEN
00148          INFO = -2
00149       ELSE IF( NRHS.LT.0 ) THEN
00150          INFO = -3
00151       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00152          INFO = -9
00153       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00154          INFO = -11
00155       END IF
00156       IF( INFO.NE.0 ) THEN
00157          CALL XERBLA( 'CPTRFS', -INFO )
00158          RETURN
00159       END IF
00160 *
00161 *     Quick return if possible
00162 *
00163       IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
00164          DO 10 J = 1, NRHS
00165             FERR( J ) = ZERO
00166             BERR( J ) = ZERO
00167    10    CONTINUE
00168          RETURN
00169       END IF
00170 *
00171 *     NZ = maximum number of nonzero elements in each row of A, plus 1
00172 *
00173       NZ = 4
00174       EPS = SLAMCH( 'Epsilon' )
00175       SAFMIN = SLAMCH( 'Safe minimum' )
00176       SAFE1 = NZ*SAFMIN
00177       SAFE2 = SAFE1 / EPS
00178 *
00179 *     Do for each right hand side
00180 *
00181       DO 100 J = 1, NRHS
00182 *
00183          COUNT = 1
00184          LSTRES = THREE
00185    20    CONTINUE
00186 *
00187 *        Loop until stopping criterion is satisfied.
00188 *
00189 *        Compute residual R = B - A * X.  Also compute
00190 *        abs(A)*abs(x) + abs(b) for use in the backward error bound.
00191 *
00192          IF( UPPER ) THEN
00193             IF( N.EQ.1 ) THEN
00194                BI = B( 1, J )
00195                DX = D( 1 )*X( 1, J )
00196                WORK( 1 ) = BI - DX
00197                RWORK( 1 ) = CABS1( BI ) + CABS1( DX )
00198             ELSE
00199                BI = B( 1, J )
00200                DX = D( 1 )*X( 1, J )
00201                EX = E( 1 )*X( 2, J )
00202                WORK( 1 ) = BI - DX - EX
00203                RWORK( 1 ) = CABS1( BI ) + CABS1( DX ) +
00204      $                      CABS1( E( 1 ) )*CABS1( X( 2, J ) )
00205                DO 30 I = 2, N - 1
00206                   BI = B( I, J )
00207                   CX = CONJG( E( I-1 ) )*X( I-1, J )
00208                   DX = D( I )*X( I, J )
00209                   EX = E( I )*X( I+1, J )
00210                   WORK( I ) = BI - CX - DX - EX
00211                   RWORK( I ) = CABS1( BI ) +
00212      $                         CABS1( E( I-1 ) )*CABS1( X( I-1, J ) ) +
00213      $                         CABS1( DX ) + CABS1( E( I ) )*
00214      $                         CABS1( X( I+1, J ) )
00215    30          CONTINUE
00216                BI = B( N, J )
00217                CX = CONJG( E( N-1 ) )*X( N-1, J )
00218                DX = D( N )*X( N, J )
00219                WORK( N ) = BI - CX - DX
00220                RWORK( N ) = CABS1( BI ) + CABS1( E( N-1 ) )*
00221      $                      CABS1( X( N-1, J ) ) + CABS1( DX )
00222             END IF
00223          ELSE
00224             IF( N.EQ.1 ) THEN
00225                BI = B( 1, J )
00226                DX = D( 1 )*X( 1, J )
00227                WORK( 1 ) = BI - DX
00228                RWORK( 1 ) = CABS1( BI ) + CABS1( DX )
00229             ELSE
00230                BI = B( 1, J )
00231                DX = D( 1 )*X( 1, J )
00232                EX = CONJG( E( 1 ) )*X( 2, J )
00233                WORK( 1 ) = BI - DX - EX
00234                RWORK( 1 ) = CABS1( BI ) + CABS1( DX ) +
00235      $                      CABS1( E( 1 ) )*CABS1( X( 2, J ) )
00236                DO 40 I = 2, N - 1
00237                   BI = B( I, J )
00238                   CX = E( I-1 )*X( I-1, J )
00239                   DX = D( I )*X( I, J )
00240                   EX = CONJG( E( I ) )*X( I+1, J )
00241                   WORK( I ) = BI - CX - DX - EX
00242                   RWORK( I ) = CABS1( BI ) +
00243      $                         CABS1( E( I-1 ) )*CABS1( X( I-1, J ) ) +
00244      $                         CABS1( DX ) + CABS1( E( I ) )*
00245      $                         CABS1( X( I+1, J ) )
00246    40          CONTINUE
00247                BI = B( N, J )
00248                CX = E( N-1 )*X( N-1, J )
00249                DX = D( N )*X( N, J )
00250                WORK( N ) = BI - CX - DX
00251                RWORK( N ) = CABS1( BI ) + CABS1( E( N-1 ) )*
00252      $                      CABS1( X( N-1, J ) ) + CABS1( DX )
00253             END IF
00254          END IF
00255 *
00256 *        Compute componentwise relative backward error from formula
00257 *
00258 *        max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
00259 *
00260 *        where abs(Z) is the componentwise absolute value of the matrix
00261 *        or vector Z.  If the i-th component of the denominator is less
00262 *        than SAFE2, then SAFE1 is added to the i-th components of the
00263 *        numerator and denominator before dividing.
00264 *
00265          S = ZERO
00266          DO 50 I = 1, N
00267             IF( RWORK( I ).GT.SAFE2 ) THEN
00268                S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
00269             ELSE
00270                S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
00271      $             ( RWORK( I )+SAFE1 ) )
00272             END IF
00273    50    CONTINUE
00274          BERR( J ) = S
00275 *
00276 *        Test stopping criterion. Continue iterating if
00277 *           1) The residual BERR(J) is larger than machine epsilon, and
00278 *           2) BERR(J) decreased by at least a factor of 2 during the
00279 *              last iteration, and
00280 *           3) At most ITMAX iterations tried.
00281 *
00282          IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
00283      $       COUNT.LE.ITMAX ) THEN
00284 *
00285 *           Update solution and try again.
00286 *
00287             CALL CPTTRS( UPLO, N, 1, DF, EF, WORK, N, INFO )
00288             CALL CAXPY( N, CMPLX( ONE ), WORK, 1, X( 1, J ), 1 )
00289             LSTRES = BERR( J )
00290             COUNT = COUNT + 1
00291             GO TO 20
00292          END IF
00293 *
00294 *        Bound error from formula
00295 *
00296 *        norm(X - XTRUE) / norm(X) .le. FERR =
00297 *        norm( abs(inv(A))*
00298 *           ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
00299 *
00300 *        where
00301 *          norm(Z) is the magnitude of the largest component of Z
00302 *          inv(A) is the inverse of A
00303 *          abs(Z) is the componentwise absolute value of the matrix or
00304 *             vector Z
00305 *          NZ is the maximum number of nonzeros in any row of A, plus 1
00306 *          EPS is machine epsilon
00307 *
00308 *        The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
00309 *        is incremented by SAFE1 if the i-th component of
00310 *        abs(A)*abs(X) + abs(B) is less than SAFE2.
00311 *
00312          DO 60 I = 1, N
00313             IF( RWORK( I ).GT.SAFE2 ) THEN
00314                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
00315             ELSE
00316                RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
00317      $                      SAFE1
00318             END IF
00319    60    CONTINUE
00320          IX = ISAMAX( N, RWORK, 1 )
00321          FERR( J ) = RWORK( IX )
00322 *
00323 *        Estimate the norm of inv(A).
00324 *
00325 *        Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
00326 *
00327 *           m(i,j) =  abs(A(i,j)), i = j,
00328 *           m(i,j) = -abs(A(i,j)), i .ne. j,
00329 *
00330 *        and e = [ 1, 1, ..., 1 ]'.  Note M(A) = M(L)*D*M(L)'.
00331 *
00332 *        Solve M(L) * x = e.
00333 *
00334          RWORK( 1 ) = ONE
00335          DO 70 I = 2, N
00336             RWORK( I ) = ONE + RWORK( I-1 )*ABS( EF( I-1 ) )
00337    70    CONTINUE
00338 *
00339 *        Solve D * M(L)' * x = b.
00340 *
00341          RWORK( N ) = RWORK( N ) / DF( N )
00342          DO 80 I = N - 1, 1, -1
00343             RWORK( I ) = RWORK( I ) / DF( I ) +
00344      $                   RWORK( I+1 )*ABS( EF( I ) )
00345    80    CONTINUE
00346 *
00347 *        Compute norm(inv(A)) = max(x(i)), 1<=i<=n.
00348 *
00349          IX = ISAMAX( N, RWORK, 1 )
00350          FERR( J ) = FERR( J )*ABS( RWORK( IX ) )
00351 *
00352 *        Normalize error.
00353 *
00354          LSTRES = ZERO
00355          DO 90 I = 1, N
00356             LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
00357    90    CONTINUE
00358          IF( LSTRES.NE.ZERO )
00359      $      FERR( J ) = FERR( J ) / LSTRES
00360 *
00361   100 CONTINUE
00362 *
00363       RETURN
00364 *
00365 *     End of CPTRFS
00366 *
00367       END