LAPACK 3.3.1
Linear Algebra PACKage

zpbsvx.f

Go to the documentation of this file.
00001       SUBROUTINE ZPBSVX( FACT, UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB,
00002      $                   EQUED, S, B, LDB, X, LDX, RCOND, FERR, BERR,
00003      $                   WORK, RWORK, INFO )
00004 *
00005 *  -- LAPACK driver routine (version 3.3.1) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *  -- April 2011                                                      --
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          EQUED, FACT, UPLO
00012       INTEGER            INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS
00013       DOUBLE PRECISION   RCOND
00014 *     ..
00015 *     .. Array Arguments ..
00016       DOUBLE PRECISION   BERR( * ), FERR( * ), RWORK( * ), S( * )
00017       COMPLEX*16         AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
00018      $                   WORK( * ), X( LDX, * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  ZPBSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
00025 *  compute the solution to a complex system of linear equations
00026 *     A * X = B,
00027 *  where A is an N-by-N Hermitian positive definite band matrix and X
00028 *  and B are N-by-NRHS matrices.
00029 *
00030 *  Error bounds on the solution and a condition estimate are also
00031 *  provided.
00032 *
00033 *  Description
00034 *  ===========
00035 *
00036 *  The following steps are performed:
00037 *
00038 *  1. If FACT = 'E', real scaling factors are computed to equilibrate
00039 *     the system:
00040 *        diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
00041 *     Whether or not the system will be equilibrated depends on the
00042 *     scaling of the matrix A, but if equilibration is used, A is
00043 *     overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
00044 *
00045 *  2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
00046 *     factor the matrix A (after equilibration if FACT = 'E') as
00047 *        A = U**H * U,  if UPLO = 'U', or
00048 *        A = L * L**H,  if UPLO = 'L',
00049 *     where U is an upper triangular band matrix, and L is a lower
00050 *     triangular band matrix.
00051 *
00052 *  3. If the leading i-by-i principal minor is not positive definite,
00053 *     then the routine returns with INFO = i. Otherwise, the factored
00054 *     form of A is used to estimate the condition number of the matrix
00055 *     A.  If the reciprocal of the condition number is less than machine
00056 *     precision, INFO = N+1 is returned as a warning, but the routine
00057 *     still goes on to solve for X and compute error bounds as
00058 *     described below.
00059 *
00060 *  4. The system of equations is solved for X using the factored form
00061 *     of A.
00062 *
00063 *  5. Iterative refinement is applied to improve the computed solution
00064 *     matrix and calculate error bounds and backward error estimates
00065 *     for it.
00066 *
00067 *  6. If equilibration was used, the matrix X is premultiplied by
00068 *     diag(S) so that it solves the original system before
00069 *     equilibration.
00070 *
00071 *  Arguments
00072 *  =========
00073 *
00074 *  FACT    (input) CHARACTER*1
00075 *          Specifies whether or not the factored form of the matrix A is
00076 *          supplied on entry, and if not, whether the matrix A should be
00077 *          equilibrated before it is factored.
00078 *          = 'F':  On entry, AFB contains the factored form of A.
00079 *                  If EQUED = 'Y', the matrix A has been equilibrated
00080 *                  with scaling factors given by S.  AB and AFB will not
00081 *                  be modified.
00082 *          = 'N':  The matrix A will be copied to AFB and factored.
00083 *          = 'E':  The matrix A will be equilibrated if necessary, then
00084 *                  copied to AFB and factored.
00085 *
00086 *  UPLO    (input) CHARACTER*1
00087 *          = 'U':  Upper triangle of A is stored;
00088 *          = 'L':  Lower triangle of A is stored.
00089 *
00090 *  N       (input) INTEGER
00091 *          The number of linear equations, i.e., the order of the
00092 *          matrix A.  N >= 0.
00093 *
00094 *  KD      (input) INTEGER
00095 *          The number of superdiagonals of the matrix A if UPLO = 'U',
00096 *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
00097 *
00098 *  NRHS    (input) INTEGER
00099 *          The number of right-hand sides, i.e., the number of columns
00100 *          of the matrices B and X.  NRHS >= 0.
00101 *
00102 *  AB      (input/output) COMPLEX*16 array, dimension (LDAB,N)
00103 *          On entry, the upper or lower triangle of the Hermitian band
00104 *          matrix A, stored in the first KD+1 rows of the array, except
00105 *          if FACT = 'F' and EQUED = 'Y', then A must contain the
00106 *          equilibrated matrix diag(S)*A*diag(S).  The j-th column of A
00107 *          is stored in the j-th column of the array AB as follows:
00108 *          if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j;
00109 *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(N,j+KD).
00110 *          See below for further details.
00111 *
00112 *          On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
00113 *          diag(S)*A*diag(S).
00114 *
00115 *  LDAB    (input) INTEGER
00116 *          The leading dimension of the array A.  LDAB >= KD+1.
00117 *
00118 *  AFB     (input or output) COMPLEX*16 array, dimension (LDAFB,N)
00119 *          If FACT = 'F', then AFB is an input argument and on entry
00120 *          contains the triangular factor U or L from the Cholesky
00121 *          factorization A = U**H *U or A = L*L**H of the band matrix
00122 *          A, in the same storage format as A (see AB).  If EQUED = 'Y',
00123 *          then AFB is the factored form of the equilibrated matrix A.
00124 *
00125 *          If FACT = 'N', then AFB is an output argument and on exit
00126 *          returns the triangular factor U or L from the Cholesky
00127 *          factorization A = U**H *U or A = L*L**H.
00128 *
00129 *          If FACT = 'E', then AFB is an output argument and on exit
00130 *          returns the triangular factor U or L from the Cholesky
00131 *          factorization A = U**H *U or A = L*L**H of the equilibrated
00132 *          matrix A (see the description of A for the form of the
00133 *          equilibrated matrix).
00134 *
00135 *  LDAFB   (input) INTEGER
00136 *          The leading dimension of the array AFB.  LDAFB >= KD+1.
00137 *
00138 *  EQUED   (input or output) CHARACTER*1
00139 *          Specifies the form of equilibration that was done.
00140 *          = 'N':  No equilibration (always true if FACT = 'N').
00141 *          = 'Y':  Equilibration was done, i.e., A has been replaced by
00142 *                  diag(S) * A * diag(S).
00143 *          EQUED is an input argument if FACT = 'F'; otherwise, it is an
00144 *          output argument.
00145 *
00146 *  S       (input or output) DOUBLE PRECISION array, dimension (N)
00147 *          The scale factors for A; not accessed if EQUED = 'N'.  S is
00148 *          an input argument if FACT = 'F'; otherwise, S is an output
00149 *          argument.  If FACT = 'F' and EQUED = 'Y', each element of S
00150 *          must be positive.
00151 *
00152 *  B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
00153 *          On entry, the N-by-NRHS right hand side matrix B.
00154 *          On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
00155 *          B is overwritten by diag(S) * B.
00156 *
00157 *  LDB     (input) INTEGER
00158 *          The leading dimension of the array B.  LDB >= max(1,N).
00159 *
00160 *  X       (output) COMPLEX*16 array, dimension (LDX,NRHS)
00161 *          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
00162 *          the original system of equations.  Note that if EQUED = 'Y',
00163 *          A and B are modified on exit, and the solution to the
00164 *          equilibrated system is inv(diag(S))*X.
00165 *
00166 *  LDX     (input) INTEGER
00167 *          The leading dimension of the array X.  LDX >= max(1,N).
00168 *
00169 *  RCOND   (output) DOUBLE PRECISION
00170 *          The estimate of the reciprocal condition number of the matrix
00171 *          A after equilibration (if done).  If RCOND is less than the
00172 *          machine precision (in particular, if RCOND = 0), the matrix
00173 *          is singular to working precision.  This condition is
00174 *          indicated by a return code of INFO > 0.
00175 *
00176 *  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
00177 *          The estimated forward error bound for each solution vector
00178 *          X(j) (the j-th column of the solution matrix X).
00179 *          If XTRUE is the true solution corresponding to X(j), FERR(j)
00180 *          is an estimated upper bound for the magnitude of the largest
00181 *          element in (X(j) - XTRUE) divided by the magnitude of the
00182 *          largest element in X(j).  The estimate is as reliable as
00183 *          the estimate for RCOND, and is almost always a slight
00184 *          overestimate of the true error.
00185 *
00186 *  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
00187 *          The componentwise relative backward error of each solution
00188 *          vector X(j) (i.e., the smallest relative change in
00189 *          any element of A or B that makes X(j) an exact solution).
00190 *
00191 *  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
00192 *
00193 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00194 *
00195 *  INFO    (output) INTEGER
00196 *          = 0: successful exit
00197 *          < 0: if INFO = -i, the i-th argument had an illegal value
00198 *          > 0: if INFO = i, and i is
00199 *                <= N:  the leading minor of order i of A is
00200 *                       not positive definite, so the factorization
00201 *                       could not be completed, and the solution has not
00202 *                       been computed. RCOND = 0 is returned.
00203 *                = N+1: U is nonsingular, but RCOND is less than machine
00204 *                       precision, meaning that the matrix is singular
00205 *                       to working precision.  Nevertheless, the
00206 *                       solution and error bounds are computed because
00207 *                       there are a number of situations where the
00208 *                       computed solution can be more accurate than the
00209 *                       value of RCOND would suggest.
00210 *
00211 *  Further Details
00212 *  ===============
00213 *
00214 *  The band storage scheme is illustrated by the following example, when
00215 *  N = 6, KD = 2, and UPLO = 'U':
00216 *
00217 *  Two-dimensional storage of the Hermitian matrix A:
00218 *
00219 *     a11  a12  a13
00220 *          a22  a23  a24
00221 *               a33  a34  a35
00222 *                    a44  a45  a46
00223 *                         a55  a56
00224 *     (aij=conjg(aji))         a66
00225 *
00226 *  Band storage of the upper triangle of A:
00227 *
00228 *      *    *   a13  a24  a35  a46
00229 *      *   a12  a23  a34  a45  a56
00230 *     a11  a22  a33  a44  a55  a66
00231 *
00232 *  Similarly, if UPLO = 'L' the format of A is as follows:
00233 *
00234 *     a11  a22  a33  a44  a55  a66
00235 *     a21  a32  a43  a54  a65   *
00236 *     a31  a42  a53  a64   *    *
00237 *
00238 *  Array elements marked * are not used by the routine.
00239 *
00240 *  =====================================================================
00241 *
00242 *     .. Parameters ..
00243       DOUBLE PRECISION   ZERO, ONE
00244       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00245 *     ..
00246 *     .. Local Scalars ..
00247       LOGICAL            EQUIL, NOFACT, RCEQU, UPPER
00248       INTEGER            I, INFEQU, J, J1, J2
00249       DOUBLE PRECISION   AMAX, ANORM, BIGNUM, SCOND, SMAX, SMIN, SMLNUM
00250 *     ..
00251 *     .. External Functions ..
00252       LOGICAL            LSAME
00253       DOUBLE PRECISION   DLAMCH, ZLANHB
00254       EXTERNAL           LSAME, DLAMCH, ZLANHB
00255 *     ..
00256 *     .. External Subroutines ..
00257       EXTERNAL           XERBLA, ZCOPY, ZLACPY, ZLAQHB, ZPBCON, ZPBEQU,
00258      $                   ZPBRFS, ZPBTRF, ZPBTRS
00259 *     ..
00260 *     .. Intrinsic Functions ..
00261       INTRINSIC          MAX, MIN
00262 *     ..
00263 *     .. Executable Statements ..
00264 *
00265       INFO = 0
00266       NOFACT = LSAME( FACT, 'N' )
00267       EQUIL = LSAME( FACT, 'E' )
00268       UPPER = LSAME( UPLO, 'U' )
00269       IF( NOFACT .OR. EQUIL ) THEN
00270          EQUED = 'N'
00271          RCEQU = .FALSE.
00272       ELSE
00273          RCEQU = LSAME( EQUED, 'Y' )
00274          SMLNUM = DLAMCH( 'Safe minimum' )
00275          BIGNUM = ONE / SMLNUM
00276       END IF
00277 *
00278 *     Test the input parameters.
00279 *
00280       IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.LSAME( FACT, 'F' ) )
00281      $     THEN
00282          INFO = -1
00283       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00284          INFO = -2
00285       ELSE IF( N.LT.0 ) THEN
00286          INFO = -3
00287       ELSE IF( KD.LT.0 ) THEN
00288          INFO = -4
00289       ELSE IF( NRHS.LT.0 ) THEN
00290          INFO = -5
00291       ELSE IF( LDAB.LT.KD+1 ) THEN
00292          INFO = -7
00293       ELSE IF( LDAFB.LT.KD+1 ) THEN
00294          INFO = -9
00295       ELSE IF( LSAME( FACT, 'F' ) .AND. .NOT.
00296      $         ( RCEQU .OR. LSAME( EQUED, 'N' ) ) ) THEN
00297          INFO = -10
00298       ELSE
00299          IF( RCEQU ) THEN
00300             SMIN = BIGNUM
00301             SMAX = ZERO
00302             DO 10 J = 1, N
00303                SMIN = MIN( SMIN, S( J ) )
00304                SMAX = MAX( SMAX, S( J ) )
00305    10       CONTINUE
00306             IF( SMIN.LE.ZERO ) THEN
00307                INFO = -11
00308             ELSE IF( N.GT.0 ) THEN
00309                SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
00310             ELSE
00311                SCOND = ONE
00312             END IF
00313          END IF
00314          IF( INFO.EQ.0 ) THEN
00315             IF( LDB.LT.MAX( 1, N ) ) THEN
00316                INFO = -13
00317             ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00318                INFO = -15
00319             END IF
00320          END IF
00321       END IF
00322 *
00323       IF( INFO.NE.0 ) THEN
00324          CALL XERBLA( 'ZPBSVX', -INFO )
00325          RETURN
00326       END IF
00327 *
00328       IF( EQUIL ) THEN
00329 *
00330 *        Compute row and column scalings to equilibrate the matrix A.
00331 *
00332          CALL ZPBEQU( UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, INFEQU )
00333          IF( INFEQU.EQ.0 ) THEN
00334 *
00335 *           Equilibrate the matrix.
00336 *
00337             CALL ZLAQHB( UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED )
00338             RCEQU = LSAME( EQUED, 'Y' )
00339          END IF
00340       END IF
00341 *
00342 *     Scale the right-hand side.
00343 *
00344       IF( RCEQU ) THEN
00345          DO 30 J = 1, NRHS
00346             DO 20 I = 1, N
00347                B( I, J ) = S( I )*B( I, J )
00348    20       CONTINUE
00349    30    CONTINUE
00350       END IF
00351 *
00352       IF( NOFACT .OR. EQUIL ) THEN
00353 *
00354 *        Compute the Cholesky factorization A = U**H *U or A = L*L**H.
00355 *
00356          IF( UPPER ) THEN
00357             DO 40 J = 1, N
00358                J1 = MAX( J-KD, 1 )
00359                CALL ZCOPY( J-J1+1, AB( KD+1-J+J1, J ), 1,
00360      $                     AFB( KD+1-J+J1, J ), 1 )
00361    40       CONTINUE
00362          ELSE
00363             DO 50 J = 1, N
00364                J2 = MIN( J+KD, N )
00365                CALL ZCOPY( J2-J+1, AB( 1, J ), 1, AFB( 1, J ), 1 )
00366    50       CONTINUE
00367          END IF
00368 *
00369          CALL ZPBTRF( UPLO, N, KD, AFB, LDAFB, INFO )
00370 *
00371 *        Return if INFO is non-zero.
00372 *
00373          IF( INFO.GT.0 )THEN
00374             RCOND = ZERO
00375             RETURN
00376          END IF
00377       END IF
00378 *
00379 *     Compute the norm of the matrix A.
00380 *
00381       ANORM = ZLANHB( '1', UPLO, N, KD, AB, LDAB, RWORK )
00382 *
00383 *     Compute the reciprocal of the condition number of A.
00384 *
00385       CALL ZPBCON( UPLO, N, KD, AFB, LDAFB, ANORM, RCOND, WORK, RWORK,
00386      $             INFO )
00387 *
00388 *     Compute the solution matrix X.
00389 *
00390       CALL ZLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
00391       CALL ZPBTRS( UPLO, N, KD, NRHS, AFB, LDAFB, X, LDX, INFO )
00392 *
00393 *     Use iterative refinement to improve the computed solution and
00394 *     compute error bounds and backward error estimates for it.
00395 *
00396       CALL ZPBRFS( UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, B, LDB, X,
00397      $             LDX, FERR, BERR, WORK, RWORK, INFO )
00398 *
00399 *     Transform the solution matrix X to a solution of the original
00400 *     system.
00401 *
00402       IF( RCEQU ) THEN
00403          DO 70 J = 1, NRHS
00404             DO 60 I = 1, N
00405                X( I, J ) = S( I )*X( I, J )
00406    60       CONTINUE
00407    70    CONTINUE
00408          DO 80 J = 1, NRHS
00409             FERR( J ) = FERR( J ) / SCOND
00410    80    CONTINUE
00411       END IF
00412 *
00413 *     Set INFO = N+1 if the matrix is singular to working precision.
00414 *
00415       IF( RCOND.LT.DLAMCH( 'Epsilon' ) )
00416      $   INFO = N + 1
00417 *
00418       RETURN
00419 *
00420 *     End of ZPBSVX
00421 *
00422       END
 All Files Functions