 LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ 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.EQ.0 and ITER.GE.0, see description below), then A is unchanged, if double precision factorization has been used (INFO.EQ.0 and ITER.LT.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.EQ.0 and ITER.GE.0) or the double precision factorization (if INFO.EQ.0 and ITER.LT.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.```
Date
June 2016

Definition at line 203 of file zcgesv.f.

203 *
204 * -- LAPACK driver routine (version 3.8.0) --
205 * -- LAPACK is a software package provided by Univ. of Tennessee, --
206 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
207 * June 2016
208 *
209 * .. Scalar Arguments ..
210  INTEGER info, iter, lda, ldb, ldx, n, nrhs
211 * ..
212 * .. Array Arguments ..
213  INTEGER ipiv( * )
214  DOUBLE PRECISION rwork( * )
215  COMPLEX swork( * )
216  COMPLEX*16 a( lda, * ), b( ldb, * ), work( n, * ),
217  \$ x( ldx, * )
218 * ..
219 *
220 * =====================================================================
221 *
222 * .. Parameters ..
223  LOGICAL doitref
224  parameter( doitref = .true. )
225 *
226  INTEGER itermax
227  parameter( itermax = 30 )
228 *
229  DOUBLE PRECISION bwdmax
230  parameter( bwdmax = 1.0e+00 )
231 *
232  COMPLEX*16 negone, one
233  parameter( negone = ( -1.0d+00, 0.0d+00 ),
234  \$ one = ( 1.0d+00, 0.0d+00 ) )
235 *
236 * .. Local Scalars ..
237  INTEGER i, iiter, ptsa, ptsx
238  DOUBLE PRECISION anrm, cte, eps, rnrm, xnrm
239  COMPLEX*16 zdum
240 *
241 * .. External Subroutines ..
242  EXTERNAL cgetrs, cgetrf, clag2z, xerbla, zaxpy, zgemm,
244 * ..
245 * .. External Functions ..
246  INTEGER izamax
247  DOUBLE PRECISION dlamch, zlange
248  EXTERNAL izamax, dlamch, zlange
249 * ..
250 * .. Intrinsic Functions ..
251  INTRINSIC abs, dble, max, sqrt
252 * ..
253 * .. Statement Functions ..
254  DOUBLE PRECISION cabs1
255 * ..
256 * .. Statement Function definitions ..
257  cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
258 * ..
259 * .. Executable Statements ..
260 *
261  info = 0
262  iter = 0
263 *
264 * Test the input parameters.
265 *
266  IF( n.LT.0 ) THEN
267  info = -1
268  ELSE IF( nrhs.LT.0 ) THEN
269  info = -2
270  ELSE IF( lda.LT.max( 1, n ) ) THEN
271  info = -4
272  ELSE IF( ldb.LT.max( 1, n ) ) THEN
273  info = -7
274  ELSE IF( ldx.LT.max( 1, n ) ) THEN
275  info = -9
276  END IF
277  IF( info.NE.0 ) THEN
278  CALL xerbla( 'ZCGESV', -info )
279  RETURN
280  END IF
281 *
282 * Quick return if (N.EQ.0).
283 *
284  IF( n.EQ.0 )
285  \$ RETURN
286 *
287 * Skip single precision iterative refinement if a priori slower
288 * than double precision factorization.
289 *
290  IF( .NOT.doitref ) THEN
291  iter = -1
292  GO TO 40
293  END IF
294 *
295 * Compute some constants.
296 *
297  anrm = zlange( 'I', n, n, a, lda, rwork )
298  eps = dlamch( 'Epsilon' )
299  cte = anrm*eps*sqrt( dble( n ) )*bwdmax
300 *
301 * Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
302 *
303  ptsa = 1
304  ptsx = ptsa + n*n
305 *
306 * Convert B from double precision to single precision and store the
307 * result in SX.
308 *
309  CALL zlag2c( n, nrhs, b, ldb, swork( ptsx ), n, info )
310 *
311  IF( info.NE.0 ) THEN
312  iter = -2
313  GO TO 40
314  END IF
315 *
316 * Convert A from double precision to single precision and store the
317 * result in SA.
318 *
319  CALL zlag2c( n, n, a, lda, swork( ptsa ), n, info )
320 *
321  IF( info.NE.0 ) THEN
322  iter = -2
323  GO TO 40
324  END IF
325 *
326 * Compute the LU factorization of SA.
327 *
328  CALL cgetrf( n, n, swork( ptsa ), n, ipiv, info )
329 *
330  IF( info.NE.0 ) THEN
331  iter = -3
332  GO TO 40
333  END IF
334 *
335 * Solve the system SA*SX = SB.
336 *
337  CALL cgetrs( 'No transpose', n, nrhs, swork( ptsa ), n, ipiv,
338  \$ swork( ptsx ), n, info )
339 *
340 * Convert SX back to double precision
341 *
342  CALL clag2z( n, nrhs, swork( ptsx ), n, x, ldx, info )
343 *
344 * Compute R = B - AX (R is WORK).
345 *
346  CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
347 *
348  CALL zgemm( 'No Transpose', 'No Transpose', n, nrhs, n, negone, a,
349  \$ lda, x, ldx, one, work, n )
350 *
351 * Check whether the NRHS normwise backward errors satisfy the
352 * stopping criterion. If yes, set ITER=0 and return.
353 *
354  DO i = 1, nrhs
355  xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
356  rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
357  IF( rnrm.GT.xnrm*cte )
358  \$ GO TO 10
359  END DO
360 *
361 * If we are here, the NRHS normwise backward errors satisfy the
362 * stopping criterion. We are good to exit.
363 *
364  iter = 0
365  RETURN
366 *
367  10 CONTINUE
368 *
369  DO 30 iiter = 1, itermax
370 *
371 * Convert R (in WORK) from double precision to single precision
372 * and store the result in SX.
373 *
374  CALL zlag2c( n, nrhs, work, n, swork( ptsx ), n, info )
375 *
376  IF( info.NE.0 ) THEN
377  iter = -2
378  GO TO 40
379  END IF
380 *
381 * Solve the system SA*SX = SR.
382 *
383  CALL cgetrs( 'No transpose', n, nrhs, swork( ptsa ), n, ipiv,
384  \$ swork( ptsx ), n, info )
385 *
386 * Convert SX back to double precision and update the current
387 * iterate.
388 *
389  CALL clag2z( n, nrhs, swork( ptsx ), n, work, n, info )
390 *
391  DO i = 1, nrhs
392  CALL zaxpy( n, one, work( 1, i ), 1, x( 1, i ), 1 )
393  END DO
394 *
395 * Compute R = B - AX (R is WORK).
396 *
397  CALL zlacpy( 'All', n, nrhs, b, ldb, work, n )
398 *
399  CALL zgemm( 'No Transpose', 'No Transpose', n, nrhs, n, negone,
400  \$ a, lda, x, ldx, one, work, n )
401 *
402 * Check whether the NRHS normwise backward errors satisfy the
403 * stopping criterion. If yes, set ITER=IITER>0 and return.
404 *
405  DO i = 1, nrhs
406  xnrm = cabs1( x( izamax( n, x( 1, i ), 1 ), i ) )
407  rnrm = cabs1( work( izamax( n, work( 1, i ), 1 ), i ) )
408  IF( rnrm.GT.xnrm*cte )
409  \$ GO TO 20
410  END DO
411 *
412 * If we are here, the NRHS normwise backward errors satisfy the
413 * stopping criterion, we are good to exit.
414 *
415  iter = iiter
416 *
417  RETURN
418 *
419  20 CONTINUE
420 *
421  30 CONTINUE
422 *
423 * If we are at this place of the code, this is because we have
424 * performed ITER=ITERMAX iterations and never satisified the stopping
425 * criterion, set up the ITER flag accordingly and follow up on double
426 * precision routine.
427 *
428  iter = -itermax - 1
429 *
430  40 CONTINUE
431 *
432 * Single-precision iterative refinement failed to converge to a
433 * satisfactory solution, so we resort to double precision.
434 *
435  CALL zgetrf( n, n, a, lda, ipiv, info )
436 *
437  IF( info.NE.0 )
438  \$ RETURN
439 *
440  CALL zlacpy( 'All', n, nrhs, b, ldb, x, ldx )
441  CALL zgetrs( 'No transpose', n, nrhs, a, lda, ipiv, x, ldx,
442  \$ info )
443 *
444  RETURN
445 *
446 * End of ZCGESV.
447 *
subroutine zgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZGETRS
Definition: zgetrs.f:123
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
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:109
subroutine cgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
CGETRS
Definition: cgetrs.f:123
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
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:105
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:73
subroutine zgetrf(M, N, A, LDA, IPIV, INFO)
ZGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
Definition: zgetrf.f:102
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:117
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cgetrf(M, N, A, LDA, IPIV, INFO)
CGETRF
Definition: cgetrf.f:110
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:90
Here is the call graph for this function:
Here is the caller graph for this function: