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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsysvx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack