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.EQ.0 and ITER.GE.0, see description below), then A is
110*>          unchanged, if double precision factorization has been used
111*>          (INFO.EQ.0 and ITER.LT.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 sucessfully 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*> \date November 2011
195*
196*> \ingroup doublePOsolve
197*
198*  =====================================================================
199      SUBROUTINE DSPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
200     $                   SWORK, ITER, INFO )
201*
202*  -- LAPACK driver routine (version 3.4.0) --
203*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
204*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
205*     November 2011
206*
207*     .. Scalar Arguments ..
208      CHARACTER          UPLO
209      INTEGER            INFO, ITER, LDA, LDB, LDX, N, NRHS
210*     ..
211*     .. Array Arguments ..
212      REAL               SWORK( * )
213      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( N, * ),
214     $                   X( LDX, * )
215*     ..
216*
217*  =====================================================================
218*
219*     .. Parameters ..
220      LOGICAL            DOITREF
221      PARAMETER          ( DOITREF = .TRUE. )
222*
223      INTEGER            ITERMAX
224      PARAMETER          ( ITERMAX = 30 )
225*
226      DOUBLE PRECISION   BWDMAX
227      PARAMETER          ( BWDMAX = 1.0E+00 )
228*
229      DOUBLE PRECISION   NEGONE, ONE
230      PARAMETER          ( NEGONE = -1.0D+0, ONE = 1.0D+0 )
231*
232*     .. Local Scalars ..
233      INTEGER            I, IITER, PTSA, PTSX
234      DOUBLE PRECISION   ANRM, CTE, EPS, RNRM, XNRM
235*
236*     .. External Subroutines ..
237      EXTERNAL           DAXPY, DSYMM, DLACPY, DLAT2S, DLAG2S, SLAG2D,
238     $                   SPOTRF, SPOTRS, XERBLA
239*     ..
240*     .. External Functions ..
241      INTEGER            IDAMAX
242      DOUBLE PRECISION   DLAMCH, DLANSY
243      LOGICAL            LSAME
244      EXTERNAL           IDAMAX, DLAMCH, DLANSY, LSAME
245*     ..
246*     .. Intrinsic Functions ..
247      INTRINSIC          ABS, DBLE, MAX, SQRT
248*     ..
249*     .. Executable Statements ..
250*
251      INFO = 0
252      ITER = 0
253*
254*     Test the input parameters.
255*
256      IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
257         INFO = -1
258      ELSE IF( N.LT.0 ) THEN
259         INFO = -2
260      ELSE IF( NRHS.LT.0 ) THEN
261         INFO = -3
262      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
263         INFO = -5
264      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
265         INFO = -7
266      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
267         INFO = -9
268      END IF
269      IF( INFO.NE.0 ) THEN
270         CALL XERBLA( 'DSPOSV', -INFO )
271         RETURN
272      END IF
273*
274*     Quick return if (N.EQ.0).
275*
276      IF( N.EQ.0 )
277     $   RETURN
278*
279*     Skip single precision iterative refinement if a priori slower
280*     than double precision factorization.
281*
282      IF( .NOT.DOITREF ) THEN
283         ITER = -1
284         GO TO 40
285      END IF
286*
287*     Compute some constants.
288*
289      ANRM = DLANSY( 'I', UPLO, N, A, LDA, WORK )
290      EPS = DLAMCH( 'Epsilon' )
291      CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX
292*
293*     Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
294*
295      PTSA = 1
296      PTSX = PTSA + N*N
297*
298*     Convert B from double precision to single precision and store the
299*     result in SX.
300*
301      CALL DLAG2S( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO )
302*
303      IF( INFO.NE.0 ) THEN
304         ITER = -2
305         GO TO 40
306      END IF
307*
308*     Convert A from double precision to single precision and store the
309*     result in SA.
310*
311      CALL DLAT2S( UPLO, N, A, LDA, SWORK( PTSA ), N, INFO )
312*
313      IF( INFO.NE.0 ) THEN
314         ITER = -2
315         GO TO 40
316      END IF
317*
318*     Compute the Cholesky factorization of SA.
319*
320      CALL SPOTRF( UPLO, N, SWORK( PTSA ), N, INFO )
321*
322      IF( INFO.NE.0 ) THEN
323         ITER = -3
324         GO TO 40
325      END IF
326*
327*     Solve the system SA*SX = SB.
328*
329      CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
330     $             INFO )
331*
332*     Convert SX back to double precision
333*
334      CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
335*
336*     Compute R = B - AX (R is WORK).
337*
338      CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
339*
340      CALL DSYMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
341     $            WORK, N )
342*
343*     Check whether the NRHS normwise backward errors satisfy the
344*     stopping criterion. If yes, set ITER=0 and return.
345*
346      DO I = 1, NRHS
347         XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
348         RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
349         IF( RNRM.GT.XNRM*CTE )
350     $      GO TO 10
351      END DO
352*
353*     If we are here, the NRHS normwise backward errors satisfy the
354*     stopping criterion. We are good to exit.
355*
356      ITER = 0
357      RETURN
358*
359   10 CONTINUE
360*
361      DO 30 IITER = 1, ITERMAX
362*
363*        Convert R (in WORK) from double precision to single precision
364*        and store the result in SX.
365*
366         CALL DLAG2S( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO )
367*
368         IF( INFO.NE.0 ) THEN
369            ITER = -2
370            GO TO 40
371         END IF
372*
373*        Solve the system SA*SX = SR.
374*
375         CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N,
376     $                INFO )
377*
378*        Convert SX back to double precision and update the current
379*        iterate.
380*
381         CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
382*
383         DO I = 1, NRHS
384            CALL DAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 )
385         END DO
386*
387*        Compute R = B - AX (R is WORK).
388*
389         CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N )
390*
391         CALL DSYMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE,
392     $               WORK, N )
393*
394*        Check whether the NRHS normwise backward errors satisfy the
395*        stopping criterion. If yes, set ITER=IITER>0 and return.
396*
397         DO I = 1, NRHS
398            XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) )
399            RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) )
400            IF( RNRM.GT.XNRM*CTE )
401     $         GO TO 20
402         END DO
403*
404*        If we are here, the NRHS normwise backward errors satisfy the
405*        stopping criterion, we are good to exit.
406*
407         ITER = IITER
408*
409         RETURN
410*
411   20    CONTINUE
412*
413   30 CONTINUE
414*
415*     If we are at this place of the code, this is because we have
416*     performed ITER=ITERMAX iterations and never satisified the
417*     stopping criterion, set up the ITER flag accordingly and follow
418*     up on double precision routine.
419*
420      ITER = -ITERMAX - 1
421*
422   40 CONTINUE
423*
424*     Single-precision iterative refinement failed to converge to a
425*     satisfactory solution, so we resort to double precision.
426*
427      CALL DPOTRF( UPLO, N, A, LDA, INFO )
428*
429      IF( INFO.NE.0 )
430     $   RETURN
431*
432      CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX )
433      CALL DPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO )
434*
435      RETURN
436*
437*     End of DSPOSV.
438*
439      END
440