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