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