001:       SUBROUTINE CPBSVX( FACT, UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB,
002:      $                   EQUED, S, B, LDB, X, LDX, RCOND, FERR, BERR,
003:      $                   WORK, RWORK, INFO )
004: *
005: *  -- LAPACK driver routine (version 3.2) --
006: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
007: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
008: *     November 2006
009: *
010: *     .. Scalar Arguments ..
011:       CHARACTER          EQUED, FACT, UPLO
012:       INTEGER            INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS
013:       REAL               RCOND
014: *     ..
015: *     .. Array Arguments ..
016:       REAL               BERR( * ), FERR( * ), RWORK( * ), S( * )
017:       COMPLEX            AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
018:      $                   WORK( * ), X( LDX, * )
019: *     ..
020: *
021: *  Purpose
022: *  =======
023: *
024: *  CPBSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
025: *  compute the solution to a complex system of linear equations
026: *     A * X = B,
027: *  where A is an N-by-N Hermitian positive definite band matrix and X
028: *  and B are N-by-NRHS matrices.
029: *
030: *  Error bounds on the solution and a condition estimate are also
031: *  provided.
032: *
033: *  Description
034: *  ===========
035: *
036: *  The following steps are performed:
037: *
038: *  1. If FACT = 'E', real scaling factors are computed to equilibrate
039: *     the system:
040: *        diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
041: *     Whether or not the system will be equilibrated depends on the
042: *     scaling of the matrix A, but if equilibration is used, A is
043: *     overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
044: *
045: *  2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
046: *     factor the matrix A (after equilibration if FACT = 'E') as
047: *        A = U**H * U,  if UPLO = 'U', or
048: *        A = L * L**H,  if UPLO = 'L',
049: *     where U is an upper triangular band matrix, and L is a lower
050: *     triangular band matrix.
051: *
052: *  3. If the leading i-by-i principal minor is not positive definite,
053: *     then the routine returns with INFO = i. Otherwise, the factored
054: *     form of A is used to estimate the condition number of the matrix
055: *     A.  If the reciprocal of the condition number is less than machine
056: *     precision, INFO = N+1 is returned as a warning, but the routine
057: *     still goes on to solve for X and compute error bounds as
058: *     described below.
059: *
060: *  4. The system of equations is solved for X using the factored form
061: *     of A.
062: *
063: *  5. Iterative refinement is applied to improve the computed solution
064: *     matrix and calculate error bounds and backward error estimates
065: *     for it.
066: *
067: *  6. If equilibration was used, the matrix X is premultiplied by
068: *     diag(S) so that it solves the original system before
069: *     equilibration.
070: *
071: *  Arguments
072: *  =========
073: *
074: *  FACT    (input) CHARACTER*1
075: *          Specifies whether or not the factored form of the matrix A is
076: *          supplied on entry, and if not, whether the matrix A should be
077: *          equilibrated before it is factored.
078: *          = 'F':  On entry, AFB contains the factored form of A.
079: *                  If EQUED = 'Y', the matrix A has been equilibrated
080: *                  with scaling factors given by S.  AB and AFB will not
081: *                  be modified.
082: *          = 'N':  The matrix A will be copied to AFB and factored.
083: *          = 'E':  The matrix A will be equilibrated if necessary, then
084: *                  copied to AFB and factored.
085: *
086: *  UPLO    (input) CHARACTER*1
087: *          = 'U':  Upper triangle of A is stored;
088: *          = 'L':  Lower triangle of A is stored.
089: *
090: *  N       (input) INTEGER
091: *          The number of linear equations, i.e., the order of the
092: *          matrix A.  N >= 0.
093: *
094: *  KD      (input) INTEGER
095: *          The number of superdiagonals of the matrix A if UPLO = 'U',
096: *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
097: *
098: *  NRHS    (input) INTEGER
099: *          The number of right-hand sides, i.e., the number of columns
100: *          of the matrices B and X.  NRHS >= 0.
101: *
102: *  AB      (input/output) COMPLEX array, dimension (LDAB,N)
103: *          On entry, the upper or lower triangle of the Hermitian band
104: *          matrix A, stored in the first KD+1 rows of the array, except
105: *          if FACT = 'F' and EQUED = 'Y', then A must contain the
106: *          equilibrated matrix diag(S)*A*diag(S).  The j-th column of A
107: *          is stored in the j-th column of the array AB as follows:
108: *          if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j;
109: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(N,j+KD).
110: *          See below for further details.
111: *
112: *          On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
113: *          diag(S)*A*diag(S).
114: *
115: *  LDAB    (input) INTEGER
116: *          The leading dimension of the array A.  LDAB >= KD+1.
117: *
118: *  AFB     (input or output) COMPLEX array, dimension (LDAFB,N)
119: *          If FACT = 'F', then AFB is an input argument and on entry
120: *          contains the triangular factor U or L from the Cholesky
121: *          factorization A = U**H*U or A = L*L**H of the band matrix
122: *          A, in the same storage format as A (see AB).  If EQUED = 'Y',
123: *          then AFB is the factored form of the equilibrated matrix A.
124: *
125: *          If FACT = 'N', then AFB is an output argument and on exit
126: *          returns the triangular factor U or L from the Cholesky
127: *          factorization A = U**H*U or A = L*L**H.
128: *
129: *          If FACT = 'E', then AFB is an output argument and on exit
130: *          returns the triangular factor U or L from the Cholesky
131: *          factorization A = U**H*U or A = L*L**H of the equilibrated
132: *          matrix A (see the description of A for the form of the
133: *          equilibrated matrix).
134: *
135: *  LDAFB   (input) INTEGER
136: *          The leading dimension of the array AFB.  LDAFB >= KD+1.
137: *
138: *  EQUED   (input or output) CHARACTER*1
139: *          Specifies the form of equilibration that was done.
140: *          = 'N':  No equilibration (always true if FACT = 'N').
141: *          = 'Y':  Equilibration was done, i.e., A has been replaced by
142: *                  diag(S) * A * diag(S).
143: *          EQUED is an input argument if FACT = 'F'; otherwise, it is an
144: *          output argument.
145: *
146: *  S       (input or output) REAL array, dimension (N)
147: *          The scale factors for A; not accessed if EQUED = 'N'.  S is
148: *          an input argument if FACT = 'F'; otherwise, S is an output
149: *          argument.  If FACT = 'F' and EQUED = 'Y', each element of S
150: *          must be positive.
151: *
152: *  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
153: *          On entry, the N-by-NRHS right hand side matrix B.
154: *          On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
155: *          B is overwritten by diag(S) * B.
156: *
157: *  LDB     (input) INTEGER
158: *          The leading dimension of the array B.  LDB >= max(1,N).
159: *
160: *  X       (output) COMPLEX array, dimension (LDX,NRHS)
161: *          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
162: *          the original system of equations.  Note that if EQUED = 'Y',
163: *          A and B are modified on exit, and the solution to the
164: *          equilibrated system is inv(diag(S))*X.
165: *
166: *  LDX     (input) INTEGER
167: *          The leading dimension of the array X.  LDX >= max(1,N).
168: *
169: *  RCOND   (output) REAL
170: *          The estimate of the reciprocal condition number of the matrix
171: *          A after equilibration (if done).  If RCOND is less than the
172: *          machine precision (in particular, if RCOND = 0), the matrix
173: *          is singular to working precision.  This condition is
174: *          indicated by a return code of INFO > 0.
175: *
176: *  FERR    (output) REAL array, dimension (NRHS)
177: *          The estimated forward error bound for each solution vector
178: *          X(j) (the j-th column of the solution matrix X).
179: *          If XTRUE is the true solution corresponding to X(j), FERR(j)
180: *          is an estimated upper bound for the magnitude of the largest
181: *          element in (X(j) - XTRUE) divided by the magnitude of the
182: *          largest element in X(j).  The estimate is as reliable as
183: *          the estimate for RCOND, and is almost always a slight
184: *          overestimate of the true error.
185: *
186: *  BERR    (output) REAL array, dimension (NRHS)
187: *          The componentwise relative backward error of each solution
188: *          vector X(j) (i.e., the smallest relative change in
189: *          any element of A or B that makes X(j) an exact solution).
190: *
191: *  WORK    (workspace) COMPLEX array, dimension (2*N)
192: *
193: *  RWORK   (workspace) REAL array, dimension (N)
194: *
195: *  INFO    (output) INTEGER
196: *          = 0: successful exit
197: *          < 0: if INFO = -i, the i-th argument had an illegal value
198: *          > 0: if INFO = i, and i is
199: *                <= N:  the leading minor of order i of A is
200: *                       not positive definite, so the factorization
201: *                       could not be completed, and the solution has not
202: *                       been computed. RCOND = 0 is returned.
203: *                = N+1: U is nonsingular, but RCOND is less than machine
204: *                       precision, meaning that the matrix is singular
205: *                       to working precision.  Nevertheless, the
206: *                       solution and error bounds are computed because
207: *                       there are a number of situations where the
208: *                       computed solution can be more accurate than the
209: *                       value of RCOND would suggest.
210: *
211: *  Further Details
212: *  ===============
213: *
214: *  The band storage scheme is illustrated by the following example, when
215: *  N = 6, KD = 2, and UPLO = 'U':
216: *
217: *  Two-dimensional storage of the Hermitian matrix A:
218: *
219: *     a11  a12  a13
220: *          a22  a23  a24
221: *               a33  a34  a35
222: *                    a44  a45  a46
223: *                         a55  a56
224: *     (aij=conjg(aji))         a66
225: *
226: *  Band storage of the upper triangle of A:
227: *
228: *      *    *   a13  a24  a35  a46
229: *      *   a12  a23  a34  a45  a56
230: *     a11  a22  a33  a44  a55  a66
231: *
232: *  Similarly, if UPLO = 'L' the format of A is as follows:
233: *
234: *     a11  a22  a33  a44  a55  a66
235: *     a21  a32  a43  a54  a65   *
236: *     a31  a42  a53  a64   *    *
237: *
238: *  Array elements marked * are not used by the routine.
239: *
240: *  =====================================================================
241: *
242: *     .. Parameters ..
243:       REAL               ZERO, ONE
244:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
245: *     ..
246: *     .. Local Scalars ..
247:       LOGICAL            EQUIL, NOFACT, RCEQU, UPPER
248:       INTEGER            I, INFEQU, J, J1, J2
249:       REAL               AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
250: *     ..
251: *     .. External Functions ..
252:       LOGICAL            LSAME
253:       REAL               CLANHB, SLAMCH
254:       EXTERNAL           LSAME, CLANHB, SLAMCH
255: *     ..
256: *     .. External Subroutines ..
257:       EXTERNAL           CCOPY, CLACPY, CLAQHB, CPBCON, CPBEQU, CPBRFS,
258:      $                   CPBTRF, CPBTRS, XERBLA
259: *     ..
260: *     .. Intrinsic Functions ..
261:       INTRINSIC          MAX, MIN
262: *     ..
263: *     .. Executable Statements ..
264: *
265:       INFO = 0
266:       NOFACT = LSAME( FACT, 'N' )
267:       EQUIL = LSAME( FACT, 'E' )
268:       UPPER = LSAME( UPLO, 'U' )
269:       IF( NOFACT .OR. EQUIL ) THEN
270:          EQUED = 'N'
271:          RCEQU = .FALSE.
272:       ELSE
273:          RCEQU = LSAME( EQUED, 'Y' )
274:          SMLNUM = SLAMCH( 'Safe minimum' )
275:          BIGNUM = ONE / SMLNUM
276:       END IF
277: *
278: *     Test the input parameters.
279: *
280:       IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
281:      $     THEN
282:          INFO = -1
283:       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
284:          INFO = -2
285:       ELSE IF( N.LT.0 ) THEN
286:          INFO = -3
287:       ELSE IF( KD.LT.0 ) THEN
288:          INFO = -4
289:       ELSE IF( NRHS.LT.0 ) THEN
290:          INFO = -5
291:       ELSE IF( LDAB.LT.KD+1 ) THEN
292:          INFO = -7
293:       ELSE IF( LDAFB.LT.KD+1 ) THEN
294:          INFO = -9
295:       ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
296:      $         ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
297:          INFO = -10
298:       ELSE
299:          IF( RCEQU ) THEN
300:             SMIN = BIGNUM
301:             SMAX = ZERO
302:             DO 10 J = 1, N
303:                SMIN = MIN( SMIN, S( J ) )
304:                SMAX = MAX( SMAX, S( J ) )
305:    10       CONTINUE
306:             IF( SMIN.LE.ZERO ) THEN
307:                INFO = -11
308:             ELSE IF( N.GT.0 ) THEN
309:                SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
310:             ELSE
311:                SCOND = ONE
312:             END IF
313:          END IF
314:          IF( INFO.EQ.0 ) THEN
315:             IF( LDB.LT.MAX( 1, N ) ) THEN
316:                INFO = -13
317:             ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
318:                INFO = -15
319:             END IF
320:          END IF
321:       END IF
322: *
323:       IF( INFO.NE.0 ) THEN
324:          CALL XERBLA( 'CPBSVX', -INFO )
325:          RETURN
326:       END IF
327: *
328:       IF( EQUIL ) THEN
329: *
330: *        Compute row and column scalings to equilibrate the matrix A.
331: *
332:          CALL CPBEQU( UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, INFEQU )
333:          IF( INFEQU.EQ.0 ) THEN
334: *
335: *           Equilibrate the matrix.
336: *
337:             CALL CLAQHB( UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED )
338:             RCEQU = LSAME( EQUED, 'Y' )
339:          END IF
340:       END IF
341: *
342: *     Scale the right-hand side.
343: *
344:       IF( RCEQU ) THEN
345:          DO 30 J = 1, NRHS
346:             DO 20 I = 1, N
347:                B( I, J ) = S( I )*B( I, J )
348:    20       CONTINUE
349:    30    CONTINUE
350:       END IF
351: *
352:       IF( NOFACT .OR. EQUIL ) THEN
353: *
354: *        Compute the Cholesky factorization A = U'*U or A = L*L'.
355: *
356:          IF( UPPER ) THEN
357:             DO 40 J = 1, N
358:                J1 = MAX( J-KD, 1 )
359:                CALL CCOPY( J-J1+1, AB( KD+1-J+J1, J ), 1,
360:      $                     AFB( KD+1-J+J1, J ), 1 )
361:    40       CONTINUE
362:          ELSE
363:             DO 50 J = 1, N
364:                J2 = MIN( J+KD, N )
365:                CALL CCOPY( J2-J+1, AB( 1, J ), 1, AFB( 1, J ), 1 )
366:    50       CONTINUE
367:          END IF
368: *
369:          CALL CPBTRF( UPLO, N, KD, AFB, LDAFB, INFO )
370: *
371: *        Return if INFO is non-zero.
372: *
373:          IF( INFO.GT.0 )THEN
374:             RCOND = ZERO
375:             RETURN
376:          END IF
377:       END IF
378: *
379: *     Compute the norm of the matrix A.
380: *
381:       ANORM = CLANHB( '1', UPLO, N, KD, AB, LDAB, RWORK )
382: *
383: *     Compute the reciprocal of the condition number of A.
384: *
385:       CALL CPBCON( UPLO, N, KD, AFB, LDAFB, ANORM, RCOND, WORK, RWORK,
386:      $             INFO )
387: *
388: *     Compute the solution matrix X.
389: *
390:       CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
391:       CALL CPBTRS( UPLO, N, KD, NRHS, AFB, LDAFB, X, LDX, INFO )
392: *
393: *     Use iterative refinement to improve the computed solution and
394: *     compute error bounds and backward error estimates for it.
395: *
396:       CALL CPBRFS( UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, B, LDB, X,
397:      $             LDX, FERR, BERR, WORK, RWORK, INFO )
398: *
399: *     Transform the solution matrix X to a solution of the original
400: *     system.
401: *
402:       IF( RCEQU ) THEN
403:          DO 70 J = 1, NRHS
404:             DO 60 I = 1, N
405:                X( I, J ) = S( I )*X( I, J )
406:    60       CONTINUE
407:    70    CONTINUE
408:          DO 80 J = 1, NRHS
409:             FERR( J ) = FERR( J ) / SCOND
410:    80    CONTINUE
411:       END IF
412: *
413: *     Set INFO = N+1 if the matrix is singular to working precision.
414: *
415:       IF( RCOND.LT.SLAMCH( 'Epsilon' ) )
416:      $   INFO = N + 1
417: *
418:       RETURN
419: *
420: *     End of CPBSVX
421: *
422:       END
423: