001:       SUBROUTINE ZCGESV( N, NRHS, A, LDA, IPIV, B, LDB, X, LDX, WORK,
002:      +                   SWORK, RWORK, ITER, INFO )
003: *
004: *  -- LAPACK PROTOTYPE driver routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     January 2007
007: *
008: *     ..
009: *     .. Scalar Arguments ..
010:       INTEGER            INFO, ITER, LDA, LDB, LDX, N, NRHS
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            IPIV( * )
014:       DOUBLE PRECISION   RWORK( * )
015:       COMPLEX            SWORK( * )
016:       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( N, * ),
017:      +                   X( LDX, * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  ZCGESV computes the solution to a complex system of linear equations
024: *     A * X = B,
025: *  where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
026: *
027: *  ZCGESV first attempts to factorize the matrix in COMPLEX and use this
028: *  factorization within an iterative refinement procedure to produce a
029: *  solution with COMPLEX*16 normwise backward error quality (see below).
030: *  If the approach fails the method switches to a COMPLEX*16
031: *  factorization and solve.
032: *
033: *  The iterative refinement is not going to be a winning strategy if
034: *  the ratio COMPLEX performance over COMPLEX*16 performance is too
035: *  small. A reasonable strategy should take the number of right-hand
036: *  sides and the size of the matrix into account. This might be done
037: *  with a call to ILAENV in the future. Up to now, we always try
038: *  iterative refinement.
039: *
040: *  The iterative refinement process is stopped if
041: *      ITER > ITERMAX
042: *  or for all the RHS we have:
043: *      RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX
044: *  where
045: *      o ITER is the number of the current iteration in the iterative
046: *        refinement process
047: *      o RNRM is the infinity-norm of the residual
048: *      o XNRM is the infinity-norm of the solution
049: *      o ANRM is the infinity-operator-norm of the matrix A
050: *      o EPS is the machine epsilon returned by DLAMCH('Epsilon')
051: *  The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00
052: *  respectively.
053: *
054: *  Arguments
055: *  =========
056: *
057: *  N       (input) INTEGER
058: *          The number of linear equations, i.e., the order of the
059: *          matrix A.  N >= 0.
060: *
061: *  NRHS    (input) INTEGER
062: *          The number of right hand sides, i.e., the number of columns
063: *          of the matrix B.  NRHS >= 0.
064: *
065: *  A       (input or input/ouptut) COMPLEX*16 array,
066: *          dimension (LDA,N)
067: *          On entry, the N-by-N coefficient matrix A.
068: *          On exit, if iterative refinement has been successfully used
069: *          (INFO.EQ.0 and ITER.GE.0, see description below), then A is
070: *          unchanged, if double precision factorization has been used
071: *          (INFO.EQ.0 and ITER.LT.0, see description below), then the
072: *          array A contains the factors L and U from the factorization
073: *          A = P*L*U; the unit diagonal elements of L are not stored.
074: *
075: *  LDA     (input) INTEGER
076: *          The leading dimension of the array A.  LDA >= max(1,N).
077: *
078: *  IPIV    (output) INTEGER array, dimension (N)
079: *          The pivot indices that define the permutation matrix P;
080: *          row i of the matrix was interchanged with row IPIV(i).
081: *          Corresponds either to the single precision factorization
082: *          (if INFO.EQ.0 and ITER.GE.0) or the double precision
083: *          factorization (if INFO.EQ.0 and ITER.LT.0).
084: *
085: *  B       (input) COMPLEX*16 array, dimension (LDB,NRHS)
086: *          The N-by-NRHS right hand side matrix B.
087: *
088: *  LDB     (input) INTEGER
089: *          The leading dimension of the array B.  LDB >= max(1,N).
090: *
091: *  X       (output) COMPLEX*16 array, dimension (LDX,NRHS)
092: *          If INFO = 0, the N-by-NRHS solution matrix X.
093: *
094: *  LDX     (input) INTEGER
095: *          The leading dimension of the array X.  LDX >= max(1,N).
096: *
097: *  WORK    (workspace) COMPLEX*16 array, dimension (N*NRHS)
098: *          This array is used to hold the residual vectors.
099: *
100: *  SWORK   (workspace) COMPLEX array, dimension (N*(N+NRHS))
101: *          This array is used to use the single precision matrix and the
102: *          right-hand sides or solutions in single precision.
103: *
104: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
105: *
106: *  ITER    (output) INTEGER
107: *          < 0: iterative refinement has failed, COMPLEX*16
108: *               factorization has been performed
109: *               -1 : the routine fell back to full precision for
110: *                    implementation- or machine-specific reasons
111: *               -2 : narrowing the precision induced an overflow,
112: *                    the routine fell back to full precision
113: *               -3 : failure of CGETRF
114: *               -31: stop the iterative refinement after the 30th
115: *                    iterations
116: *          > 0: iterative refinement has been sucessfully used.
117: *               Returns the number of iterations
118: *
119: *  INFO    (output) INTEGER
120: *          = 0:  successful exit
121: *          < 0:  if INFO = -i, the i-th argument had an illegal value
122: *          > 0:  if INFO = i, U(i,i) computed in COMPLEX*16 is exactly
123: *                zero.  The factorization has been completed, but the
124: *                factor U is exactly singular, so the solution
125: *                could not be computed.
126: *
127: *  =========
128: *
129: *     .. Parameters ..
130:       LOGICAL            DOITREF
131:       PARAMETER          ( DOITREF = .TRUE. )
132: *
133:       INTEGER            ITERMAX
134:       PARAMETER          ( ITERMAX = 30 )
135: *
136:       DOUBLE PRECISION   BWDMAX
137:       PARAMETER          ( BWDMAX = 1.0E+00 )
138: *
139:       COMPLEX*16         NEGONE, ONE
140:       PARAMETER          ( NEGONE = ( -1.0D+00, 0.0D+00 ),
141:      +                   ONE = ( 1.0D+00, 0.0D+00 ) )
142: *
143: *     .. Local Scalars ..
144:       INTEGER            I, IITER, PTSA, PTSX
145:       DOUBLE PRECISION   ANRM, CTE, EPS, RNRM, XNRM
146:       COMPLEX*16         ZDUM
147: *
148: *     .. External Subroutines ..
149:       EXTERNAL           CGETRS, CGETRF, CLAG2Z, XERBLA, ZAXPY, ZGEMM,
150:      +                   ZLACPY, ZLAG2C
151: *     ..
152: *     .. External Functions ..
153:       INTEGER            IZAMAX
154:       DOUBLE PRECISION   DLAMCH, ZLANGE
155:       EXTERNAL           IZAMAX, DLAMCH, ZLANGE
156: *     ..
157: *     .. Intrinsic Functions ..
158:       INTRINSIC          ABS, DBLE, MAX, SQRT
159: *     ..
160: *     .. Statement Functions ..
161:       DOUBLE PRECISION   CABS1
162: *     ..
163: *     .. Statement Function definitions ..
164:       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
165: *     ..
166: *     .. Executable Statements ..
167: *
168:       INFO = 0
169:       ITER = 0
170: *
171: *     Test the input parameters.
172: *
173:       IF( N.LT.0 ) THEN
174:          INFO = -1
175:       ELSE IF( NRHS.LT.0 ) THEN
176:          INFO = -2
177:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
178:          INFO = -4
179:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
180:          INFO = -7
181:       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
182:          INFO = -9
183:       END IF
184:       IF( INFO.NE.0 ) THEN
185:          CALL XERBLA( 'ZCGESV', -INFO )
186:          RETURN
187:       END IF
188: *
189: *     Quick return if (N.EQ.0).
190: *
191:       IF( N.EQ.0 )
192:      +   RETURN
193: *
194: *     Skip single precision iterative refinement if a priori slower
195: *     than double precision factorization.
196: *
197:       IF( .NOT.DOITREF ) THEN
198:          ITER = -1
199:          GO TO 40
200:       END IF
201: *
202: *     Compute some constants.
203: *
204:       ANRM = ZLANGE( 'I', N, N, A, LDA, RWORK )
205:       EPS = DLAMCH( 'Epsilon' )
206:       CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX
207: *
208: *     Set the indices PTSA, PTSX for referencing SA and SX in SWORK.
209: *
210:       PTSA = 1
211:       PTSX = PTSA + N*N
212: *
213: *     Convert B from double precision to single precision and store the
214: *     result in SX.
215: *
216:       CALL ZLAG2C( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO )
217: *
218:       IF( INFO.NE.0 ) THEN
219:          ITER = -2
220:          GO TO 40
221:       END IF
222: *
223: *     Convert A from double precision to single precision and store the
224: *     result in SA.
225: *
226:       CALL ZLAG2C( N, N, A, LDA, SWORK( PTSA ), N, INFO )
227: *
228:       IF( INFO.NE.0 ) THEN
229:          ITER = -2
230:          GO TO 40
231:       END IF
232: *
233: *     Compute the LU factorization of SA.
234: *
235:       CALL CGETRF( N, N, SWORK( PTSA ), N, IPIV, INFO )
236: *
237:       IF( INFO.NE.0 ) THEN
238:          ITER = -3
239:          GO TO 40
240:       END IF
241: *
242: *     Solve the system SA*SX = SB.
243: *
244:       CALL CGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV,
245:      +             SWORK( PTSX ), N, INFO )
246: *
247: *     Convert SX back to double precision
248: *
249:       CALL CLAG2Z( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO )
250: *
251: *     Compute R = B - AX (R is WORK).
252: *
253:       CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N )
254: *
255:       CALL ZGEMM( 'No Transpose', 'No Transpose', N, NRHS, N, NEGONE, A,
256:      +            LDA, X, LDX, ONE, WORK, N )
257: *
258: *     Check whether the NRHS normwise backward errors satisfy the
259: *     stopping criterion. If yes, set ITER=0 and return.
260: *
261:       DO I = 1, NRHS
262:          XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) )
263:          RNRM = CABS1( WORK( IZAMAX( N, WORK( 1, I ), 1 ), I ) )
264:          IF( RNRM.GT.XNRM*CTE )
265:      +      GO TO 10
266:       END DO
267: *
268: *     If we are here, the NRHS normwise backward errors satisfy the
269: *     stopping criterion. We are good to exit.
270: *
271:       ITER = 0
272:       RETURN
273: *
274:    10 CONTINUE
275: *
276:       DO 30 IITER = 1, ITERMAX
277: *
278: *        Convert R (in WORK) from double precision to single precision
279: *        and store the result in SX.
280: *
281:          CALL ZLAG2C( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO )
282: *
283:          IF( INFO.NE.0 ) THEN
284:             ITER = -2
285:             GO TO 40
286:          END IF
287: *
288: *        Solve the system SA*SX = SR.
289: *
290:          CALL CGETRS( 'No transpose', N, NRHS, SWORK( PTSA ), N, IPIV,
291:      +                SWORK( PTSX ), N, INFO )
292: *
293: *        Convert SX back to double precision and update the current
294: *        iterate.
295: *
296:          CALL CLAG2Z( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO )
297: *
298:          DO I = 1, NRHS
299:             CALL ZAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 )
300:          END DO
301: *
302: *        Compute R = B - AX (R is WORK).
303: *
304:          CALL ZLACPY( 'All', N, NRHS, B, LDB, WORK, N )
305: *
306:          CALL ZGEMM( 'No Transpose', 'No Transpose', N, NRHS, N, NEGONE,
307:      +               A, LDA, X, LDX, ONE, WORK, N )
308: *
309: *        Check whether the NRHS normwise backward errors satisfy the
310: *        stopping criterion. If yes, set ITER=IITER>0 and return.
311: *
312:          DO I = 1, NRHS
313:             XNRM = CABS1( X( IZAMAX( N, X( 1, I ), 1 ), I ) )
314:             RNRM = CABS1( WORK( IZAMAX( N, WORK( 1, I ), 1 ), I ) )
315:             IF( RNRM.GT.XNRM*CTE )
316:      +         GO TO 20
317:          END DO
318: *
319: *        If we are here, the NRHS normwise backward errors satisfy the
320: *        stopping criterion, we are good to exit.
321: *
322:          ITER = IITER
323: *
324:          RETURN
325: *
326:    20    CONTINUE
327: *
328:    30 CONTINUE
329: *
330: *     If we are at this place of the code, this is because we have
331: *     performed ITER=ITERMAX iterations and never satisified the stopping
332: *     criterion, set up the ITER flag accordingly and follow up on double
333: *     precision routine.
334: *
335:       ITER = -ITERMAX - 1
336: *
337:    40 CONTINUE
338: *
339: *     Single-precision iterative refinement failed to converge to a
340: *     satisfactory solution, so we resort to double precision.
341: *
342:       CALL ZGETRF( N, N, A, LDA, IPIV, INFO )
343: *
344:       IF( INFO.NE.0 )
345:      +   RETURN
346: *
347:       CALL ZLACPY( 'All', N, NRHS, B, LDB, X, LDX )
348:       CALL ZGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, X, LDX,
349:      +             INFO )
350: *
351:       RETURN
352: *
353: *     End of ZCGESV.
354: *
355:       END
356: