001:       SUBROUTINE DSPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
002:      +                   SWORK, ITER, INFO )
003: *
004: *  -- LAPACK PROTOTYPE driver routine (version 3.1.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
006: *     May 2007
007: *
008: *     ..
009: *     .. Scalar Arguments ..
010:       CHARACTER          UPLO
011:       INTEGER            INFO, ITER, LDA, LDB, LDX, N, NRHS
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               SWORK( * )
015:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( N, * ),
016:      +                   X( LDX, * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  DSPOSV computes the solution to a real system of linear equations
023: *     A * X = B,
024: *  where A is an N-by-N symmetric positive definite matrix and X and B
025: *  are N-by-NRHS matrices.
026: *
027: *  DSPOSV first attempts to factorize the matrix in SINGLE PRECISION
028: *  and use this factorization within an iterative refinement procedure
029: *  to produce a solution with DOUBLE PRECISION normwise backward error
030: *  quality (see below). If the approach fails the method switches to a
031: *  DOUBLE PRECISION factorization and solve.
032: *
033: *  The iterative refinement is not going to be a winning strategy if
034: *  the ratio SINGLE PRECISION performance over DOUBLE PRECISION
035: *  performance is too small. A reasonable strategy should take the
036: *  number of right-hand sides and the size of the matrix into account.
037: *  This might be done with a call to ILAENV in the future. Up to now, we
038: *  always try iterative refinement.
039: *
040: *  The iterative refinement process is stopped if
041: *      ITER > ITERMAX
042: *  or for all the RHS we have:
043: *      RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
044: *  where
045: *      o ITER is the number of the current iteration in the iterative
046: *        refinement process
047: *      o RNRM is the infinity-norm of the residual
048: *      o XNRM is the infinity-norm of the solution
049: *      o ANRM is the infinity-operator-norm of the matrix A
050: *      o EPS is the machine epsilon returned by DLAMCH('Epsilon')
051: *  The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
052: *  respectively.
053: *
054: *  Arguments
055: *  =========
056: *
057: *  UPLO    (input) CHARACTER
058: *          = 'U':  Upper triangle of A is stored;
059: *          = 'L':  Lower triangle of A is stored.
060: *
061: *  N       (input) INTEGER
062: *          The number of linear equations, i.e., the order of the
063: *          matrix A.  N >= 0.
064: *
065: *  NRHS    (input) INTEGER
066: *          The number of right hand sides, i.e., the number of columns
067: *          of the matrix B.  NRHS >= 0.
068: *
069: *  A       (input or input/ouptut) DOUBLE PRECISION array,
070: *          dimension (LDA,N)
071: *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
072: *          N-by-N upper triangular part of A contains the upper
073: *          triangular part of the matrix A, and the strictly lower
074: *          triangular part of A is not referenced.  If UPLO = 'L', the
075: *          leading N-by-N lower triangular part of A contains the lower
076: *          triangular part of the matrix A, and the strictly upper
077: *          triangular part of A is not referenced.
078: *          On exit, if iterative refinement has been successfully used
079: *          (INFO.EQ.0 and ITER.GE.0, see description below), then A is
080: *          unchanged, if double precision factorization has been used
081: *          (INFO.EQ.0 and ITER.LT.0, see description below), then the
082: *          array A contains the factor U or L from the Cholesky
083: *          factorization A = U**T*U or A = L*L**T.
084: *
085: *
086: *  LDA     (input) INTEGER
087: *          The leading dimension of the array A.  LDA >= max(1,N).
088: *
089: *  B       (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
090: *          The N-by-NRHS right hand side matrix B.
091: *
092: *  LDB     (input) INTEGER
093: *          The leading dimension of the array B.  LDB >= max(1,N).
094: *
095: *  X       (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
096: *          If INFO = 0, the N-by-NRHS solution matrix X.
097: *
098: *  LDX     (input) INTEGER
099: *          The leading dimension of the array X.  LDX >= max(1,N).
100: *
101: *  WORK    (workspace) DOUBLE PRECISION array, dimension (N*NRHS)
102: *          This array is used to hold the residual vectors.
103: *
104: *  SWORK   (workspace) REAL array, dimension (N*(N+NRHS))
105: *          This array is used to use the single precision matrix and the
106: *          right-hand sides or solutions in single precision.
107: *
108: *  ITER    (output) INTEGER
109: *          < 0: iterative refinement has failed, double precision
110: *               factorization has been performed
111: *               -1 : the routine fell back to full precision for
112: *                    implementation- or machine-specific reasons
113: *               -2 : narrowing the precision induced an overflow,
114: *                    the routine fell back to full precision
115: *               -3 : failure of SPOTRF
116: *               -31: stop the iterative refinement after the 30th
117: *                    iterations
118: *          > 0: iterative refinement has been sucessfully used.
119: *               Returns the number of iterations
120: *
121: *  INFO    (output) INTEGER
122: *          = 0:  successful exit
123: *          < 0:  if INFO = -i, the i-th argument had an illegal value
124: *          > 0:  if INFO = i, the leading minor of order i of (DOUBLE
125: *                PRECISION) A is not positive definite, so the
126: *                factorization could not be completed, and the solution
127: *                has not been computed.
128: *
129: *  =========
130: *
131: *     .. Parameters ..
132:       LOGICAL            DOITREF
133:       PARAMETER          ( DOITREF = .TRUE. )
134: *
135:       INTEGER            ITERMAX
136:       PARAMETER          ( ITERMAX = 30 )
137: *
138:       DOUBLE PRECISION   BWDMAX
139:       PARAMETER          ( BWDMAX = 1.0E+00 )
140: *
141:       DOUBLE PRECISION   NEGONE, ONE
142:       PARAMETER          ( NEGONE = -1.0D+0, ONE = 1.0D+0 )
143: *
144: *     .. Local Scalars ..
145:       INTEGER            I, IITER, PTSA, PTSX
146:       DOUBLE PRECISION   ANRM, CTE, EPS, RNRM, XNRM
147: *
148: *     .. External Subroutines ..
149:       EXTERNAL           DAXPY, DSYMM, DLACPY, DLAT2S, DLAG2S, SLAG2D,
150:      +                   SPOTRF, SPOTRS, XERBLA
151: *     ..
152: *     .. External Functions ..
153:       INTEGER            IDAMAX
154:       DOUBLE PRECISION   DLAMCH, DLANSY
155:       LOGICAL            LSAME
156:       EXTERNAL           IDAMAX, DLAMCH, DLANSY, LSAME
157: *     ..
158: *     .. Intrinsic Functions ..
159:       INTRINSIC          ABS, DBLE, MAX, SQRT
160: *     ..
161: *     .. Executable Statements ..
162: *
163:       INFO = 0
164:       ITER = 0
165: *
166: *     Test the input parameters.
167: *
168:       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
169:          INFO = -1
170:       ELSE IF( N.LT.0 ) THEN
171:          INFO = -2
172:       ELSE IF( NRHS.LT.0 ) THEN
173:          INFO = -3
174:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
175:          INFO = -5
176:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
177:          INFO = -7
178:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
179:          INFO = -9
180:       END IF
181:       IF( INFO.NE.0 ) THEN
182:          CALL XERBLA( 'DSPOSV', -INFO )
183:          RETURN
184:       END IF
185: *
186: *     Quick return if (N.EQ.0).
187: *
188:       IF( N.EQ.0 )
189:      +   RETURN
190: *
191: *     Skip single precision iterative refinement if a priori slower
192: *     than double precision factorization.
193: *
194:       IF( .NOT.DOITREF ) THEN
195:          ITER = -1
196:          GO TO 40
197:       END IF
198: *
199: *     Compute some constants.
200: *
201:       ANRM = DLANSY( 'I', UPLO, N, A, LDA, WORK )
202:       EPS = DLAMCH( 'Epsilon' )
203:       CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX
204: *
205: *     Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
206: *
207:       PTSA = 1
208:       PTSX = PTSA + N*N
209: *
210: *     Convert B from double precision to single precision and store the
211: *     result in SX.
212: *
213:       CALL DLAG2S( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO )
214: *
215:       IF( INFO.NE.0 ) THEN
216:          ITER = -2
217:          GO TO 40
218:       END IF
219: *
220: *     Convert A from double precision to single precision and store the
221: *     result in SA.
222: *
223:       CALL DLAT2S( UPLO, N, A, LDA, SWORK( PTSA ), N, INFO )
224: *
225:       IF( INFO.NE.0 ) THEN
226:          ITER = -2
227:          GO TO 40
228:       END IF
229: *
230: *     Compute the Cholesky factorization of SA.
231: *
232:       CALL SPOTRF( UPLO, N, SWORK( PTSA ), N, INFO )
233: *
234:       IF( INFO.NE.0 ) THEN
235:          ITER = -3
236:          GO TO 40
237:       END IF
238: *
239: *     Solve the system SA*SX = SB.
240: *
241:       CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
242:      +             INFO )
243: *
244: *     Convert SX back to double precision
245: *
246:       CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
247: *
248: *     Compute R = B - AX (R is WORK).
249: *
250:       CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
251: *
252:       CALL DSYMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
253:      +            WORK, N )
254: *
255: *     Check whether the NRHS normwise backward errors satisfy the
256: *     stopping criterion. If yes, set ITER=0 and return.
257: *
258:       DO I = 1, NRHS
259:          XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
260:          RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
261:          IF( RNRM.GT.XNRM*CTE )
262:      +      GO TO 10
263:       END DO
264: *
265: *     If we are here, the NRHS normwise backward errors satisfy the
266: *     stopping criterion. We are good to exit.
267: *
268:       ITER = 0
269:       RETURN
270: *
271:    10 CONTINUE
272: *
273:       DO 30 IITER = 1, ITERMAX
274: *
275: *        Convert R (in WORK) from double precision to single precision
276: *        and store the result in SX.
277: *
278:          CALL DLAG2S( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO )
279: *
280:          IF( INFO.NE.0 ) THEN
281:             ITER = -2
282:             GO TO 40
283:          END IF
284: *
285: *        Solve the system SA*SX = SR.
286: *
287:          CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
288:      +                INFO )
289: *
290: *        Convert SX back to double precision and update the current
291: *        iterate.
292: *
293:          CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
294: *
295:          DO I = 1, NRHS
296:             CALL DAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 )
297:          END DO
298: *
299: *        Compute R = B - AX (R is WORK).
300: *
301:          CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
302: *
303:          CALL DSYMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
304:      +               WORK, N )
305: *
306: *        Check whether the NRHS normwise backward errors satisfy the
307: *        stopping criterion. If yes, set ITER=IITER>0 and return.
308: *
309:          DO I = 1, NRHS
310:             XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
311:             RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
312:             IF( RNRM.GT.XNRM*CTE )
313:      +         GO TO 20
314:          END DO
315: *
316: *        If we are here, the NRHS normwise backward errors satisfy the
317: *        stopping criterion, we are good to exit.
318: *
319:          ITER = IITER
320: *
321:          RETURN
322: *
323:    20    CONTINUE
324: *
325:    30 CONTINUE
326: *
327: *     If we are at this place of the code, this is because we have
328: *     performed ITER=ITERMAX iterations and never satisified the
329: *     stopping criterion, set up the ITER flag accordingly and follow
330: *     up on double precision routine.
331: *
332:       ITER = -ITERMAX - 1
333: *
334:    40 CONTINUE
335: *
336: *     Single-precision iterative refinement failed to converge to a
337: *     satisfactory solution, so we resort to double precision.
338: *
339:       CALL DPOTRF( UPLO, N, A, LDA, INFO )
340: *
341:       IF( INFO.NE.0 )
342:      +   RETURN
343: *
344:       CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX )
345:       CALL DPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO )
346: *
347:       RETURN
348: *
349: *     End of DSPOSV.
350: *
351:       END
352: