LAPACK 3.3.1 Linear Algebra PACKage

# cppsvx.f

Go to the documentation of this file.
```00001       SUBROUTINE CPPSVX( FACT, UPLO, N, NRHS, AP, AFP, EQUED, S, B, LDB,
00002      \$                   X, LDX, RCOND, FERR, BERR, WORK, RWORK, 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       REAL               BERR( * ), FERR( * ), RWORK( * ), S( * )
00016       COMPLEX            AFP( * ), AP( * ), B( LDB, * ), WORK( * ),
00017      \$                   X( LDX, * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  CPPSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
00024 *  compute the solution to a complex system of linear equations
00025 *     A * X = B,
00026 *  where A is an N-by-N Hermitian 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**H * U ,  if UPLO = 'U', or
00047 *        A = L * L**H,  if UPLO = 'L',
00048 *     where U is an upper triangular matrix, L is a lower triangular
00049 *     matrix, and **H indicates conjugate transpose.
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) COMPLEX array, dimension (N*(N+1)/2)
00098 *          On entry, the upper or lower triangle of the Hermitian 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) COMPLEX array, dimension (N*(N+1)/2)
00112 *          If FACT = 'F', then AFP is an input argument and on entry
00113 *          contains the triangular factor U or L from the Cholesky
00114 *          factorization A = U**H*U or A = L*L**H, in the same storage
00115 *          format as A.  If EQUED .ne. 'N', then AFP is the factored
00116 *          form of the equilibrated matrix A.
00117 *
00118 *          If FACT = 'N', then AFP is an output argument and on exit
00119 *          returns the triangular factor U or L from the Cholesky
00120 *          factorization A = U**H * U or A = L * L**H of the original
00121 *          matrix A.
00122 *
00123 *          If FACT = 'E', then AFP is an output argument and on exit
00124 *          returns the triangular factor U or L from the Cholesky
00125 *          factorization A = U**H*U or A = L*L**H of the equilibrated
00126 *          matrix A (see the description of AP for the form of the
00127 *          equilibrated matrix).
00128 *
00129 *  EQUED   (input or output) CHARACTER*1
00130 *          Specifies the form of equilibration that was done.
00131 *          = 'N':  No equilibration (always true if FACT = 'N').
00132 *          = 'Y':  Equilibration was done, i.e., A has been replaced by
00133 *                  diag(S) * A * diag(S).
00134 *          EQUED is an input argument if FACT = 'F'; otherwise, it is an
00135 *          output argument.
00136 *
00137 *  S       (input or output) REAL array, dimension (N)
00138 *          The scale factors for A; not accessed if EQUED = 'N'.  S is
00139 *          an input argument if FACT = 'F'; otherwise, S is an output
00140 *          argument.  If FACT = 'F' and EQUED = 'Y', each element of S
00141 *          must be positive.
00142 *
00143 *  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
00144 *          On entry, the N-by-NRHS right hand side matrix B.
00145 *          On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
00146 *          B is overwritten by diag(S) * B.
00147 *
00148 *  LDB     (input) INTEGER
00149 *          The leading dimension of the array B.  LDB >= max(1,N).
00150 *
00151 *  X       (output) COMPLEX array, dimension (LDX,NRHS)
00152 *          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
00153 *          the original system of equations.  Note that if EQUED = 'Y',
00154 *          A and B are modified on exit, and the solution to the
00155 *          equilibrated system is inv(diag(S))*X.
00156 *
00157 *  LDX     (input) INTEGER
00158 *          The leading dimension of the array X.  LDX >= max(1,N).
00159 *
00160 *  RCOND   (output) REAL
00161 *          The estimate of the reciprocal condition number of the matrix
00162 *          A after equilibration (if done).  If RCOND is less than the
00163 *          machine precision (in particular, if RCOND = 0), the matrix
00164 *          is singular to working precision.  This condition is
00165 *          indicated by a return code of INFO > 0.
00166 *
00167 *  FERR    (output) REAL array, dimension (NRHS)
00168 *          The estimated forward error bound for each solution vector
00169 *          X(j) (the j-th column of the solution matrix X).
00170 *          If XTRUE is the true solution corresponding to X(j), FERR(j)
00171 *          is an estimated upper bound for the magnitude of the largest
00172 *          element in (X(j) - XTRUE) divided by the magnitude of the
00173 *          largest element in X(j).  The estimate is as reliable as
00174 *          the estimate for RCOND, and is almost always a slight
00175 *          overestimate of the true error.
00176 *
00177 *  BERR    (output) REAL array, dimension (NRHS)
00178 *          The componentwise relative backward error of each solution
00179 *          vector X(j) (i.e., the smallest relative change in
00180 *          any element of A or B that makes X(j) an exact solution).
00181 *
00182 *  WORK    (workspace) COMPLEX array, dimension (2*N)
00183 *
00184 *  RWORK   (workspace) REAL array, dimension (N)
00185 *
00186 *  INFO    (output) INTEGER
00187 *          = 0:  successful exit
00188 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00189 *          > 0:  if INFO = i, and i is
00190 *                <= N:  the leading minor of order i of A is
00191 *                       not positive definite, so the factorization
00192 *                       could not be completed, and the solution has not
00193 *                       been computed. RCOND = 0 is returned.
00194 *                = N+1: U is nonsingular, but RCOND is less than machine
00195 *                       precision, meaning that the matrix is singular
00196 *                       to working precision.  Nevertheless, the
00197 *                       solution and error bounds are computed because
00198 *                       there are a number of situations where the
00199 *                       computed solution can be more accurate than the
00200 *                       value of RCOND would suggest.
00201 *
00202 *  Further Details
00203 *  ===============
00204 *
00205 *  The packed storage scheme is illustrated by the following example
00206 *  when N = 4, UPLO = 'U':
00207 *
00208 *  Two-dimensional storage of the Hermitian matrix A:
00209 *
00210 *     a11 a12 a13 a14
00211 *         a22 a23 a24
00212 *             a33 a34     (aij = conjg(aji))
00213 *                 a44
00214 *
00215 *  Packed storage of the upper triangle of A:
00216 *
00217 *  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
00218 *
00219 *  =====================================================================
00220 *
00221 *     .. Parameters ..
00222       REAL               ZERO, ONE
00223       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00224 *     ..
00225 *     .. Local Scalars ..
00226       LOGICAL            EQUIL, NOFACT, RCEQU
00227       INTEGER            I, INFEQU, J
00228       REAL               AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
00229 *     ..
00230 *     .. External Functions ..
00231       LOGICAL            LSAME
00232       REAL               CLANHP, SLAMCH
00233       EXTERNAL           LSAME, CLANHP, SLAMCH
00234 *     ..
00235 *     .. External Subroutines ..
00236       EXTERNAL           CCOPY, CLACPY, CLAQHP, CPPCON, CPPEQU, CPPRFS,
00237      \$                   CPPTRF, CPPTRS, XERBLA
00238 *     ..
00239 *     .. Intrinsic Functions ..
00240       INTRINSIC          MAX, MIN
00241 *     ..
00242 *     .. Executable Statements ..
00243 *
00244       INFO = 0
00245       NOFACT = LSAME( FACT, 'N' )
00246       EQUIL = LSAME( FACT, 'E' )
00247       IF( NOFACT .OR. EQUIL ) THEN
00248          EQUED = 'N'
00249          RCEQU = .FALSE.
00250       ELSE
00251          RCEQU = LSAME( EQUED, 'Y' )
00252          SMLNUM = SLAMCH( 'Safe minimum' )
00253          BIGNUM = ONE / SMLNUM
00254       END IF
00255 *
00256 *     Test the input parameters.
00257 *
00258       IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
00259      \$     THEN
00260          INFO = -1
00261       ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
00262      \$          THEN
00263          INFO = -2
00264       ELSE IF( N.LT.0 ) THEN
00265          INFO = -3
00266       ELSE IF( NRHS.LT.0 ) THEN
00267          INFO = -4
00268       ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
00269      \$         ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
00270          INFO = -7
00271       ELSE
00272          IF( RCEQU ) THEN
00273             SMIN = BIGNUM
00274             SMAX = ZERO
00275             DO 10 J = 1, N
00276                SMIN = MIN( SMIN, S( J ) )
00277                SMAX = MAX( SMAX, S( J ) )
00278    10       CONTINUE
00279             IF( SMIN.LE.ZERO ) THEN
00280                INFO = -8
00281             ELSE IF( N.GT.0 ) THEN
00282                SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
00283             ELSE
00284                SCOND = ONE
00285             END IF
00286          END IF
00287          IF( INFO.EQ.0 ) THEN
00288             IF( LDB.LT.MAX( 1, N ) ) THEN
00289                INFO = -10
00290             ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00291                INFO = -12
00292             END IF
00293          END IF
00294       END IF
00295 *
00296       IF( INFO.NE.0 ) THEN
00297          CALL XERBLA( 'CPPSVX', -INFO )
00298          RETURN
00299       END IF
00300 *
00301       IF( EQUIL ) THEN
00302 *
00303 *        Compute row and column scalings to equilibrate the matrix A.
00304 *
00305          CALL CPPEQU( UPLO, N, AP, S, SCOND, AMAX, INFEQU )
00306          IF( INFEQU.EQ.0 ) THEN
00307 *
00308 *           Equilibrate the matrix.
00309 *
00310             CALL CLAQHP( UPLO, N, AP, S, SCOND, AMAX, EQUED )
00311             RCEQU = LSAME( EQUED, 'Y' )
00312          END IF
00313       END IF
00314 *
00315 *     Scale the right-hand side.
00316 *
00317       IF( RCEQU ) THEN
00318          DO 30 J = 1, NRHS
00319             DO 20 I = 1, N
00320                B( I, J ) = S( I )*B( I, J )
00321    20       CONTINUE
00322    30    CONTINUE
00323       END IF
00324 *
00325       IF( NOFACT .OR. EQUIL ) THEN
00326 *
00327 *        Compute the Cholesky factorization A = U**H * U or A = L * L**H.
00328 *
00329          CALL CCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 )
00330          CALL CPPTRF( UPLO, N, AFP, INFO )
00331 *
00332 *        Return if INFO is non-zero.
00333 *
00334          IF( INFO.GT.0 )THEN
00335             RCOND = ZERO
00336             RETURN
00337          END IF
00338       END IF
00339 *
00340 *     Compute the norm of the matrix A.
00341 *
00342       ANORM = CLANHP( 'I', UPLO, N, AP, RWORK )
00343 *
00344 *     Compute the reciprocal of the condition number of A.
00345 *
00346       CALL CPPCON( UPLO, N, AFP, ANORM, RCOND, WORK, RWORK, INFO )
00347 *
00348 *     Compute the solution matrix X.
00349 *
00350       CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
00351       CALL CPPTRS( UPLO, N, NRHS, AFP, X, LDX, INFO )
00352 *
00353 *     Use iterative refinement to improve the computed solution and
00354 *     compute error bounds and backward error estimates for it.
00355 *
00356       CALL CPPRFS( UPLO, N, NRHS, AP, AFP, B, LDB, X, LDX, FERR, BERR,
00357      \$             WORK, RWORK, INFO )
00358 *
00359 *     Transform the solution matrix X to a solution of the original
00360 *     system.
00361 *
00362       IF( RCEQU ) THEN
00363          DO 50 J = 1, NRHS
00364             DO 40 I = 1, N
00365                X( I, J ) = S( I )*X( I, J )
00366    40       CONTINUE
00367    50    CONTINUE
00368          DO 60 J = 1, NRHS
00369             FERR( J ) = FERR( J ) / SCOND
00370    60    CONTINUE
00371       END IF
00372 *
00373 *     Set INFO = N+1 if the matrix is singular to working precision.
00374 *
00375       IF( RCOND.LT.SLAMCH( 'Epsilon' ) )
00376      \$   INFO = N + 1
00377 *
00378       RETURN
00379 *
00380 *     End of CPPSVX
00381 *
00382       END
```