1*> \brief <b> DSPSVX computes the solution to system of linear equations A * X = B for OTHER 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 DSPSVX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dspsvx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dspsvx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dspsvx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DSPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X,
22*                          LDX, RCOND, FERR, BERR, WORK, IWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          FACT, UPLO
26*       INTEGER            INFO, LDB, LDX, N, NRHS
27*       DOUBLE PRECISION   RCOND
28*       ..
29*       .. Array Arguments ..
30*       INTEGER            IPIV( * ), IWORK( * )
31*       DOUBLE PRECISION   AFP( * ), AP( * ), B( LDB, * ), BERR( * ),
32*      $                   FERR( * ), WORK( * ), X( LDX, * )
33*       ..
34*
35*
36*> \par Purpose:
37*  =============
38*>
39*> \verbatim
40*>
41*> DSPSVX uses the diagonal pivoting factorization A = U*D*U**T or
42*> A = L*D*L**T to compute the solution to a real system of linear
43*> equations A * X = B, where A is an N-by-N symmetric matrix stored
44*> in packed format and X and B are N-by-NRHS matrices.
45*>
46*> Error bounds on the solution and a condition estimate are also
47*> provided.
48*> \endverbatim
49*
50*> \par Description:
51*  =================
52*>
53*> \verbatim
54*>
55*> The following steps are performed:
56*>
57*> 1. If FACT = 'N', the diagonal pivoting method is used to factor A as
58*>       A = U * D * U**T,  if UPLO = 'U', or
59*>       A = L * D * L**T,  if UPLO = 'L',
60*>    where U (or L) is a product of permutation and unit upper (lower)
61*>    triangular matrices and D is symmetric and block diagonal with
62*>    1-by-1 and 2-by-2 diagonal blocks.
63*>
64*> 2. If some D(i,i)=0, so that D is exactly singular, then the routine
65*>    returns with INFO = i. Otherwise, the factored form of A is used
66*>    to estimate the condition number of the matrix A.  If the
67*>    reciprocal of the condition number is less than machine precision,
68*>    INFO = N+1 is returned as a warning, but the routine still goes on
69*>    to solve for X and compute error bounds as described below.
70*>
71*> 3. The system of equations is solved for X using the factored form
72*>    of A.
73*>
74*> 4. Iterative refinement is applied to improve the computed solution
75*>    matrix and calculate error bounds and backward error estimates
76*>    for it.
77*> \endverbatim
78*
79*  Arguments:
80*  ==========
81*
82*> \param[in] FACT
83*> \verbatim
84*>          FACT is CHARACTER*1
85*>          Specifies whether or not the factored form of A has been
86*>          supplied on entry.
87*>          = 'F':  On entry, AFP and IPIV contain the factored form of
88*>                  A.  AP, AFP and IPIV will not be modified.
89*>          = 'N':  The matrix A will be copied to AFP and factored.
90*> \endverbatim
91*>
92*> \param[in] UPLO
93*> \verbatim
94*>          UPLO is CHARACTER*1
95*>          = 'U':  Upper triangle of A is stored;
96*>          = 'L':  Lower triangle of A is stored.
97*> \endverbatim
98*>
99*> \param[in] N
100*> \verbatim
101*>          N is INTEGER
102*>          The number of linear equations, i.e., the order of the
103*>          matrix A.  N >= 0.
104*> \endverbatim
105*>
106*> \param[in] NRHS
107*> \verbatim
108*>          NRHS is INTEGER
109*>          The number of right hand sides, i.e., the number of columns
110*>          of the matrices B and X.  NRHS >= 0.
111*> \endverbatim
112*>
113*> \param[in] AP
114*> \verbatim
115*>          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
116*>          The upper or lower triangle of the symmetric matrix A, packed
117*>          columnwise in a linear array.  The j-th column of A is stored
118*>          in the array AP as follows:
119*>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
120*>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
121*>          See below for further details.
122*> \endverbatim
123*>
124*> \param[in,out] AFP
125*> \verbatim
126*>          AFP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
127*>          If FACT = 'F', then AFP is an input argument and on entry
128*>          contains the block diagonal matrix D and the multipliers used
129*>          to obtain the factor U or L from the factorization
130*>          A = U*D*U**T or A = L*D*L**T as computed by DSPTRF, stored as
131*>          a packed triangular matrix in the same storage format as A.
132*>
133*>          If FACT = 'N', then AFP is an output argument and on exit
134*>          contains the block diagonal matrix D and the multipliers used
135*>          to obtain the factor U or L from the factorization
136*>          A = U*D*U**T or A = L*D*L**T as computed by DSPTRF, stored as
137*>          a packed triangular matrix in the same storage format as A.
138*> \endverbatim
139*>
140*> \param[in,out] IPIV
141*> \verbatim
142*>          IPIV is INTEGER array, dimension (N)
143*>          If FACT = 'F', then IPIV is an input argument and on entry
144*>          contains details of the interchanges and the block structure
145*>          of D, as determined by DSPTRF.
146*>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
147*>          interchanged and D(k,k) is a 1-by-1 diagonal block.
148*>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
149*>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
150*>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
151*>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
152*>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
153*>
154*>          If FACT = 'N', then IPIV is an output argument and on exit
155*>          contains details of the interchanges and the block structure
156*>          of D, as determined by DSPTRF.
157*> \endverbatim
158*>
159*> \param[in] B
160*> \verbatim
161*>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
162*>          The N-by-NRHS right hand side matrix B.
163*> \endverbatim
164*>
165*> \param[in] LDB
166*> \verbatim
167*>          LDB is INTEGER
168*>          The leading dimension of the array B.  LDB >= max(1,N).
169*> \endverbatim
170*>
171*> \param[out] X
172*> \verbatim
173*>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
174*>          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
175*> \endverbatim
176*>
177*> \param[in] LDX
178*> \verbatim
179*>          LDX is INTEGER
180*>          The leading dimension of the array X.  LDX >= max(1,N).
181*> \endverbatim
182*>
183*> \param[out] RCOND
184*> \verbatim
185*>          RCOND is DOUBLE PRECISION
186*>          The estimate of the reciprocal condition number of the matrix
187*>          A.  If RCOND is less than the machine precision (in
188*>          particular, if RCOND = 0), the matrix is singular to working
189*>          precision.  This condition is indicated by a return code of
190*>          INFO > 0.
191*> \endverbatim
192*>
193*> \param[out] FERR
194*> \verbatim
195*>          FERR is DOUBLE PRECISION array, dimension (NRHS)
196*>          The estimated forward error bound for each solution vector
197*>          X(j) (the j-th column of the solution matrix X).
198*>          If XTRUE is the true solution corresponding to X(j), FERR(j)
199*>          is an estimated upper bound for the magnitude of the largest
200*>          element in (X(j) - XTRUE) divided by the magnitude of the
201*>          largest element in X(j).  The estimate is as reliable as
202*>          the estimate for RCOND, and is almost always a slight
203*>          overestimate of the true error.
204*> \endverbatim
205*>
206*> \param[out] BERR
207*> \verbatim
208*>          BERR is DOUBLE PRECISION array, dimension (NRHS)
209*>          The componentwise relative backward error of each solution
210*>          vector X(j) (i.e., the smallest relative change in
211*>          any element of A or B that makes X(j) an exact solution).
212*> \endverbatim
213*>
214*> \param[out] WORK
215*> \verbatim
216*>          WORK is DOUBLE PRECISION array, dimension (3*N)
217*> \endverbatim
218*>
219*> \param[out] IWORK
220*> \verbatim
221*>          IWORK is INTEGER array, dimension (N)
222*> \endverbatim
223*>
224*> \param[out] INFO
225*> \verbatim
226*>          INFO is INTEGER
227*>          = 0: successful exit
228*>          < 0: if INFO = -i, the i-th argument had an illegal value
229*>          > 0:  if INFO = i, and i is
230*>                <= N:  D(i,i) is exactly zero.  The factorization
231*>                       has been completed but the factor D is exactly
232*>                       singular, so the solution and error bounds could
233*>                       not be computed. RCOND = 0 is returned.
234*>                = N+1: D is nonsingular, but RCOND is less than machine
235*>                       precision, meaning that the matrix is singular
236*>                       to working precision.  Nevertheless, the
237*>                       solution and error bounds are computed because
238*>                       there are a number of situations where the
239*>                       computed solution can be more accurate than the
240*>                       value of RCOND would suggest.
241*> \endverbatim
242*
243*  Authors:
244*  ========
245*
246*> \author Univ. of Tennessee
247*> \author Univ. of California Berkeley
248*> \author Univ. of Colorado Denver
249*> \author NAG Ltd.
250*
251*> \ingroup doubleOTHERsolve
252*
253*> \par Further Details:
254*  =====================
255*>
256*> \verbatim
257*>
258*>  The packed storage scheme is illustrated by the following example
259*>  when N = 4, UPLO = 'U':
260*>
261*>  Two-dimensional storage of the symmetric matrix A:
262*>
263*>     a11 a12 a13 a14
264*>         a22 a23 a24
265*>             a33 a34     (aij = aji)
266*>                 a44
267*>
268*>  Packed storage of the upper triangle of A:
269*>
270*>  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
271*> \endverbatim
272*>
273*  =====================================================================
274      SUBROUTINE DSPSVX( FACT, UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X,
275     $                   LDX, RCOND, FERR, BERR, WORK, IWORK, INFO )
276*
277*  -- LAPACK driver routine --
278*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
279*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
280*
281*     .. Scalar Arguments ..
282      CHARACTER          FACT, UPLO
283      INTEGER            INFO, LDB, LDX, N, NRHS
284      DOUBLE PRECISION   RCOND
285*     ..
286*     .. Array Arguments ..
287      INTEGER            IPIV( * ), IWORK( * )
288      DOUBLE PRECISION   AFP( * ), AP( * ), B( LDB, * ), BERR( * ),
289     $                   FERR( * ), WORK( * ), X( LDX, * )
290*     ..
291*
292*  =====================================================================
293*
294*     .. Parameters ..
295      DOUBLE PRECISION   ZERO
296      PARAMETER          ( ZERO = 0.0D+0 )
297*     ..
298*     .. Local Scalars ..
299      LOGICAL            NOFACT
300      DOUBLE PRECISION   ANORM
301*     ..
302*     .. External Functions ..
303      LOGICAL            LSAME
304      DOUBLE PRECISION   DLAMCH, DLANSP
305      EXTERNAL           LSAME, DLAMCH, DLANSP
306*     ..
307*     .. External Subroutines ..
308      EXTERNAL           DCOPY, DLACPY, DSPCON, DSPRFS, DSPTRF, DSPTRS,
309     $                   XERBLA
310*     ..
311*     .. Intrinsic Functions ..
312      INTRINSIC          MAX
313*     ..
314*     .. Executable Statements ..
315*
316*     Test the input parameters.
317*
318      INFO = 0
319      NOFACT = LSAME( FACT, 'N' )
320      IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
321         INFO = -1
322      ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
323     $          THEN
324         INFO = -2
325      ELSE IF( N.LT.0 ) THEN
326         INFO = -3
327      ELSE IF( NRHS.LT.0 ) THEN
328         INFO = -4
329      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
330         INFO = -9
331      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
332         INFO = -11
333      END IF
334      IF( INFO.NE.0 ) THEN
335         CALL XERBLA( 'DSPSVX', -INFO )
336         RETURN
337      END IF
338*
339      IF( NOFACT ) THEN
340*
341*        Compute the factorization A = U*D*U**T or A = L*D*L**T.
342*
343         CALL DCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 )
344         CALL DSPTRF( UPLO, N, AFP, IPIV, INFO )
345*
346*        Return if INFO is non-zero.
347*
348         IF( INFO.GT.0 )THEN
349            RCOND = ZERO
350            RETURN
351         END IF
352      END IF
353*
354*     Compute the norm of the matrix A.
355*
356      ANORM = DLANSP( 'I', UPLO, N, AP, WORK )
357*
358*     Compute the reciprocal of the condition number of A.
359*
360      CALL DSPCON( UPLO, N, AFP, IPIV, ANORM, RCOND, WORK, IWORK, INFO )
361*
362*     Compute the solution vectors X.
363*
364      CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
365      CALL DSPTRS( UPLO, N, NRHS, AFP, IPIV, X, LDX, INFO )
366*
367*     Use iterative refinement to improve the computed solutions and
368*     compute error bounds and backward error estimates for them.
369*
370      CALL DSPRFS( UPLO, N, NRHS, AP, AFP, IPIV, B, LDB, X, LDX, FERR,
371     $             BERR, WORK, IWORK, INFO )
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 DSPSVX
381*
382      END
383