LAPACK  3.10.1 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 = 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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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:187
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
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
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 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 cgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
CGETRS
Definition: cgetrs.f:121
subroutine cgetrf(M, N, A, LDA, IPIV, INFO)
CGETRF
Definition: cgetrf.f:108
subroutine zgetrf(M, N, A, LDA, IPIV, INFO)
ZGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
Definition: zgetrf.f:102
Here is the call graph for this function:
Here is the caller graph for this function: