LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ zcgesv()

 subroutine zcgesv ( integer n, integer nrhs, complex*16, dimension( lda, * ) a, integer lda, integer, dimension( * ) ipiv, complex*16, dimension( ldb, * ) b, integer ldb, complex*16, dimension( ldx, * ) x, integer ldx, complex*16, dimension( n, * ) work, complex, dimension( * ) swork, double precision, dimension( * ) rwork, integer iter, integer info )

ZCGESV computes the solution to system of linear equations A * X = B for GE matrices (mixed precision with iterative refinement)

Purpose:
``` ZCGESV computes the solution to a complex system of linear equations
A * X = B,
where A is an N-by-N matrix and X and B are N-by-NRHS matrices.

ZCGESV first attempts to factorize the matrix in COMPLEX and use this
factorization within an iterative refinement procedure to produce a
solution with COMPLEX*16 normwise backward error quality (see below).
If the approach fails the method switches to a COMPLEX*16
factorization and solve.

The iterative refinement is not going to be a winning strategy if
the ratio COMPLEX performance over COMPLEX*16 performance is too
small. A reasonable strategy should take the number of right-hand
sides and the size of the matrix into account. This might be done
with a call to ILAENV in the future. Up to now, we always try
iterative refinement.

The iterative refinement process is stopped if
ITER > ITERMAX
or for all the RHS we have:
RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
where
o ITER is the number of the current iteration in the iterative
refinement process
o RNRM is the infinity-norm of the residual
o XNRM is the infinity-norm of the solution
o ANRM is the infinity-operator-norm of the matrix A
o EPS is the machine epsilon returned by DLAMCH('Epsilon')
The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
respectively.```
Parameters
 [in] N ``` N is INTEGER The number of linear equations, i.e., the order of the matrix A. N >= 0.``` [in] NRHS ``` NRHS is INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0.``` [in,out] A ``` A is COMPLEX*16 array, dimension (LDA,N) On entry, the N-by-N coefficient matrix A. On exit, if iterative refinement has been successfully used (INFO = 0 and ITER >= 0, see description below), then A is unchanged, if double precision factorization has been used (INFO = 0 and ITER < 0, see description below), then the array A contains the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [out] IPIV ``` IPIV is INTEGER array, dimension (N) The pivot indices that define the permutation matrix P; row i of the matrix was interchanged with row IPIV(i). Corresponds either to the single precision factorization (if INFO = 0 and ITER >= 0) or the double precision factorization (if INFO = 0 and ITER < 0).``` [in] B ``` B is COMPLEX*16 array, dimension (LDB,NRHS) The N-by-NRHS right hand side matrix B.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= max(1,N).``` [out] X ``` X is COMPLEX*16 array, dimension (LDX,NRHS) If INFO = 0, the N-by-NRHS solution matrix X.``` [in] LDX ``` LDX is INTEGER The leading dimension of the array X. LDX >= max(1,N).``` [out] WORK ``` WORK is COMPLEX*16 array, dimension (N,NRHS) This array is used to hold the residual vectors.``` [out] SWORK ``` SWORK is COMPLEX array, dimension (N*(N+NRHS)) This array is used to use the single precision matrix and the right-hand sides or solutions in single precision.``` [out] RWORK ` RWORK is DOUBLE PRECISION array, dimension (N)` [out] ITER ``` ITER is INTEGER < 0: iterative refinement has failed, COMPLEX*16 factorization has been performed -1 : the routine fell back to full precision for implementation- or machine-specific reasons -2 : narrowing the precision induced an overflow, the routine fell back to full precision -3 : failure of CGETRF -31: stop the iterative refinement after the 30th iterations > 0: iterative refinement has been successfully used. Returns the number of iterations``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, U(i,i) computed in COMPLEX*16 is exactly zero. The factorization has been completed, but the factor U is exactly singular, so the solution could not be computed.```

Definition at line 199 of file zcgesv.f.

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*
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 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
integer function izamax(n, zx, incx)
IZAMAX
Definition izamax.f:71
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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition zlange.f:115
Here is the call graph for this function:
Here is the caller graph for this function: