LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
chesvx.f
Go to the documentation of this file.
1 *> \brief <b> CHESVX 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
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chesvx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chesvx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chesvx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CHESVX( 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 * REAL RCOND
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IPIV( * )
32 * REAL BERR( * ), FERR( * ), RWORK( * )
33 * COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
34 * \$ WORK( * ), X( LDX, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> CHESVX 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 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 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 CHETRF.
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 CHETRF.
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 CHETRF.
171 *> \endverbatim
172 *>
173 *> \param[in] B
174 *> \verbatim
175 *> B is COMPLEX 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 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 REAL
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 REAL 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 REAL 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 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 CHETRF.
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 REAL 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 *> \date April 2012
280 *
281 *> \ingroup complexHEsolve
282 *
283 * =====================================================================
284  SUBROUTINE chesvx( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B,
285  \$ ldb, x, ldx, rcond, ferr, berr, work, lwork,
286  \$ rwork, info )
287 *
288 * -- LAPACK driver routine (version 3.4.1) --
289 * -- LAPACK is a software package provided by Univ. of Tennessee, --
290 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
291 * April 2012
292 *
293 * .. Scalar Arguments ..
294  CHARACTER fact, uplo
295  INTEGER info, lda, ldaf, ldb, ldx, lwork, n, nrhs
296  REAL rcond
297 * ..
298 * .. Array Arguments ..
299  INTEGER ipiv( * )
300  REAL berr( * ), ferr( * ), rwork( * )
301  COMPLEX a( lda, * ), af( ldaf, * ), b( ldb, * ),
302  \$ work( * ), x( ldx, * )
303 * ..
304 *
305 * =====================================================================
306 *
307 * .. Parameters ..
308  REAL zero
309  parameter( zero = 0.0e+0 )
310 * ..
311 * .. Local Scalars ..
312  LOGICAL lquery, nofact
313  INTEGER lwkopt, nb
314  REAL anorm
315 * ..
316 * .. External Functions ..
317  LOGICAL lsame
318  INTEGER ilaenv
319  REAL clanhe, slamch
320  EXTERNAL ilaenv, lsame, clanhe, slamch
321 * ..
322 * .. External Subroutines ..
323  EXTERNAL checon, cherfs, chetrf, chetrs, clacpy, xerbla
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, 'CHETRF', 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( 'CHESVX', -info )
367  return
368  ELSE IF( lquery ) THEN
369  return
370  END IF
371 *
372  IF( nofact ) THEN
373 *
374 * Compute the factorization A = U*D*U**H or A = L*D*L**H.
375 *
376  CALL clacpy( uplo, n, n, a, lda, af, ldaf )
377  CALL chetrf( 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 = clanhe( 'I', uplo, n, a, lda, rwork )
390 *
391 * Compute the reciprocal of the condition number of A.
392 *
393  CALL checon( uplo, n, af, ldaf, ipiv, anorm, rcond, work, info )
394 *
395 * Compute the solution vectors X.
396 *
397  CALL clacpy( 'Full', n, nrhs, b, ldb, x, ldx )
398  CALL chetrs( 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 cherfs( 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.slamch( 'Epsilon' ) )
409  \$ info = n + 1
410 *
411  work( 1 ) = lwkopt
412 *
413  return
414 *
415 * End of CHESVX
416 *
417  END