 LAPACK 3.3.1 Linear Algebra PACKage

# sppsvx.f

Go to the documentation of this file.
```00001       SUBROUTINE SPPSVX( FACT, UPLO, N, NRHS, AP, AFP, EQUED, S, B, LDB,
00002      \$                   X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO )
00003 *
00004 *  -- LAPACK driver routine (version 3.3.1) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *  -- April 2011                                                      --
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          EQUED, FACT, UPLO
00011       INTEGER            INFO, LDB, LDX, N, NRHS
00012       REAL               RCOND
00013 *     ..
00014 *     .. Array Arguments ..
00015       INTEGER            IWORK( * )
00016       REAL               AFP( * ), AP( * ), B( LDB, * ), BERR( * ),
00017      \$                   FERR( * ), S( * ), WORK( * ), X( LDX, * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  SPPSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
00024 *  compute the solution to a real system of linear equations
00025 *     A * X = B,
00026 *  where A is an N-by-N symmetric positive definite matrix stored in
00027 *  packed format and X and B are N-by-NRHS matrices.
00028 *
00029 *  Error bounds on the solution and a condition estimate are also
00030 *  provided.
00031 *
00032 *  Description
00033 *  ===========
00034 *
00035 *  The following steps are performed:
00036 *
00037 *  1. If FACT = 'E', real scaling factors are computed to equilibrate
00038 *     the system:
00039 *        diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
00040 *     Whether or not the system will be equilibrated depends on the
00041 *     scaling of the matrix A, but if equilibration is used, A is
00042 *     overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
00043 *
00044 *  2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
00045 *     factor the matrix A (after equilibration if FACT = 'E') as
00046 *        A = U**T* U,  if UPLO = 'U', or
00047 *        A = L * L**T,  if UPLO = 'L',
00048 *     where U is an upper triangular matrix and L is a lower triangular
00049 *     matrix.
00050 *
00051 *  3. If the leading i-by-i principal minor is not positive definite,
00052 *     then the routine returns with INFO = i. Otherwise, the factored
00053 *     form of A is used to estimate the condition number of the matrix
00054 *     A.  If the reciprocal of the condition number is less than machine
00055 *     precision, INFO = N+1 is returned as a warning, but the routine
00056 *     still goes on to solve for X and compute error bounds as
00057 *     described below.
00058 *
00059 *  4. The system of equations is solved for X using the factored form
00060 *     of A.
00061 *
00062 *  5. Iterative refinement is applied to improve the computed solution
00063 *     matrix and calculate error bounds and backward error estimates
00064 *     for it.
00065 *
00066 *  6. If equilibration was used, the matrix X is premultiplied by
00067 *     diag(S) so that it solves the original system before
00068 *     equilibration.
00069 *
00070 *  Arguments
00071 *  =========
00072 *
00073 *  FACT    (input) CHARACTER*1
00074 *          Specifies whether or not the factored form of the matrix A is
00075 *          supplied on entry, and if not, whether the matrix A should be
00076 *          equilibrated before it is factored.
00077 *          = 'F':  On entry, AFP contains the factored form of A.
00078 *                  If EQUED = 'Y', the matrix A has been equilibrated
00079 *                  with scaling factors given by S.  AP and AFP will not
00080 *                  be modified.
00081 *          = 'N':  The matrix A will be copied to AFP and factored.
00082 *          = 'E':  The matrix A will be equilibrated if necessary, then
00083 *                  copied to AFP and factored.
00084 *
00085 *  UPLO    (input) CHARACTER*1
00086 *          = 'U':  Upper triangle of A is stored;
00087 *          = 'L':  Lower triangle of A is stored.
00088 *
00089 *  N       (input) INTEGER
00090 *          The number of linear equations, i.e., the order of the
00091 *          matrix A.  N >= 0.
00092 *
00093 *  NRHS    (input) INTEGER
00094 *          The number of right hand sides, i.e., the number of columns
00095 *          of the matrices B and X.  NRHS >= 0.
00096 *
00097 *  AP      (input/output) REAL array, dimension (N*(N+1)/2)
00098 *          On entry, the upper or lower triangle of the symmetric matrix
00099 *          A, packed columnwise in a linear array, except if FACT = 'F'
00100 *          and EQUED = 'Y', then A must contain the equilibrated matrix
00101 *          diag(S)*A*diag(S).  The j-th column of A is stored in the
00102 *          array AP as follows:
00103 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
00104 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
00105 *          See below for further details.  A is not modified if
00106 *          FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
00107 *
00108 *          On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
00109 *          diag(S)*A*diag(S).
00110 *
00111 *  AFP     (input or output) REAL array, dimension
00112 *                            (N*(N+1)/2)
00113 *          If FACT = 'F', then AFP is an input argument and on entry
00114 *          contains the triangular factor U or L from the Cholesky
00115 *          factorization A = U**T*U or A = L*L**T, in the same storage
00116 *          format as A.  If EQUED .ne. 'N', then AFP is the factored
00117 *          form of the equilibrated matrix A.
00118 *
00119 *          If FACT = 'N', then AFP is an output argument and on exit
00120 *          returns the triangular factor U or L from the Cholesky
00121 *          factorization A = U**T * U or A = L * L**T of the original
00122 *          matrix A.
00123 *
00124 *          If FACT = 'E', then AFP is an output argument and on exit
00125 *          returns the triangular factor U or L from the Cholesky
00126 *          factorization A = U**T * U or A = L * L**T of the equilibrated
00127 *          matrix A (see the description of AP for the form of the
00128 *          equilibrated matrix).
00129 *
00130 *  EQUED   (input or output) CHARACTER*1
00131 *          Specifies the form of equilibration that was done.
00132 *          = 'N':  No equilibration (always true if FACT = 'N').
00133 *          = 'Y':  Equilibration was done, i.e., A has been replaced by
00134 *                  diag(S) * A * diag(S).
00135 *          EQUED is an input argument if FACT = 'F'; otherwise, it is an
00136 *          output argument.
00137 *
00138 *  S       (input or output) REAL array, dimension (N)
00139 *          The scale factors for A; not accessed if EQUED = 'N'.  S is
00140 *          an input argument if FACT = 'F'; otherwise, S is an output
00141 *          argument.  If FACT = 'F' and EQUED = 'Y', each element of S
00142 *          must be positive.
00143 *
00144 *  B       (input/output) REAL array, dimension (LDB,NRHS)
00145 *          On entry, the N-by-NRHS right hand side matrix B.
00146 *          On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
00147 *          B is overwritten by diag(S) * B.
00148 *
00149 *  LDB     (input) INTEGER
00150 *          The leading dimension of the array B.  LDB >= max(1,N).
00151 *
00152 *  X       (output) REAL array, dimension (LDX,NRHS)
00153 *          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
00154 *          the original system of equations.  Note that if EQUED = 'Y',
00155 *          A and B are modified on exit, and the solution to the
00156 *          equilibrated system is inv(diag(S))*X.
00157 *
00158 *  LDX     (input) INTEGER
00159 *          The leading dimension of the array X.  LDX >= max(1,N).
00160 *
00161 *  RCOND   (output) REAL
00162 *          The estimate of the reciprocal condition number of the matrix
00163 *          A after equilibration (if done).  If RCOND is less than the
00164 *          machine precision (in particular, if RCOND = 0), the matrix
00165 *          is singular to working precision.  This condition is
00166 *          indicated by a return code of INFO > 0.
00167 *
00168 *  FERR    (output) REAL array, dimension (NRHS)
00169 *          The estimated forward error bound for each solution vector
00170 *          X(j) (the j-th column of the solution matrix X).
00171 *          If XTRUE is the true solution corresponding to X(j), FERR(j)
00172 *          is an estimated upper bound for the magnitude of the largest
00173 *          element in (X(j) - XTRUE) divided by the magnitude of the
00174 *          largest element in X(j).  The estimate is as reliable as
00175 *          the estimate for RCOND, and is almost always a slight
00176 *          overestimate of the true error.
00177 *
00178 *  BERR    (output) REAL array, dimension (NRHS)
00179 *          The componentwise relative backward error of each solution
00180 *          vector X(j) (i.e., the smallest relative change in
00181 *          any element of A or B that makes X(j) an exact solution).
00182 *
00183 *  WORK    (workspace) REAL array, dimension (3*N)
00184 *
00185 *  IWORK   (workspace) INTEGER array, dimension (N)
00186 *
00187 *  INFO    (output) INTEGER
00188 *          = 0:  successful exit
00189 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00190 *          > 0:  if INFO = i, and i is
00191 *                <= N:  the leading minor of order i of A is
00192 *                       not positive definite, so the factorization
00193 *                       could not be completed, and the solution has not
00194 *                       been computed. RCOND = 0 is returned.
00195 *                = N+1: U is nonsingular, but RCOND is less than machine
00196 *                       precision, meaning that the matrix is singular
00197 *                       to working precision.  Nevertheless, the
00198 *                       solution and error bounds are computed because
00199 *                       there are a number of situations where the
00200 *                       computed solution can be more accurate than the
00201 *                       value of RCOND would suggest.
00202 *
00203 *  Further Details
00204 *  ===============
00205 *
00206 *  The packed storage scheme is illustrated by the following example
00207 *  when N = 4, UPLO = 'U':
00208 *
00209 *  Two-dimensional storage of the symmetric matrix A:
00210 *
00211 *     a11 a12 a13 a14
00212 *         a22 a23 a24
00213 *             a33 a34     (aij = conjg(aji))
00214 *                 a44
00215 *
00216 *  Packed storage of the upper triangle of A:
00217 *
00218 *  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
00219 *
00220 *  =====================================================================
00221 *
00222 *     .. Parameters ..
00223       REAL               ZERO, ONE
00224       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00225 *     ..
00226 *     .. Local Scalars ..
00227       LOGICAL            EQUIL, NOFACT, RCEQU
00228       INTEGER            I, INFEQU, J
00229       REAL               AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
00230 *     ..
00231 *     .. External Functions ..
00232       LOGICAL            LSAME
00233       REAL               SLAMCH, SLANSP
00234       EXTERNAL           LSAME, SLAMCH, SLANSP
00235 *     ..
00236 *     .. External Subroutines ..
00237       EXTERNAL           SCOPY, SLACPY, SLAQSP, SPPCON, SPPEQU, SPPRFS,
00238      \$                   SPPTRF, SPPTRS, XERBLA
00239 *     ..
00240 *     .. Intrinsic Functions ..
00241       INTRINSIC          MAX, MIN
00242 *     ..
00243 *     .. Executable Statements ..
00244 *
00245       INFO = 0
00246       NOFACT = LSAME( FACT, 'N' )
00247       EQUIL = LSAME( FACT, 'E' )
00248       IF( NOFACT .OR. EQUIL ) THEN
00249          EQUED = 'N'
00250          RCEQU = .FALSE.
00251       ELSE
00252          RCEQU = LSAME( EQUED, 'Y' )
00253          SMLNUM = SLAMCH( 'Safe minimum' )
00254          BIGNUM = ONE / SMLNUM
00255       END IF
00256 *
00257 *     Test the input parameters.
00258 *
00259       IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
00260      \$     THEN
00261          INFO = -1
00262       ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
00263      \$          THEN
00264          INFO = -2
00265       ELSE IF( N.LT.0 ) THEN
00266          INFO = -3
00267       ELSE IF( NRHS.LT.0 ) THEN
00268          INFO = -4
00269       ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
00270      \$         ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
00271          INFO = -7
00272       ELSE
00273          IF( RCEQU ) THEN
00274             SMIN = BIGNUM
00275             SMAX = ZERO
00276             DO 10 J = 1, N
00277                SMIN = MIN( SMIN, S( J ) )
00278                SMAX = MAX( SMAX, S( J ) )
00279    10       CONTINUE
00280             IF( SMIN.LE.ZERO ) THEN
00281                INFO = -8
00282             ELSE IF( N.GT.0 ) THEN
00283                SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
00284             ELSE
00285                SCOND = ONE
00286             END IF
00287          END IF
00288          IF( INFO.EQ.0 ) THEN
00289             IF( LDB.LT.MAX( 1, N ) ) THEN
00290                INFO = -10
00291             ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00292                INFO = -12
00293             END IF
00294          END IF
00295       END IF
00296 *
00297       IF( INFO.NE.0 ) THEN
00298          CALL XERBLA( 'SPPSVX', -INFO )
00299          RETURN
00300       END IF
00301 *
00302       IF( EQUIL ) THEN
00303 *
00304 *        Compute row and column scalings to equilibrate the matrix A.
00305 *
00306          CALL SPPEQU( UPLO, N, AP, S, SCOND, AMAX, INFEQU )
00307          IF( INFEQU.EQ.0 ) THEN
00308 *
00309 *           Equilibrate the matrix.
00310 *
00311             CALL SLAQSP( UPLO, N, AP, S, SCOND, AMAX, EQUED )
00312             RCEQU = LSAME( EQUED, 'Y' )
00313          END IF
00314       END IF
00315 *
00316 *     Scale the right-hand side.
00317 *
00318       IF( RCEQU ) THEN
00319          DO 30 J = 1, NRHS
00320             DO 20 I = 1, N
00321                B( I, J ) = S( I )*B( I, J )
00322    20       CONTINUE
00323    30    CONTINUE
00324       END IF
00325 *
00326       IF( NOFACT .OR. EQUIL ) THEN
00327 *
00328 *        Compute the Cholesky factorization A = U**T * U or A = L * L**T.
00329 *
00330          CALL SCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 )
00331          CALL SPPTRF( UPLO, N, AFP, INFO )
00332 *
00333 *        Return if INFO is non-zero.
00334 *
00335          IF( INFO.GT.0 )THEN
00336             RCOND = ZERO
00337             RETURN
00338          END IF
00339       END IF
00340 *
00341 *     Compute the norm of the matrix A.
00342 *
00343       ANORM = SLANSP( 'I', UPLO, N, AP, WORK )
00344 *
00345 *     Compute the reciprocal of the condition number of A.
00346 *
00347       CALL SPPCON( UPLO, N, AFP, ANORM, RCOND, WORK, IWORK, INFO )
00348 *
00349 *     Compute the solution matrix X.
00350 *
00351       CALL SLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
00352       CALL SPPTRS( UPLO, N, NRHS, AFP, X, LDX, INFO )
00353 *
00354 *     Use iterative refinement to improve the computed solution and
00355 *     compute error bounds and backward error estimates for it.
00356 *
00357       CALL SPPRFS( UPLO, N, NRHS, AP, AFP, B, LDB, X, LDX, FERR, BERR,
00358      \$             WORK, IWORK, INFO )
00359 *
00360 *     Transform the solution matrix X to a solution of the original
00361 *     system.
00362 *
00363       IF( RCEQU ) THEN
00364          DO 50 J = 1, NRHS
00365             DO 40 I = 1, N
00366                X( I, J ) = S( I )*X( I, J )
00367    40       CONTINUE
00368    50    CONTINUE
00369          DO 60 J = 1, NRHS
00370             FERR( J ) = FERR( J ) / SCOND
00371    60    CONTINUE
00372       END IF
00373 *
00374 *     Set INFO = N+1 if the matrix is singular to working precision.
00375 *
00376       IF( RCOND.LT.SLAMCH( 'Epsilon' ) )
00377      \$   INFO = N + 1
00378 *
00379       RETURN
00380 *
00381 *     End of SPPSVX
00382 *
00383       END
```