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