File:
[local] /
rpl /
lapack /
lapack /
dsysvx.f
Revision
1.14:
download - view:
text,
annotated -
select for diffs -
revision graph
Mon Jan 27 09:28:28 2014 UTC (10 years, 7 months ago) by
bertrand
Branches:
MAIN
CVS tags:
rpl-4_1_24,
rpl-4_1_23,
rpl-4_1_22,
rpl-4_1_21,
rpl-4_1_20,
rpl-4_1_19,
rpl-4_1_18,
rpl-4_1_17,
HEAD
Cohérence.
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: *> \date April 2012
279: *
280: *> \ingroup doubleSYsolve
281: *
282: * =====================================================================
283: SUBROUTINE DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
284: $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
285: $ IWORK, INFO )
286: *
287: * -- LAPACK driver routine (version 3.4.1) --
288: * -- LAPACK is a software package provided by Univ. of Tennessee, --
289: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
290: * April 2012
291: *
292: * .. Scalar Arguments ..
293: CHARACTER FACT, UPLO
294: INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
295: DOUBLE PRECISION RCOND
296: * ..
297: * .. Array Arguments ..
298: INTEGER IPIV( * ), IWORK( * )
299: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
300: $ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
301: * ..
302: *
303: * =====================================================================
304: *
305: * .. Parameters ..
306: DOUBLE PRECISION ZERO
307: PARAMETER ( ZERO = 0.0D+0 )
308: * ..
309: * .. Local Scalars ..
310: LOGICAL LQUERY, NOFACT
311: INTEGER LWKOPT, NB
312: DOUBLE PRECISION ANORM
313: * ..
314: * .. External Functions ..
315: LOGICAL LSAME
316: INTEGER ILAENV
317: DOUBLE PRECISION DLAMCH, DLANSY
318: EXTERNAL LSAME, ILAENV, DLAMCH, DLANSY
319: * ..
320: * .. External Subroutines ..
321: EXTERNAL DLACPY, DSYCON, DSYRFS, DSYTRF, DSYTRS, XERBLA
322: * ..
323: * .. Intrinsic Functions ..
324: INTRINSIC MAX
325: * ..
326: * .. Executable Statements ..
327: *
328: * Test the input parameters.
329: *
330: INFO = 0
331: NOFACT = LSAME( FACT, 'N' )
332: LQUERY = ( LWORK.EQ.-1 )
333: IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
334: INFO = -1
335: ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
336: $ THEN
337: INFO = -2
338: ELSE IF( N.LT.0 ) THEN
339: INFO = -3
340: ELSE IF( NRHS.LT.0 ) THEN
341: INFO = -4
342: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
343: INFO = -6
344: ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
345: INFO = -8
346: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
347: INFO = -11
348: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
349: INFO = -13
350: ELSE IF( LWORK.LT.MAX( 1, 3*N ) .AND. .NOT.LQUERY ) THEN
351: INFO = -18
352: END IF
353: *
354: IF( INFO.EQ.0 ) THEN
355: LWKOPT = MAX( 1, 3*N )
356: IF( NOFACT ) THEN
357: NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 )
358: LWKOPT = MAX( LWKOPT, N*NB )
359: END IF
360: WORK( 1 ) = LWKOPT
361: END IF
362: *
363: IF( INFO.NE.0 ) THEN
364: CALL XERBLA( 'DSYSVX', -INFO )
365: RETURN
366: ELSE IF( LQUERY ) THEN
367: RETURN
368: END IF
369: *
370: IF( NOFACT ) THEN
371: *
372: * Compute the factorization A = U*D*U**T or A = L*D*L**T.
373: *
374: CALL DLACPY( UPLO, N, N, A, LDA, AF, LDAF )
375: CALL DSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO )
376: *
377: * Return if INFO is non-zero.
378: *
379: IF( INFO.GT.0 )THEN
380: RCOND = ZERO
381: RETURN
382: END IF
383: END IF
384: *
385: * Compute the norm of the matrix A.
386: *
387: ANORM = DLANSY( 'I', UPLO, N, A, LDA, WORK )
388: *
389: * Compute the reciprocal of the condition number of A.
390: *
391: CALL DSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, IWORK,
392: $ INFO )
393: *
394: * Compute the solution vectors X.
395: *
396: CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
397: CALL DSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
398: *
399: * Use iterative refinement to improve the computed solutions and
400: * compute error bounds and backward error estimates for them.
401: *
402: CALL DSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X,
403: $ LDX, FERR, BERR, WORK, IWORK, INFO )
404: *
405: * Set INFO = N+1 if the matrix is singular to working precision.
406: *
407: IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
408: $ INFO = N + 1
409: *
410: WORK( 1 ) = LWKOPT
411: *
412: RETURN
413: *
414: * End of DSYSVX
415: *
416: END
CVSweb interface <joel.bertrand@systella.fr>