Annotation of rpl/lapack/lapack/zhesvx.f, revision 1.14
1.9 bertrand 1: *> \brief <b> ZHESVX computes the solution to system of linear equations A * X = B for HE 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 ZHESVX + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhesvx.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhesvx.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhesvx.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHESVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
22: * LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
23: * RWORK, 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( * )
32: * DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
33: * COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
34: * $ WORK( * ), X( LDX, * )
35: * ..
36: *
37: *
38: *> \par Purpose:
39: * =============
40: *>
41: *> \verbatim
42: *>
43: *> ZHESVX 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 Hermitian 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**H, if UPLO = 'U', or
62: *> A = L * D * L**H, if UPLO = 'L',
63: *> where U (or L) is a product of permutation and unit upper (lower)
64: *> triangular matrices, and D is Hermitian 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 Hermitian 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**H or A = L*D*L**H as computed by ZHETRF.
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**H or A = L*D*L**H.
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 ZHETRF.
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 ZHETRF.
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 ZHETRF.
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: *
274: *> \author Univ. of Tennessee
275: *> \author Univ. of California Berkeley
276: *> \author Univ. of Colorado Denver
277: *> \author NAG Ltd.
278: *
1.11 bertrand 279: *> \date April 2012
1.9 bertrand 280: *
281: *> \ingroup complex16HEsolve
282: *
283: * =====================================================================
1.1 bertrand 284: SUBROUTINE ZHESVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
285: $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
286: $ RWORK, INFO )
287: *
1.11 bertrand 288: * -- LAPACK driver routine (version 3.4.1) --
1.1 bertrand 289: * -- LAPACK is a software package provided by Univ. of Tennessee, --
290: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.11 bertrand 291: * April 2012
1.1 bertrand 292: *
293: * .. Scalar Arguments ..
294: CHARACTER FACT, UPLO
295: INTEGER INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
296: DOUBLE PRECISION RCOND
297: * ..
298: * .. Array Arguments ..
299: INTEGER IPIV( * )
300: DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
301: COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
302: $ WORK( * ), X( LDX, * )
303: * ..
304: *
305: * =====================================================================
306: *
307: * .. Parameters ..
308: DOUBLE PRECISION ZERO
309: PARAMETER ( ZERO = 0.0D+0 )
310: * ..
311: * .. Local Scalars ..
312: LOGICAL LQUERY, NOFACT
313: INTEGER LWKOPT, NB
314: DOUBLE PRECISION ANORM
315: * ..
316: * .. External Functions ..
317: LOGICAL LSAME
318: INTEGER ILAENV
319: DOUBLE PRECISION DLAMCH, ZLANHE
320: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANHE
321: * ..
322: * .. External Subroutines ..
323: EXTERNAL XERBLA, ZHECON, ZHERFS, ZHETRF, ZHETRS, ZLACPY
324: * ..
325: * .. Intrinsic Functions ..
326: INTRINSIC MAX
327: * ..
328: * .. Executable Statements ..
329: *
330: * Test the input parameters.
331: *
332: INFO = 0
333: NOFACT = LSAME( FACT, 'N' )
334: LQUERY = ( LWORK.EQ.-1 )
335: IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
336: INFO = -1
337: ELSE IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) )
338: $ THEN
339: INFO = -2
340: ELSE IF( N.LT.0 ) THEN
341: INFO = -3
342: ELSE IF( NRHS.LT.0 ) THEN
343: INFO = -4
344: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
345: INFO = -6
346: ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
347: INFO = -8
348: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
349: INFO = -11
350: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
351: INFO = -13
352: ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
353: INFO = -18
354: END IF
355: *
356: IF( INFO.EQ.0 ) THEN
357: LWKOPT = MAX( 1, 2*N )
358: IF( NOFACT ) THEN
359: NB = ILAENV( 1, 'ZHETRF', UPLO, N, -1, -1, -1 )
360: LWKOPT = MAX( LWKOPT, N*NB )
361: END IF
362: WORK( 1 ) = LWKOPT
363: END IF
364: *
365: IF( INFO.NE.0 ) THEN
366: CALL XERBLA( 'ZHESVX', -INFO )
367: RETURN
368: ELSE IF( LQUERY ) THEN
369: RETURN
370: END IF
371: *
372: IF( NOFACT ) THEN
373: *
1.8 bertrand 374: * Compute the factorization A = U*D*U**H or A = L*D*L**H.
1.1 bertrand 375: *
376: CALL ZLACPY( UPLO, N, N, A, LDA, AF, LDAF )
377: CALL ZHETRF( UPLO, N, AF, LDAF, IPIV, WORK, LWORK, INFO )
378: *
379: * Return if INFO is non-zero.
380: *
381: IF( INFO.GT.0 )THEN
382: RCOND = ZERO
383: RETURN
384: END IF
385: END IF
386: *
387: * Compute the norm of the matrix A.
388: *
389: ANORM = ZLANHE( 'I', UPLO, N, A, LDA, RWORK )
390: *
391: * Compute the reciprocal of the condition number of A.
392: *
393: CALL ZHECON( UPLO, N, AF, LDAF, IPIV, ANORM, RCOND, WORK, INFO )
394: *
395: * Compute the solution vectors X.
396: *
397: CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
398: CALL ZHETRS( UPLO, N, NRHS, AF, LDAF, IPIV, X, LDX, INFO )
399: *
400: * Use iterative refinement to improve the computed solutions and
401: * compute error bounds and backward error estimates for them.
402: *
403: CALL ZHERFS( UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X,
404: $ LDX, FERR, BERR, WORK, RWORK, INFO )
405: *
406: * Set INFO = N+1 if the matrix is singular to working precision.
407: *
408: IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
409: $ INFO = N + 1
410: *
411: WORK( 1 ) = LWKOPT
412: *
413: RETURN
414: *
415: * End of ZHESVX
416: *
417: END
CVSweb interface <joel.bertrand@systella.fr>