Annotation of rpl/lapack/lapack/zsysvx.f, revision 1.1.1.1
1.1 bertrand 1: SUBROUTINE ZSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
2: $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
3: $ RWORK, INFO )
4: *
5: * -- LAPACK driver routine (version 3.2) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * November 2006
9: *
10: * .. Scalar Arguments ..
11: CHARACTER FACT, UPLO
12: INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
13: DOUBLE PRECISION RCOND
14: * ..
15: * .. Array Arguments ..
16: INTEGER IPIV( * )
17: DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
18: COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
19: $ WORK( * ), X( LDX, * )
20: * ..
21: *
22: * Purpose
23: * =======
24: *
25: * ZSYSVX uses the diagonal pivoting factorization to compute the
26: * solution to a complex system of linear equations A * X = B,
27: * where A is an N-by-N symmetric matrix and X and B are N-by-NRHS
28: * matrices.
29: *
30: * Error bounds on the solution and a condition estimate are also
31: * provided.
32: *
33: * Description
34: * ===========
35: *
36: * The following steps are performed:
37: *
38: * 1. If FACT = 'N', the diagonal pivoting method is used to factor A.
39: * The form of the factorization is
40: * A = U * D * U**T, if UPLO = 'U', or
41: * A = L * D * L**T, if UPLO = 'L',
42: * where U (or L) is a product of permutation and unit upper (lower)
43: * triangular matrices, and D is symmetric and block diagonal with
44: * 1-by-1 and 2-by-2 diagonal blocks.
45: *
46: * 2. If some D(i,i)=0, so that D is exactly singular, then the routine
47: * returns with INFO = i. Otherwise, the factored form of A is used
48: * to estimate the condition number of the matrix A. If the
49: * reciprocal of the condition number is less than machine precision,
50: * INFO = N+1 is returned as a warning, but the routine still goes on
51: * to solve for X and compute error bounds as described below.
52: *
53: * 3. The system of equations is solved for X using the factored form
54: * of A.
55: *
56: * 4. Iterative refinement is applied to improve the computed solution
57: * matrix and calculate error bounds and backward error estimates
58: * for it.
59: *
60: * Arguments
61: * =========
62: *
63: * FACT (input) CHARACTER*1
64: * Specifies whether or not the factored form of A has been
65: * supplied on entry.
66: * = 'F': On entry, AF and IPIV contain the factored form
67: * of A. A, AF and IPIV will not be modified.
68: * = 'N': The matrix A will be copied to AF and factored.
69: *
70: * UPLO (input) CHARACTER*1
71: * = 'U': Upper triangle of A is stored;
72: * = 'L': Lower triangle of A is stored.
73: *
74: * N (input) INTEGER
75: * The number of linear equations, i.e., the order of the
76: * matrix A. N >= 0.
77: *
78: * NRHS (input) INTEGER
79: * The number of right hand sides, i.e., the number of columns
80: * of the matrices B and X. NRHS >= 0.
81: *
82: * A (input) COMPLEX*16 array, dimension (LDA,N)
83: * The symmetric matrix A. If UPLO = 'U', the leading N-by-N
84: * upper triangular part of A contains the upper triangular part
85: * of the matrix A, and the strictly lower triangular part of A
86: * is not referenced. If UPLO = 'L', the leading N-by-N lower
87: * triangular part of A contains the lower triangular part of
88: * the matrix A, and the strictly upper triangular part of A is
89: * not referenced.
90: *
91: * LDA (input) INTEGER
92: * The leading dimension of the array A. LDA >= max(1,N).
93: *
94: * AF (input or output) COMPLEX*16 array, dimension (LDAF,N)
95: * If FACT = 'F', then AF is an input argument and on entry
96: * contains the block diagonal matrix D and the multipliers used
97: * to obtain the factor U or L from the factorization
98: * A = U*D*U**T or A = L*D*L**T as computed by ZSYTRF.
99: *
100: * If FACT = 'N', then AF is an output argument and on exit
101: * returns the block diagonal matrix D and the multipliers used
102: * to obtain the factor U or L from the factorization
103: * A = U*D*U**T or A = L*D*L**T.
104: *
105: * LDAF (input) INTEGER
106: * The leading dimension of the array AF. LDAF >= max(1,N).
107: *
108: * IPIV (input or output) INTEGER array, dimension (N)
109: * If FACT = 'F', then IPIV is an input argument and on entry
110: * contains details of the interchanges and the block structure
111: * of D, as determined by ZSYTRF.
112: * If IPIV(k) > 0, then rows and columns k and IPIV(k) were
113: * interchanged and D(k,k) is a 1-by-1 diagonal block.
114: * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
115: * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
116: * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
117: * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
118: * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
119: *
120: * If FACT = 'N', then IPIV is an output argument and on exit
121: * contains details of the interchanges and the block structure
122: * of D, as determined by ZSYTRF.
123: *
124: * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
125: * The N-by-NRHS right hand side matrix B.
126: *
127: * LDB (input) INTEGER
128: * The leading dimension of the array B. LDB >= max(1,N).
129: *
130: * X (output) COMPLEX*16 array, dimension (LDX,NRHS)
131: * If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
132: *
133: * LDX (input) INTEGER
134: * The leading dimension of the array X. LDX >= max(1,N).
135: *
136: * RCOND (output) DOUBLE PRECISION
137: * The estimate of the reciprocal condition number of the matrix
138: * A. If RCOND is less than the machine precision (in
139: * particular, if RCOND = 0), the matrix is singular to working
140: * precision. This condition is indicated by a return code of
141: * INFO > 0.
142: *
143: * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
144: * The estimated forward error bound for each solution vector
145: * X(j) (the j-th column of the solution matrix X).
146: * If XTRUE is the true solution corresponding to X(j), FERR(j)
147: * is an estimated upper bound for the magnitude of the largest
148: * element in (X(j) - XTRUE) divided by the magnitude of the
149: * largest element in X(j). The estimate is as reliable as
150: * the estimate for RCOND, and is almost always a slight
151: * overestimate of the true error.
152: *
153: * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
154: * The componentwise relative backward error of each solution
155: * vector X(j) (i.e., the smallest relative change in
156: * any element of A or B that makes X(j) an exact solution).
157: *
158: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
159: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
160: *
161: * LWORK (input) INTEGER
162: * The length of WORK. LWORK >= max(1,2*N), and for best
163: * performance, when FACT = 'N', LWORK >= max(1,2*N,N*NB), where
164: * NB is the optimal blocksize for ZSYTRF.
165: *
166: * If LWORK = -1, then a workspace query is assumed; the routine
167: * only calculates the optimal size of the WORK array, returns
168: * this value as the first entry of the WORK array, and no error
169: * message related to LWORK is issued by XERBLA.
170: *
171: * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
172: *
173: * INFO (output) INTEGER
174: * = 0: successful exit
175: * < 0: if INFO = -i, the i-th argument had an illegal value
176: * > 0: if INFO = i, and i is
177: * <= N: D(i,i) is exactly zero. The factorization
178: * has been completed but the factor D is exactly
179: * singular, so the solution and error bounds could
180: * not be computed. RCOND = 0 is returned.
181: * = N+1: D is nonsingular, but RCOND is less than machine
182: * precision, meaning that the matrix is singular
183: * to working precision. Nevertheless, the
184: * solution and error bounds are computed because
185: * there are a number of situations where the
186: * computed solution can be more accurate than the
187: * value of RCOND would suggest.
188: *
189: * =====================================================================
190: *
191: * .. Parameters ..
192: DOUBLE PRECISION ZERO
193: PARAMETER ( ZERO = 0.0D+0 )
194: * ..
195: * .. Local Scalars ..
196: LOGICAL LQUERY, NOFACT
197: INTEGER LWKOPT, NB
198: DOUBLE PRECISION ANORM
199: * ..
200: * .. External Functions ..
201: LOGICAL LSAME
202: INTEGER ILAENV
203: DOUBLE PRECISION DLAMCH, ZLANSY
204: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANSY
205: * ..
206: * .. External Subroutines ..
207: EXTERNAL XERBLA, ZLACPY, ZSYCON, ZSYRFS, ZSYTRF, ZSYTRS
208: * ..
209: * .. Intrinsic Functions ..
210: INTRINSIC MAX
211: * ..
212: * .. Executable Statements ..
213: *
214: * Test the input parameters.
215: *
216: INFO = 0
217: NOFACT = LSAME( FACT, 'N' )
218: LQUERY = ( LWORK.EQ.-1 )
219: IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
220: INFO = -1
221: ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
222: $ THEN
223: INFO = -2
224: ELSE IF( N.LT.0 ) THEN
225: INFO = -3
226: ELSE IF( NRHS.LT.0 ) THEN
227: INFO = -4
228: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
229: INFO = -6
230: ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
231: INFO = -8
232: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
233: INFO = -11
234: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
235: INFO = -13
236: ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
237: INFO = -18
238: END IF
239: *
240: IF( INFO.EQ.0 ) THEN
241: LWKOPT = MAX( 1, 2*N )
242: IF( NOFACT ) THEN
243: NB = ILAENV( 1, 'ZSYTRF', UPLO, N, -1, -1, -1 )
244: LWKOPT = MAX( LWKOPT, N*NB )
245: END IF
246: WORK( 1 ) = LWKOPT
247: END IF
248: *
249: IF( INFO.NE.0 ) THEN
250: CALL XERBLA( 'ZSYSVX', -INFO )
251: RETURN
252: ELSE IF( LQUERY ) THEN
253: RETURN
254: END IF
255: *
256: IF( NOFACT ) THEN
257: *
258: * Compute the factorization A = U*D*U' or A = L*D*L'.
259: *
260: CALL ZLACPY( UPLO, N, N, A, LDA, AF, LDAF )
261: CALL ZSYTRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO )
262: *
263: * Return if INFO is non-zero.
264: *
265: IF( INFO.GT.0 )THEN
266: RCOND = ZERO
267: RETURN
268: END IF
269: END IF
270: *
271: * Compute the norm of the matrix A.
272: *
273: ANORM = ZLANSY( 'I', UPLO, N, A, LDA, RWORK )
274: *
275: * Compute the reciprocal of the condition number of A.
276: *
277: CALL ZSYCON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, INFO )
278: *
279: * Compute the solution vectors X.
280: *
281: CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
282: CALL ZSYTRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
283: *
284: * Use iterative refinement to improve the computed solutions and
285: * compute error bounds and backward error estimates for them.
286: *
287: CALL ZSYRFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X,
288: $ LDX, FERR, BERR, WORK, RWORK, INFO )
289: *
290: * Set INFO = N+1 if the matrix is singular to working precision.
291: *
292: IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
293: $ INFO = N + 1
294: *
295: WORK( 1 ) = LWKOPT
296: *
297: RETURN
298: *
299: * End of ZSYSVX
300: *
301: 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