LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsposv.f
Go to the documentation of this file.
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 principal minor of order i
181*> of (DOUBLE PRECISION) A is not positive, 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 posv_mixed
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,
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
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dlag2s(m, n, a, lda, sa, ldsa, info)
DLAG2S converts a double precision matrix to a single precision matrix.
Definition dlag2s.f:108
subroutine slag2d(m, n, sa, ldsa, a, lda, info)
SLAG2D converts a single precision matrix to a double precision matrix.
Definition slag2d.f:104
subroutine dlat2s(uplo, n, a, lda, sa, ldsa, info)
DLAT2S converts a double-precision triangular matrix to a single-precision triangular matrix.
Definition dlat2s.f:111
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dsymm(side, uplo, m, n, alpha, a, lda, b, ldb, beta, c, ldc)
DSYMM
Definition dsymm.f:189
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
subroutine dsposv(uplo, n, nrhs, a, lda, b, ldb, x, ldx, work, swork, iter, info)
DSPOSV computes the solution to system of linear equations A * X = B for PO matrices
Definition dsposv.f:199
subroutine dpotrf(uplo, n, a, lda, info)
DPOTRF
Definition dpotrf.f:107
subroutine spotrf(uplo, n, a, lda, info)
SPOTRF
Definition spotrf.f:107
subroutine dpotrs(uplo, n, nrhs, a, lda, b, ldb, info)
DPOTRS
Definition dpotrs.f:110
subroutine spotrs(uplo, n, nrhs, a, lda, b, ldb, info)
SPOTRS
Definition spotrs.f:110