LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
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)

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,
240 $ zlacpy, zlag2c, zgetrf, zgetrs
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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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:103
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:107
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
cgetrf
subroutine cgetrf(m, n, a, lda, ipiv, info)
CGETRF
Definition cgetrf.f:108
zgetrf
subroutine zgetrf(m, n, a, lda, ipiv, info)
ZGETRF
Definition zgetrf.f:108
cgetrs
subroutine cgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
CGETRS
Definition cgetrs.f:121
zgetrs
subroutine zgetrs(trans, n, nrhs, a, lda, ipiv, b, ldb, info)
ZGETRS
Definition zgetrs.f:121
izamax
integer function izamax(n, zx, incx)
IZAMAX
Definition izamax.f:71
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:103
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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:115
Here is the call graph for this function:
Here is the caller graph for this function: