LAPACK  3.10.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)

Download ZCGESV + dependencies [TGZ] [ZIP] [TXT]

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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: