1: *> \brief <b> DSPOSV computes the solution to system of linear equations A * X = B for PO matrices</b>
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DSPOSV + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsposv.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsposv.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsposv.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DSPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
22: * SWORK, ITER, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER UPLO
26: * INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
27: * ..
28: * .. Array Arguments ..
29: * REAL SWORK( * )
30: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ),
31: * $ X( LDX, * )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> DSPOSV computes the solution to a real system of linear equations
41: *> A * X = B,
42: *> where A is an N-by-N symmetric positive definite matrix and X and B
43: *> are N-by-NRHS matrices.
44: *>
45: *> DSPOSV first attempts to factorize the matrix in SINGLE PRECISION
46: *> and use this factorization within an iterative refinement procedure
47: *> to produce a solution with DOUBLE PRECISION normwise backward error
48: *> quality (see below). If the approach fails the method switches to a
49: *> DOUBLE PRECISION factorization and solve.
50: *>
51: *> The iterative refinement is not going to be a winning strategy if
52: *> the ratio SINGLE PRECISION performance over DOUBLE PRECISION
53: *> performance is too small. A reasonable strategy should take the
54: *> number of right-hand sides and the size of the matrix into account.
55: *> This might be done with a call to ILAENV in the future. Up to now, we
56: *> always try iterative refinement.
57: *>
58: *> The iterative refinement process is stopped if
59: *> ITER > ITERMAX
60: *> or for all the RHS we have:
61: *> RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
62: *> where
63: *> o ITER is the number of the current iteration in the iterative
64: *> refinement process
65: *> o RNRM is the infinity-norm of the residual
66: *> o XNRM is the infinity-norm of the solution
67: *> o ANRM is the infinity-operator-norm of the matrix A
68: *> o EPS is the machine epsilon returned by DLAMCH('Epsilon')
69: *> The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
70: *> respectively.
71: *> \endverbatim
72: *
73: * Arguments:
74: * ==========
75: *
76: *> \param[in] UPLO
77: *> \verbatim
78: *> UPLO is CHARACTER*1
79: *> = 'U': Upper triangle of A is stored;
80: *> = 'L': Lower triangle of A is stored.
81: *> \endverbatim
82: *>
83: *> \param[in] N
84: *> \verbatim
85: *> N is INTEGER
86: *> The number of linear equations, i.e., the order of the
87: *> matrix A. N >= 0.
88: *> \endverbatim
89: *>
90: *> \param[in] NRHS
91: *> \verbatim
92: *> NRHS is INTEGER
93: *> The number of right hand sides, i.e., the number of columns
94: *> of the matrix B. NRHS >= 0.
95: *> \endverbatim
96: *>
97: *> \param[in,out] A
98: *> \verbatim
99: *> A is DOUBLE PRECISION array,
100: *> dimension (LDA,N)
101: *> On entry, the symmetric matrix A. If UPLO = 'U', the leading
102: *> N-by-N upper triangular part of A contains the upper
103: *> triangular part of the matrix A, and the strictly lower
104: *> triangular part of A is not referenced. If UPLO = 'L', the
105: *> leading N-by-N lower triangular part of A contains the lower
106: *> triangular part of the matrix A, and the strictly upper
107: *> triangular part of A is not referenced.
108: *> On exit, if iterative refinement has been successfully used
109: *> (INFO = 0 and ITER >= 0, see description below), then A is
110: *> unchanged, if double precision factorization has been used
111: *> (INFO = 0 and ITER < 0, see description below), then the
112: *> array A contains the factor U or L from the Cholesky
113: *> factorization A = U**T*U or A = L*L**T.
114: *> \endverbatim
115: *>
116: *> \param[in] LDA
117: *> \verbatim
118: *> LDA is INTEGER
119: *> The leading dimension of the array A. LDA >= max(1,N).
120: *> \endverbatim
121: *>
122: *> \param[in] B
123: *> \verbatim
124: *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
125: *> The N-by-NRHS right hand side matrix B.
126: *> \endverbatim
127: *>
128: *> \param[in] LDB
129: *> \verbatim
130: *> LDB is INTEGER
131: *> The leading dimension of the array B. LDB >= max(1,N).
132: *> \endverbatim
133: *>
134: *> \param[out] X
135: *> \verbatim
136: *> X is DOUBLE PRECISION array, dimension (LDX,NRHS)
137: *> If INFO = 0, the N-by-NRHS solution matrix X.
138: *> \endverbatim
139: *>
140: *> \param[in] LDX
141: *> \verbatim
142: *> LDX is INTEGER
143: *> The leading dimension of the array X. LDX >= max(1,N).
144: *> \endverbatim
145: *>
146: *> \param[out] WORK
147: *> \verbatim
148: *> WORK is DOUBLE PRECISION array, dimension (N,NRHS)
149: *> This array is used to hold the residual vectors.
150: *> \endverbatim
151: *>
152: *> \param[out] SWORK
153: *> \verbatim
154: *> SWORK is REAL array, dimension (N*(N+NRHS))
155: *> This array is used to use the single precision matrix and the
156: *> right-hand sides or solutions in single precision.
157: *> \endverbatim
158: *>
159: *> \param[out] ITER
160: *> \verbatim
161: *> ITER is INTEGER
162: *> < 0: iterative refinement has failed, double precision
163: *> factorization has been performed
164: *> -1 : the routine fell back to full precision for
165: *> implementation- or machine-specific reasons
166: *> -2 : narrowing the precision induced an overflow,
167: *> the routine fell back to full precision
168: *> -3 : failure of SPOTRF
169: *> -31: stop the iterative refinement after the 30th
170: *> iterations
171: *> > 0: iterative refinement has been successfully used.
172: *> Returns the number of iterations
173: *> \endverbatim
174: *>
175: *> \param[out] INFO
176: *> \verbatim
177: *> INFO is INTEGER
178: *> = 0: successful exit
179: *> < 0: if INFO = -i, the i-th argument had an illegal value
180: *> > 0: if INFO = i, the leading minor of order i of (DOUBLE
181: *> PRECISION) A is not positive definite, so the
182: *> factorization could not be completed, and the solution
183: *> has not been computed.
184: *> \endverbatim
185: *
186: * Authors:
187: * ========
188: *
189: *> \author Univ. of Tennessee
190: *> \author Univ. of California Berkeley
191: *> \author Univ. of Colorado Denver
192: *> \author NAG Ltd.
193: *
194: *> \ingroup doublePOsolve
195: *
196: * =====================================================================
197: SUBROUTINE DSPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
198: $ SWORK, ITER, INFO )
199: *
200: * -- LAPACK driver routine --
201: * -- LAPACK is a software package provided by Univ. of Tennessee, --
202: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
203: *
204: * .. Scalar Arguments ..
205: CHARACTER UPLO
206: INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
207: * ..
208: * .. Array Arguments ..
209: REAL SWORK( * )
210: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ),
211: $ X( LDX, * )
212: * ..
213: *
214: * =====================================================================
215: *
216: * .. Parameters ..
217: LOGICAL DOITREF
218: PARAMETER ( DOITREF = .TRUE. )
219: *
220: INTEGER ITERMAX
221: PARAMETER ( ITERMAX = 30 )
222: *
223: DOUBLE PRECISION BWDMAX
224: PARAMETER ( BWDMAX = 1.0E+00 )
225: *
226: DOUBLE PRECISION NEGONE, ONE
227: PARAMETER ( NEGONE = -1.0D+0, ONE = 1.0D+0 )
228: *
229: * .. Local Scalars ..
230: INTEGER I, IITER, PTSA, PTSX
231: DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
232: *
233: * .. External Subroutines ..
234: EXTERNAL DAXPY, DSYMM, DLACPY, DLAT2S, DLAG2S, SLAG2D,
235: $ SPOTRF, SPOTRS, DPOTRF, DPOTRS, XERBLA
236: * ..
237: * .. External Functions ..
238: INTEGER IDAMAX
239: DOUBLE PRECISION DLAMCH, DLANSY
240: LOGICAL LSAME
241: EXTERNAL IDAMAX, DLAMCH, DLANSY, LSAME
242: * ..
243: * .. Intrinsic Functions ..
244: INTRINSIC ABS, DBLE, MAX, SQRT
245: * ..
246: * .. Executable Statements ..
247: *
248: INFO = 0
249: ITER = 0
250: *
251: * Test the input parameters.
252: *
253: IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
254: INFO = -1
255: ELSE IF( N.LT.0 ) THEN
256: INFO = -2
257: ELSE IF( NRHS.LT.0 ) THEN
258: INFO = -3
259: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
260: INFO = -5
261: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
262: INFO = -7
263: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
264: INFO = -9
265: END IF
266: IF( INFO.NE.0 ) THEN
267: CALL XERBLA( 'DSPOSV', -INFO )
268: RETURN
269: END IF
270: *
271: * Quick return if (N.EQ.0).
272: *
273: IF( N.EQ.0 )
274: $ RETURN
275: *
276: * Skip single precision iterative refinement if a priori slower
277: * than double precision factorization.
278: *
279: IF( .NOT.DOITREF ) THEN
280: ITER = -1
281: GO TO 40
282: END IF
283: *
284: * Compute some constants.
285: *
286: ANRM = DLANSY( 'I', UPLO, N, A, LDA, WORK )
287: EPS = DLAMCH( 'Epsilon' )
288: CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX
289: *
290: * Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
291: *
292: PTSA = 1
293: PTSX = PTSA + N*N
294: *
295: * Convert B from double precision to single precision and store the
296: * result in SX.
297: *
298: CALL DLAG2S( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO )
299: *
300: IF( INFO.NE.0 ) THEN
301: ITER = -2
302: GO TO 40
303: END IF
304: *
305: * Convert A from double precision to single precision and store the
306: * result in SA.
307: *
308: CALL DLAT2S( UPLO, N, A, LDA, SWORK( PTSA ), N, INFO )
309: *
310: IF( INFO.NE.0 ) THEN
311: ITER = -2
312: GO TO 40
313: END IF
314: *
315: * Compute the Cholesky factorization of SA.
316: *
317: CALL SPOTRF( UPLO, N, SWORK( PTSA ), N, INFO )
318: *
319: IF( INFO.NE.0 ) THEN
320: ITER = -3
321: GO TO 40
322: END IF
323: *
324: * Solve the system SA*SX = SB.
325: *
326: CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
327: $ INFO )
328: *
329: * Convert SX back to double precision
330: *
331: CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
332: *
333: * Compute R = B - AX (R is WORK).
334: *
335: CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
336: *
337: CALL DSYMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
338: $ WORK, N )
339: *
340: * Check whether the NRHS normwise backward errors satisfy the
341: * stopping criterion. If yes, set ITER=0 and return.
342: *
343: DO I = 1, NRHS
344: XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
345: RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
346: IF( RNRM.GT.XNRM*CTE )
347: $ GO TO 10
348: END DO
349: *
350: * If we are here, the NRHS normwise backward errors satisfy the
351: * stopping criterion. We are good to exit.
352: *
353: ITER = 0
354: RETURN
355: *
356: 10 CONTINUE
357: *
358: DO 30 IITER = 1, ITERMAX
359: *
360: * Convert R (in WORK) from double precision to single precision
361: * and store the result in SX.
362: *
363: CALL DLAG2S( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO )
364: *
365: IF( INFO.NE.0 ) THEN
366: ITER = -2
367: GO TO 40
368: END IF
369: *
370: * Solve the system SA*SX = SR.
371: *
372: CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
373: $ INFO )
374: *
375: * Convert SX back to double precision and update the current
376: * iterate.
377: *
378: CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
379: *
380: DO I = 1, NRHS
381: CALL DAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 )
382: END DO
383: *
384: * Compute R = B - AX (R is WORK).
385: *
386: CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
387: *
388: CALL DSYMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
389: $ WORK, N )
390: *
391: * Check whether the NRHS normwise backward errors satisfy the
392: * stopping criterion. If yes, set ITER=IITER>0 and return.
393: *
394: DO I = 1, NRHS
395: XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
396: RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
397: IF( RNRM.GT.XNRM*CTE )
398: $ GO TO 20
399: END DO
400: *
401: * If we are here, the NRHS normwise backward errors satisfy the
402: * stopping criterion, we are good to exit.
403: *
404: ITER = IITER
405: *
406: RETURN
407: *
408: 20 CONTINUE
409: *
410: 30 CONTINUE
411: *
412: * If we are at this place of the code, this is because we have
413: * performed ITER=ITERMAX iterations and never satisfied the
414: * stopping criterion, set up the ITER flag accordingly and follow
415: * up on double precision routine.
416: *
417: ITER = -ITERMAX - 1
418: *
419: 40 CONTINUE
420: *
421: * Single-precision iterative refinement failed to converge to a
422: * satisfactory solution, so we resort to double precision.
423: *
424: CALL DPOTRF( UPLO, N, A, LDA, INFO )
425: *
426: IF( INFO.NE.0 )
427: $ RETURN
428: *
429: CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX )
430: CALL DPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO )
431: *
432: RETURN
433: *
434: * End of DSPOSV
435: *
436: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsposv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack