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