001:       SUBROUTINE SPPSVX( FACT, UPLO, N, NRHS, AP, AFP, EQUED, S, B, LDB,
002:      $                   X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          EQUED, FACT, UPLO
010:       INTEGER            INFO, LDB, LDX, N, NRHS
011:       REAL               RCOND
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IWORK( * )
015:       REAL               AFP( * ), AP( * ), B( LDB, * ), BERR( * ),
016:      $                   FERR( * ), S( * ), WORK( * ), X( LDX, * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  SPPSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
023: *  compute the solution to a real system of linear equations
024: *     A * X = B,
025: *  where A is an N-by-N symmetric positive definite matrix stored in
026: *  packed format and X and B are N-by-NRHS matrices.
027: *
028: *  Error bounds on the solution and a condition estimate are also
029: *  provided.
030: *
031: *  Description
032: *  ===========
033: *
034: *  The following steps are performed:
035: *
036: *  1. If FACT = 'E', real scaling factors are computed to equilibrate
037: *     the system:
038: *        diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
039: *     Whether or not the system will be equilibrated depends on the
040: *     scaling of the matrix A, but if equilibration is used, A is
041: *     overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
042: *
043: *  2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
044: *     factor the matrix A (after equilibration if FACT = 'E') as
045: *        A = U**T* U,  if UPLO = 'U', or
046: *        A = L * L**T,  if UPLO = 'L',
047: *     where U is an upper triangular matrix and L is a lower triangular
048: *     matrix.
049: *
050: *  3. If the leading i-by-i principal minor is not positive definite,
051: *     then the routine returns with INFO = i. Otherwise, the factored
052: *     form of A is used to estimate the condition number of the matrix
053: *     A.  If the reciprocal of the condition number is less than machine
054: *     precision, INFO = N+1 is returned as a warning, but the routine
055: *     still goes on to solve for X and compute error bounds as
056: *     described below.
057: *
058: *  4. The system of equations is solved for X using the factored form
059: *     of A.
060: *
061: *  5. Iterative refinement is applied to improve the computed solution
062: *     matrix and calculate error bounds and backward error estimates
063: *     for it.
064: *
065: *  6. If equilibration was used, the matrix X is premultiplied by
066: *     diag(S) so that it solves the original system before
067: *     equilibration.
068: *
069: *  Arguments
070: *  =========
071: *
072: *  FACT    (input) CHARACTER*1
073: *          Specifies whether or not the factored form of the matrix A is
074: *          supplied on entry, and if not, whether the matrix A should be
075: *          equilibrated before it is factored.
076: *          = 'F':  On entry, AFP contains the factored form of A.
077: *                  If EQUED = 'Y', the matrix A has been equilibrated
078: *                  with scaling factors given by S.  AP and AFP will not
079: *                  be modified.
080: *          = 'N':  The matrix A will be copied to AFP and factored.
081: *          = 'E':  The matrix A will be equilibrated if necessary, then
082: *                  copied to AFP and factored.
083: *
084: *  UPLO    (input) CHARACTER*1
085: *          = 'U':  Upper triangle of A is stored;
086: *          = 'L':  Lower triangle of A is stored.
087: *
088: *  N       (input) INTEGER
089: *          The number of linear equations, i.e., the order of the
090: *          matrix A.  N >= 0.
091: *
092: *  NRHS    (input) INTEGER
093: *          The number of right hand sides, i.e., the number of columns
094: *          of the matrices B and X.  NRHS >= 0.
095: *
096: *  AP      (input/output) REAL array, dimension (N*(N+1)/2)
097: *          On entry, the upper or lower triangle of the symmetric matrix
098: *          A, packed columnwise in a linear array, except if FACT = 'F'
099: *          and EQUED = 'Y', then A must contain the equilibrated matrix
100: *          diag(S)*A*diag(S).  The j-th column of A is stored in the
101: *          array AP as follows:
102: *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
103: *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
104: *          See below for further details.  A is not modified if
105: *          FACT = 'F' or 'N', or if FACT = 'E' and EQUED = 'N' on exit.
106: *
107: *          On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
108: *          diag(S)*A*diag(S).
109: *
110: *  AFP     (input or output) REAL array, dimension
111: *                            (N*(N+1)/2)
112: *          If FACT = 'F', then AFP is an input argument and on entry
113: *          contains the triangular factor U or L from the Cholesky
114: *          factorization A = U'*U or A = L*L', in the same storage
115: *          format as A.  If EQUED .ne. 'N', then AFP is the factored
116: *          form of the equilibrated matrix A.
117: *
118: *          If FACT = 'N', then AFP is an output argument and on exit
119: *          returns the triangular factor U or L from the Cholesky
120: *          factorization A = U'*U or A = L*L' of the original matrix A.
121: *
122: *          If FACT = 'E', then AFP is an output argument and on exit
123: *          returns the triangular factor U or L from the Cholesky
124: *          factorization A = U'*U or A = L*L' of the equilibrated
125: *          matrix A (see the description of AP for the form of the
126: *          equilibrated matrix).
127: *
128: *  EQUED   (input or output) CHARACTER*1
129: *          Specifies the form of equilibration that was done.
130: *          = 'N':  No equilibration (always true if FACT = 'N').
131: *          = 'Y':  Equilibration was done, i.e., A has been replaced by
132: *                  diag(S) * A * diag(S).
133: *          EQUED is an input argument if FACT = 'F'; otherwise, it is an
134: *          output argument.
135: *
136: *  S       (input or output) REAL array, dimension (N)
137: *          The scale factors for A; not accessed if EQUED = 'N'.  S is
138: *          an input argument if FACT = 'F'; otherwise, S is an output
139: *          argument.  If FACT = 'F' and EQUED = 'Y', each element of S
140: *          must be positive.
141: *
142: *  B       (input/output) REAL array, dimension (LDB,NRHS)
143: *          On entry, the N-by-NRHS right hand side matrix B.
144: *          On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
145: *          B is overwritten by diag(S) * B.
146: *
147: *  LDB     (input) INTEGER
148: *          The leading dimension of the array B.  LDB >= max(1,N).
149: *
150: *  X       (output) REAL array, dimension (LDX,NRHS)
151: *          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
152: *          the original system of equations.  Note that if EQUED = 'Y',
153: *          A and B are modified on exit, and the solution to the
154: *          equilibrated system is inv(diag(S))*X.
155: *
156: *  LDX     (input) INTEGER
157: *          The leading dimension of the array X.  LDX >= max(1,N).
158: *
159: *  RCOND   (output) REAL
160: *          The estimate of the reciprocal condition number of the matrix
161: *          A after equilibration (if done).  If RCOND is less than the
162: *          machine precision (in particular, if RCOND = 0), the matrix
163: *          is singular to working precision.  This condition is
164: *          indicated by a return code of INFO > 0.
165: *
166: *  FERR    (output) REAL array, dimension (NRHS)
167: *          The estimated forward error bound for each solution vector
168: *          X(j) (the j-th column of the solution matrix X).
169: *          If XTRUE is the true solution corresponding to X(j), FERR(j)
170: *          is an estimated upper bound for the magnitude of the largest
171: *          element in (X(j) - XTRUE) divided by the magnitude of the
172: *          largest element in X(j).  The estimate is as reliable as
173: *          the estimate for RCOND, and is almost always a slight
174: *          overestimate of the true error.
175: *
176: *  BERR    (output) REAL array, dimension (NRHS)
177: *          The componentwise relative backward error of each solution
178: *          vector X(j) (i.e., the smallest relative change in
179: *          any element of A or B that makes X(j) an exact solution).
180: *
181: *  WORK    (workspace) REAL array, dimension (3*N)
182: *
183: *  IWORK   (workspace) INTEGER array, dimension (N)
184: *
185: *  INFO    (output) INTEGER
186: *          = 0:  successful exit
187: *          < 0:  if INFO = -i, the i-th argument had an illegal value
188: *          > 0:  if INFO = i, and i is
189: *                <= N:  the leading minor of order i of A is
190: *                       not positive definite, so the factorization
191: *                       could not be completed, and the solution has not
192: *                       been computed. RCOND = 0 is returned.
193: *                = N+1: U is nonsingular, but RCOND is less than machine
194: *                       precision, meaning that the matrix is singular
195: *                       to working precision.  Nevertheless, the
196: *                       solution and error bounds are computed because
197: *                       there are a number of situations where the
198: *                       computed solution can be more accurate than the
199: *                       value of RCOND would suggest.
200: *
201: *  Further Details
202: *  ===============
203: *
204: *  The packed storage scheme is illustrated by the following example
205: *  when N = 4, UPLO = 'U':
206: *
207: *  Two-dimensional storage of the symmetric matrix A:
208: *
209: *     a11 a12 a13 a14
210: *         a22 a23 a24
211: *             a33 a34     (aij = conjg(aji))
212: *                 a44
213: *
214: *  Packed storage of the upper triangle of A:
215: *
216: *  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
217: *
218: *  =====================================================================
219: *
220: *     .. Parameters ..
221:       REAL               ZERO, ONE
222:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
223: *     ..
224: *     .. Local Scalars ..
225:       LOGICAL            EQUIL, NOFACT, RCEQU
226:       INTEGER            I, INFEQU, J
227:       REAL               AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
228: *     ..
229: *     .. External Functions ..
230:       LOGICAL            LSAME
231:       REAL               SLAMCH, SLANSP
232:       EXTERNAL           LSAME, SLAMCH, SLANSP
233: *     ..
234: *     .. External Subroutines ..
235:       EXTERNAL           SCOPY, SLACPY, SLAQSP, SPPCON, SPPEQU, SPPRFS,
236:      $                   SPPTRF, SPPTRS, XERBLA
237: *     ..
238: *     .. Intrinsic Functions ..
239:       INTRINSIC          MAX, MIN
240: *     ..
241: *     .. Executable Statements ..
242: *
243:       INFO = 0
244:       NOFACT = LSAME( FACT, 'N' )
245:       EQUIL = LSAME( FACT, 'E' )
246:       IF( NOFACT .OR. EQUIL ) THEN
247:          EQUED = 'N'
248:          RCEQU = .FALSE.
249:       ELSE
250:          RCEQU = LSAME( EQUED, 'Y' )
251:          SMLNUM = SLAMCH( 'Safe minimum' )
252:          BIGNUM = ONE / SMLNUM
253:       END IF
254: *
255: *     Test the input parameters.
256: *
257:       IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
258:      $     THEN
259:          INFO = -1
260:       ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
261:      $          THEN
262:          INFO = -2
263:       ELSE IF( N.LT.0 ) THEN
264:          INFO = -3
265:       ELSE IF( NRHS.LT.0 ) THEN
266:          INFO = -4
267:       ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
268:      $         ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
269:          INFO = -7
270:       ELSE
271:          IF( RCEQU ) THEN
272:             SMIN = BIGNUM
273:             SMAX = ZERO
274:             DO 10 J = 1, N
275:                SMIN = MIN( SMIN, S( J ) )
276:                SMAX = MAX( SMAX, S( J ) )
277:    10       CONTINUE
278:             IF( SMIN.LE.ZERO ) THEN
279:                INFO = -8
280:             ELSE IF( N.GT.0 ) THEN
281:                SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
282:             ELSE
283:                SCOND = ONE
284:             END IF
285:          END IF
286:          IF( INFO.EQ.0 ) THEN
287:             IF( LDB.LT.MAX( 1, N ) ) THEN
288:                INFO = -10
289:             ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
290:                INFO = -12
291:             END IF
292:          END IF
293:       END IF
294: *
295:       IF( INFO.NE.0 ) THEN
296:          CALL XERBLA( 'SPPSVX', -INFO )
297:          RETURN
298:       END IF
299: *
300:       IF( EQUIL ) THEN
301: *
302: *        Compute row and column scalings to equilibrate the matrix A.
303: *
304:          CALL SPPEQU( UPLO, N, AP, S, SCOND, AMAX, INFEQU )
305:          IF( INFEQU.EQ.0 ) THEN
306: *
307: *           Equilibrate the matrix.
308: *
309:             CALL SLAQSP( UPLO, N, AP, S, SCOND, AMAX, EQUED )
310:             RCEQU = LSAME( EQUED, 'Y' )
311:          END IF
312:       END IF
313: *
314: *     Scale the right-hand side.
315: *
316:       IF( RCEQU ) THEN
317:          DO 30 J = 1, NRHS
318:             DO 20 I = 1, N
319:                B( I, J ) = S( I )*B( I, J )
320:    20       CONTINUE
321:    30    CONTINUE
322:       END IF
323: *
324:       IF( NOFACT .OR. EQUIL ) THEN
325: *
326: *        Compute the Cholesky factorization A = U'*U or A = L*L'.
327: *
328:          CALL SCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 )
329:          CALL SPPTRF( UPLO, N, AFP, INFO )
330: *
331: *        Return if INFO is non-zero.
332: *
333:          IF( INFO.GT.0 )THEN
334:             RCOND = ZERO
335:             RETURN
336:          END IF
337:       END IF
338: *
339: *     Compute the norm of the matrix A.
340: *
341:       ANORM = SLANSP( 'I', UPLO, N, AP, WORK )
342: *
343: *     Compute the reciprocal of the condition number of A.
344: *
345:       CALL SPPCON( UPLO, N, AFP, ANORM, RCOND, WORK, IWORK, INFO )
346: *
347: *     Compute the solution matrix X.
348: *
349:       CALL SLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
350:       CALL SPPTRS( UPLO, N, NRHS, AFP, X, LDX, INFO )
351: *
352: *     Use iterative refinement to improve the computed solution and
353: *     compute error bounds and backward error estimates for it.
354: *
355:       CALL SPPRFS( UPLO, N, NRHS, AP, AFP, B, LDB, X, LDX, FERR, BERR,
356:      $             WORK, IWORK, INFO )
357: *
358: *     Transform the solution matrix X to a solution of the original
359: *     system.
360: *
361:       IF( RCEQU ) THEN
362:          DO 50 J = 1, NRHS
363:             DO 40 I = 1, N
364:                X( I, J ) = S( I )*X( I, J )
365:    40       CONTINUE
366:    50    CONTINUE
367:          DO 60 J = 1, NRHS
368:             FERR( J ) = FERR( J ) / SCOND
369:    60    CONTINUE
370:       END IF
371: *
372: *     Set INFO = N+1 if the matrix is singular to working precision.
373: *
374:       IF( RCOND.LT.SLAMCH( 'Epsilon' ) )
375:      $   INFO = N + 1
376: *
377:       RETURN
378: *
379: *     End of SPPSVX
380: *
381:       END
382: