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