LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
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*> \ingroup hesvx
280*
281* =====================================================================
282 SUBROUTINE chesvx( 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 REAL RCOND
294* ..
295* .. Array Arguments ..
296 INTEGER IPIV( * )
297 REAL BERR( * ), FERR( * ), RWORK( * )
298 COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
299 \$ work( * ), x( ldx, * )
300* ..
301*
302* =====================================================================
303*
304* .. Parameters ..
305 REAL ZERO
306 PARAMETER ( ZERO = 0.0e+0 )
307* ..
308* .. Local Scalars ..
309 LOGICAL LQUERY, NOFACT
310 INTEGER LWKOPT, NB
311 REAL ANORM
312* ..
313* .. External Functions ..
314 LOGICAL LSAME
315 INTEGER ILAENV
316 REAL CLANHE, SLAMCH, SROUNDUP_LWORK
317 EXTERNAL ilaenv, lsame, clanhe, slamch, sroundup_lwork
318* ..
319* .. External Subroutines ..
320 EXTERNAL checon, cherfs, chetrf, chetrs, clacpy, xerbla
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, 'CHETRF', uplo, n, -1, -1, -1 )
357 lwkopt = max( lwkopt, n*nb )
358 END IF
359 work( 1 ) = sroundup_lwork(lwkopt)
360 END IF
361*
362 IF( info.NE.0 ) THEN
363 CALL xerbla( 'CHESVX', -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 clacpy( uplo, n, n, a, lda, af, ldaf )
374 CALL chetrf( 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 = clanhe( 'I', uplo, n, a, lda, rwork )
387*
388* Compute the reciprocal of the condition number of A.
389*
390 CALL checon( uplo, n, af, ldaf, ipiv, anorm, rcond, work, info )
391*
392* Compute the solution vectors X.
393*
394 CALL clacpy( 'Full', n, nrhs, b, ldb, x, ldx )
395 CALL chetrs( 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 cherfs( 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.slamch( 'Epsilon' ) )
406 \$ info = n + 1
407*
408 work( 1 ) = sroundup_lwork(lwkopt)
409*
410 RETURN
411*
412* End of CHESVX
413*
414 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine checon(uplo, n, a, lda, ipiv, anorm, rcond, work, info)
CHECON
Definition checon.f:125
subroutine cherfs(uplo, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x, ldx, ferr, berr, work, rwork, info)
CHERFS
Definition cherfs.f:192
subroutine chesvx(fact, uplo, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x, ldx, rcond, ferr, berr, work, lwork, rwork, info)
CHESVX computes the solution to system of linear equations A * X = B for HE matrices
Definition chesvx.f:285
subroutine chetrf(uplo, n, a, lda, ipiv, work, lwork, info)
CHETRF
Definition chetrf.f:177
subroutine chetrs(uplo, n, nrhs, a, lda, ipiv, b, ldb, info)
CHETRS
Definition chetrs.f:120
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103