LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2016

Definition at line 203 of file zcgesv.f.

203 *
204 * -- LAPACK driver routine (version 3.6.1) --
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,
243  $ zlacpy, zlag2c
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 *
zgetrs
subroutine zgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
ZGETRS
Definition: zgetrs.f:123
zlacpy
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
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
zlag2c
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
cgetrs
subroutine cgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
CGETRS
Definition: cgetrs.f:123
zgetrf
subroutine zgetrf(M, N, A, LDA, IPIV, INFO)
ZGETRF VARIANT: Crout Level 3 BLAS version of the algorithm.
Definition: zgetrf.f:102
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zlange
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
izamax
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:53
cgetrf
subroutine cgetrf(M, N, A, LDA, IPIV, INFO)
CGETRF
Definition: cgetrf.f:110
clag2z
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
zaxpy
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: