1ZCGESV ‐ the solution to a real system of linear equations A * X
2= B, SUBROUTINE ZCGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX,
3WORK,
4 + SWORK, ITER, INFO)
5 INTEGER INFO,ITER,LDA,LDB,LDX,N,NRHS
6 INTEGER IPIV(*)
7 COMPLEX SWORK(*)
8 COMPLEX*16 A(LDA,*),B(LDB,*),WORK(N,*),X(LDX,*) ZCGESV com‐
9putes the solution to a real system of linear equations
10 A * X = B, where A is an N‐by‐N matrix and X and B are N‐by‐
11NRHS matrices.
12ZCGESV first attempts to factorize the matrix in SINGLE COMPLEX
14PRECISION and use this factorization within an iterative refine‐
15ment procedure to produce a solution with DOUBLE COMPLEX PRECI‐
16SION normwise backward error quality (see below). If the approach
17fails the method switches to a DOUBLE COMPLEX PRECISION factor‐
18ization and solve.
19The iterative refinement is not going to be a winning strategy if
21the ratio SINGLE PRECISION performance over DOUBLE PRECISION per‐
22formance is too small. A reasonable strategy should take the num‐
23ber of right‐hand sides and the size of the matrix into account.
24This might be done with a call to ILAENV in the future. Up to
25now, we always try iterative refinement.
26The iterative refinement process is stopped if
28 ITER > ITERMAX
29or for all the RHS we have:
30 RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
31where
32 o ITER is the number of the current iteration in the itera‐
33tive
34 refinement process
35 o RNRM is the infinity‐norm of the residual
36 o XNRM is the infinity‐norm of the solution
37 o ANRM is the infinity‐operator‐norm of the matrix A
38 o EPS is the machine epsilon returned by DLAMCH('Epsilon')
39The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 respec‐
40tively.
41N (input) INTEGER The number of linear equations, i.e., the
43order of the matrix A. N >= 0. NRHS (input) INTEGER The num‐
44ber of right hand sides, i.e., the number of columns of the ma‐
45trix B. NRHS >= 0. A (input or input/ouptut) COMPLEX*16
46array, dimension (LDA,N) On entry, the N‐by‐N coefficient matrix
47A. On exit, if iterative refinement has been successfully used
48(INFO.EQ.0 and ITER.GE.0, see description below), then A is un‐
49changed, if double precision factorization has been used (IN‐
50FO.EQ.0 and ITER.LT.0, see description below), then the array A
51contains the factors L and U from the factorization A = P*L*U;
52the unit diagonal elements of L are not stored. LDA (input)
53INTEGER The leading dimension of the array A. LDA >= max(1,N).
54IPIV (output) INTEGER array, dimension (N) The pivot indices
55that define the permutation matrix P; row i of the matrix was in‐
56terchanged with row IPIV(i). Corresponds either to the single
57precision factorization (if INFO.EQ.0 and ITER.GE.0) or the dou‐
58ble precision factorization (if INFO.EQ.0 and ITER.LT.0). B
59(input) COMPLEX*16 array, dimension (LDB,NRHS) The N‐by‐NRHS ma‐
60trix of right hand side matrix B. LDB (input) INTEGER The
61leading dimension of the array B. LDB >= max(1,N). X
62(output) COMPLEX*16 array, dimension (LDX,NRHS) If INFO = 0, the
63N‐by‐NRHS solution matrix X. LDX (input) INTEGER The leading
64dimension of the array X. LDX >= max(1,N). WORK (workspace)
65COMPLEX*16 array, dimension (N*NRHS) This array is used to hold
66the residual vectors. SWORK (workspace) COMPLEX array, dimen‐
67sion (N*(N+NRHS)) This array is used to use the single precision
68matrix and the right‐hand sides or solutions in single precision.
69ITER (output) INTEGER < 0: iterative refinement has failed,
70double precision factorization has been performed ‐1 : taking in‐
71to account machine parameters, N, NRHS, it is a priori not worth
72working in SINGLE PRECISION ‐2 : overflow of an entry when moving
73from double to SINGLE PRECISION ‐3 : failure of SGETRF
74‐31: stop the iterative refinement after the 30th iterations > 0:
75iterative refinement has been sucessfully used. Returns the num‐
76ber of iterations INFO (output) INTEGER = 0: successful exit
77< 0: if INFO = ‐i, the i‐th argument had an illegal value
78> 0: if INFO = i, U(i,i) computed in DOUBLE PRECISION is exactly
79zero. The factorization has been completed, but the factor U is
80exactly singular, so the solution could not be computed.
81=========
83