LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
zhesvx.f
Go to the documentation of this file.
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
136 *> AF is COMPLEX*16 array, dimension (LDAF,N)
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
156 *> IPIV is INTEGER array, dimension (N)
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 *
279 *> \ingroup complex16HEsolve
280 *
281 * =====================================================================
282  SUBROUTINE zhesvx( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
283  $ LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK,
284  $ RWORK, INFO )
285 *
286 * -- LAPACK driver routine --
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, ZLANHE
317  EXTERNAL lsame, ilaenv, dlamch, zlanhe
318 * ..
319 * .. External Subroutines ..
320  EXTERNAL xerbla, zhecon, zherfs, zhetrf, zhetrs, zlacpy
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, 'ZHETRF', 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( 'ZHESVX', -info )
364  RETURN
365  ELSE IF( lquery ) THEN
366  RETURN
367  END IF
368 *
369  IF( nofact ) THEN
370 *
371 * Compute the factorization A = U*D*U**H or A = L*D*L**H.
372 *
373  CALL zlacpy( uplo, n, n, a, lda, af, ldaf )
374  CALL zhetrf( 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 = zlanhe( 'I', uplo, n, a, lda, rwork )
387 *
388 * Compute the reciprocal of the condition number of A.
389 *
390  CALL zhecon( 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 zhetrs( 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 zherfs( 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 ZHESVX
413 *
414  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zherfs(UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
ZHERFS
Definition: zherfs.f:192
subroutine zhecon(UPLO, N, A, LDA, IPIV, ANORM, RCOND, WORK, INFO)
ZHECON
Definition: zhecon.f:125
subroutine zhetrs(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZHETRS
Definition: zhetrs.f:120
subroutine zhetrf(UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO)
ZHETRF
Definition: zhetrf.f:177
subroutine zhesvx(FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, RWORK, INFO)
ZHESVX computes the solution to system of linear equations A * X = B for HE matrices
Definition: zhesvx.f:285
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103