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