LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
zcposv.f
Go to the documentation of this file.
1*> \brief <b> ZCPOSV computes the solution to system of linear equations A * X = B for PO 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/zcposv.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zcposv.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zcposv.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZCPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
22* SWORK, RWORK, ITER, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER UPLO
26* INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION RWORK( * )
30* COMPLEX SWORK( * )
31* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ),
32* \$ X( LDX, * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> ZCPOSV computes the solution to a complex system of linear equations
42*> A * X = B,
43*> where A is an N-by-N Hermitian positive definite matrix and X and B
44*> are N-by-NRHS matrices.
45*>
46*> ZCPOSV first attempts to factorize the matrix in COMPLEX and use this
47*> factorization within an iterative refinement procedure to produce a
48*> solution with COMPLEX*16 normwise backward error quality (see below).
49*> If the approach fails the method switches to a COMPLEX*16
50*> factorization and solve.
51*>
52*> The iterative refinement is not going to be a winning strategy if
53*> the ratio COMPLEX performance over COMPLEX*16 performance is too
54*> small. A reasonable strategy should take the number of right-hand
55*> sides and the size of the matrix into account. This might be done
56*> with a call to ILAENV in the future. Up to now, we always try
57*> iterative refinement.
58*>
59*> The iterative refinement process is stopped if
60*> ITER > ITERMAX
61*> or for all the RHS we have:
62*> RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
63*> where
64*> o ITER is the number of the current iteration in the iterative
65*> refinement process
66*> o RNRM is the infinity-norm of the residual
67*> o XNRM is the infinity-norm of the solution
68*> o ANRM is the infinity-operator-norm of the matrix A
69*> o EPS is the machine epsilon returned by DLAMCH('Epsilon')
70*> The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
71*> respectively.
72*> \endverbatim
73*
74* Arguments:
75* ==========
76*
77*> \param[in] UPLO
78*> \verbatim
79*> UPLO is CHARACTER*1
80*> = 'U': Upper triangle of A is stored;
81*> = 'L': Lower triangle of A is stored.
82*> \endverbatim
83*>
84*> \param[in] N
85*> \verbatim
86*> N is INTEGER
87*> The number of linear equations, i.e., the order of the
88*> matrix A. N >= 0.
89*> \endverbatim
90*>
91*> \param[in] NRHS
92*> \verbatim
93*> NRHS is INTEGER
94*> The number of right hand sides, i.e., the number of columns
95*> of the matrix B. NRHS >= 0.
96*> \endverbatim
97*>
98*> \param[in,out] A
99*> \verbatim
100*> A is COMPLEX*16 array,
101*> dimension (LDA,N)
102*> On entry, the Hermitian matrix A. If UPLO = 'U', the leading
103*> N-by-N upper triangular part of A contains the upper
104*> triangular part of the matrix A, and the strictly lower
105*> triangular part of A is not referenced. If UPLO = 'L', the
106*> leading N-by-N lower triangular part of A contains the lower
107*> triangular part of the matrix A, and the strictly upper
108*> triangular part of A is not referenced.
109*>
110*> Note that the imaginary parts of the diagonal
111*> elements need not be set and are assumed to be zero.
112*>
113*> On exit, if iterative refinement has been successfully used
114*> (INFO = 0 and ITER >= 0, see description below), then A is
115*> unchanged, if double precision factorization has been used
116*> (INFO = 0 and ITER < 0, see description below), then the
117*> array A contains the factor U or L from the Cholesky
118*> factorization A = U**H*U or A = L*L**H.
119*> \endverbatim
120*>
121*> \param[in] LDA
122*> \verbatim
123*> LDA is INTEGER
124*> The leading dimension of the array A. LDA >= max(1,N).
125*> \endverbatim
126*>
127*> \param[in] B
128*> \verbatim
129*> B is COMPLEX*16 array, dimension (LDB,NRHS)
130*> The N-by-NRHS right hand side matrix B.
131*> \endverbatim
132*>
133*> \param[in] LDB
134*> \verbatim
135*> LDB is INTEGER
136*> The leading dimension of the array B. LDB >= max(1,N).
137*> \endverbatim
138*>
139*> \param[out] X
140*> \verbatim
141*> X is COMPLEX*16 array, dimension (LDX,NRHS)
142*> If INFO = 0, the N-by-NRHS solution matrix X.
143*> \endverbatim
144*>
145*> \param[in] LDX
146*> \verbatim
147*> LDX is INTEGER
148*> The leading dimension of the array X. LDX >= max(1,N).
149*> \endverbatim
150*>
151*> \param[out] WORK
152*> \verbatim
153*> WORK is COMPLEX*16 array, dimension (N,NRHS)
154*> This array is used to hold the residual vectors.
155*> \endverbatim
156*>
157*> \param[out] SWORK
158*> \verbatim
159*> SWORK is COMPLEX array, dimension (N*(N+NRHS))
160*> This array is used to use the single precision matrix and the
161*> right-hand sides or solutions in single precision.
162*> \endverbatim
163*>
164*> \param[out] RWORK
165*> \verbatim
166*> RWORK is DOUBLE PRECISION array, dimension (N)
167*> \endverbatim
168*>
169*> \param[out] ITER
170*> \verbatim
171*> ITER is INTEGER
172*> < 0: iterative refinement has failed, COMPLEX*16
173*> factorization has been performed
174*> -1 : the routine fell back to full precision for
175*> implementation- or machine-specific reasons
176*> -2 : narrowing the precision induced an overflow,
177*> the routine fell back to full precision
178*> -3 : failure of CPOTRF
179*> -31: stop the iterative refinement after the 30th
180*> iterations
181*> > 0: iterative refinement has been successfully used.
182*> Returns the number of iterations
183*> \endverbatim
184*>
185*> \param[out] INFO
186*> \verbatim
187*> INFO is INTEGER
188*> = 0: successful exit
189*> < 0: if INFO = -i, the i-th argument had an illegal value
190*> > 0: if INFO = i, the leading principal minor of order i
191*> of (COMPLEX*16) A is not positive, so the factorization
192*> could not be completed, and the solution has not been
193*> computed.
194*> \endverbatim
195*
196* Authors:
197* ========
198*
199*> \author Univ. of Tennessee
200*> \author Univ. of California Berkeley
201*> \author Univ. of Colorado Denver
202*> \author NAG Ltd.
203*
204*> \ingroup posv_mixed
205*
206* =====================================================================
207 SUBROUTINE zcposv( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK,
208 \$ SWORK, RWORK, ITER, INFO )
209*
210* -- LAPACK driver routine --
211* -- LAPACK is a software package provided by Univ. of Tennessee, --
212* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213*
214* .. Scalar Arguments ..
215 CHARACTER UPLO
216 INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS
217* ..
218* .. Array Arguments ..
219 DOUBLE PRECISION RWORK( * )
220 COMPLEX SWORK( * )
221 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( N, * ),
222 \$ x( ldx, * )
223* ..
224*
225* =====================================================================
226*
227* .. Parameters ..
228 LOGICAL DOITREF
229 parameter( doitref = .true. )
230*
231 INTEGER ITERMAX
232 parameter( itermax = 30 )
233*
234 DOUBLE PRECISION BWDMAX
235 parameter( bwdmax = 1.0e+00 )
236*
237 COMPLEX*16 NEGONE, ONE
238 parameter( negone = ( -1.0d+00, 0.0d+00 ),
239 \$ one = ( 1.0d+00, 0.0d+00 ) )
240*
241* .. Local Scalars ..
242 INTEGER I, IITER, PTSA, PTSX
243 DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM
244 COMPLEX*16 ZDUM
245*
246* .. External Subroutines ..
247 EXTERNAL zaxpy, zhemm, zlacpy, zlat2c, zlag2c, clag2z,
249* ..
250* .. External Functions ..
251 INTEGER IZAMAX
252 DOUBLE PRECISION DLAMCH, ZLANHE
253 LOGICAL LSAME
254 EXTERNAL izamax, dlamch, zlanhe, lsame
255* ..
256* .. Intrinsic Functions ..
257 INTRINSIC abs, dble, max, sqrt
258* .. Statement Functions ..
259 DOUBLE PRECISION CABS1
260* ..
261* .. Statement Function definitions ..
262 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
263* ..
264* .. Executable Statements ..
265*
266 info = 0
267 iter = 0
268*
269* Test the input parameters.
270*
271 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lsame( uplo, 'L' ) ) THEN
272 info = -1
273 ELSE IF( n.LT.0 ) THEN
274 info = -2
275 ELSE IF( nrhs.LT.0 ) THEN
276 info = -3
277 ELSE IF( lda.LT.max( 1, n ) ) THEN
278 info = -5
279 ELSE IF( ldb.LT.max( 1, n ) ) THEN
280 info = -7
281 ELSE IF( ldx.LT.max( 1, n ) ) THEN
282 info = -9
283 END IF
284 IF( info.NE.0 ) THEN
285 CALL xerbla( 'ZCPOSV', -info )
286 RETURN
287 END IF
288*
289* Quick return if (N.EQ.0).
290*
291 IF( n.EQ.0 )
292 \$ RETURN
293*
294* Skip single precision iterative refinement if a priori slower
295* than double precision factorization.
296*
297 IF( .NOT.doitref ) THEN
298 iter = -1
299 GO TO 40
300 END IF
301*
302* Compute some constants.
303*
304 anrm = zlanhe( 'I', uplo, n, a, lda, rwork )
305 eps = dlamch( 'Epsilon' )
306 cte = anrm*eps*sqrt( dble( n ) )*bwdmax
307*
308* Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
309*
310 ptsa = 1
311 ptsx = ptsa + n*n
312*
313* Convert B from double precision to single precision and store the
314* result in SX.
315*
316 CALL zlag2c( n, nrhs, b, ldb, swork( ptsx ), n, info )
317*
318 IF( info.NE.0 ) THEN
319 iter = -2
320 GO TO 40
321 END IF
322*
323* Convert A from double precision to single precision and store the
324* result in SA.
325*
326 CALL zlat2c( uplo, n, a, lda, swork( ptsa ), n, info )
327*
328 IF( info.NE.0 ) THEN
329 iter = -2
330 GO TO 40
331 END IF
332*
333* Compute the Cholesky factorization of SA.
334*
335 CALL cpotrf( uplo, n, swork( ptsa ), n, info )
336*
337 IF( info.NE.0 ) THEN
338 iter = -3
339 GO TO 40
340 END IF
341*
342* Solve the system SA*SX = SB.
343*
344 CALL cpotrs( uplo, n, nrhs, swork( ptsa ), n, swork( ptsx ), n,
345 \$ info )
346*
347* Convert SX back to COMPLEX*16
348*
349 CALL clag2z( n, nrhs, swork( ptsx ), n, x, ldx, info )
350*
351* Compute R = B - AX (R is WORK).
352*
353 CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
354*
355 CALL zhemm( 'Left', uplo, n, nrhs, negone, a, lda, x, ldx, one,
356 \$ work, n )
357*
358* Check whether the NRHS normwise backward errors satisfy the
359* stopping criterion. If yes, set ITER=0 and return.
360*
361 DO i = 1, nrhs
362 xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
363 rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
364 IF( rnrm.GT.xnrm*cte )
365 \$ GO TO 10
366 END DO
367*
368* If we are here, the NRHS normwise backward errors satisfy the
369* stopping criterion. We are good to exit.
370*
371 iter = 0
372 RETURN
373*
374 10 CONTINUE
375*
376 DO 30 iiter = 1, itermax
377*
378* Convert R (in WORK) from double precision to single precision
379* and store the result in SX.
380*
381 CALL zlag2c( n, nrhs, work, n, swork( ptsx ), n, info )
382*
383 IF( info.NE.0 ) THEN
384 iter = -2
385 GO TO 40
386 END IF
387*
388* Solve the system SA*SX = SR.
389*
390 CALL cpotrs( uplo, n, nrhs, swork( ptsa ), n, swork( ptsx ), n,
391 \$ info )
392*
393* Convert SX back to double precision and update the current
394* iterate.
395*
396 CALL clag2z( n, nrhs, swork( ptsx ), n, work, n, info )
397*
398 DO i = 1, nrhs
399 CALL zaxpy( n, one, work( 1, i ), 1, x( 1, i ), 1 )
400 END DO
401*
402* Compute R = B - AX (R is WORK).
403*
404 CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
405*
406 CALL zhemm( 'L', uplo, n, nrhs, negone, a, lda, x, ldx, one,
407 \$ work, n )
408*
409* Check whether the NRHS normwise backward errors satisfy the
410* stopping criterion. If yes, set ITER=IITER>0 and return.
411*
412 DO i = 1, nrhs
413 xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
414 rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
415 IF( rnrm.GT.xnrm*cte )
416 \$ GO TO 20
417 END DO
418*
419* If we are here, the NRHS normwise backward errors satisfy the
420* stopping criterion, we are good to exit.
421*
422 iter = iiter
423*
424 RETURN
425*
426 20 CONTINUE
427*
428 30 CONTINUE
429*
430* If we are at this place of the code, this is because we have
431* performed ITER=ITERMAX iterations and never satisfied the
432* stopping criterion, set up the ITER flag accordingly and follow
433* up on double precision routine.
434*
435 iter = -itermax - 1
436*
437 40 CONTINUE
438*
439* Single-precision iterative refinement failed to converge to a
440* satisfactory solution, so we resort to double precision.
441*
442 CALL zpotrf( uplo, n, a, lda, info )
443*
444 IF( info.NE.0 )
445 \$ RETURN
446*
447 CALL zlacpy( 'All', n, nrhs, b, ldb, x, ldx )
448 CALL zpotrs( uplo, n, nrhs, a, lda, x, ldx, info )
449*
450 RETURN
451*
452* End of ZCPOSV
453*
454 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine clag2z(m, n, sa, ldsa, a, lda, info)
CLAG2Z converts a complex single precision matrix to a complex double precision matrix.
Definition clag2z.f:103
subroutine zlag2c(m, n, a, lda, sa, ldsa, info)
ZLAG2C converts a complex double precision matrix to a complex single precision matrix.
Definition zlag2c.f:107
subroutine zlat2c(uplo, n, a, lda, sa, ldsa, info)
ZLAT2C converts a double complex triangular matrix to a complex triangular matrix.
Definition zlat2c.f:111
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
subroutine zhemm(side, uplo, m, n, alpha, a, lda, b, ldb, beta, c, ldc)
ZHEMM
Definition zhemm.f:191
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
subroutine zcposv(uplo, n, nrhs, a, lda, b, ldb, x, ldx, work, swork, rwork, iter, info)
ZCPOSV computes the solution to system of linear equations A * X = B for PO matrices
Definition zcposv.f:209
subroutine zpotrf(uplo, n, a, lda, info)
ZPOTRF
Definition zpotrf.f:107
subroutine cpotrf(uplo, n, a, lda, info)
CPOTRF
Definition cpotrf.f:107
subroutine cpotrs(uplo, n, nrhs, a, lda, b, ldb, info)
CPOTRS
Definition cpotrs.f:110
subroutine zpotrs(uplo, n, nrhs, a, lda, b, ldb, info)
ZPOTRS
Definition zpotrs.f:110