LAPACK 3.3.0

cptsvx.f

Go to the documentation of this file.
00001       SUBROUTINE CPTSVX( FACT, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
00002      $                   RCOND, FERR, BERR, WORK, RWORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          FACT
00011       INTEGER            INFO, LDB, LDX, N, NRHS
00012       REAL               RCOND
00013 *     ..
00014 *     .. Array Arguments ..
00015       REAL               BERR( * ), D( * ), DF( * ), FERR( * ),
00016      $                   RWORK( * )
00017       COMPLEX            B( LDB, * ), E( * ), EF( * ), WORK( * ),
00018      $                   X( LDX, * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  CPTSVX uses the factorization A = L*D*L**H to compute the solution
00025 *  to a complex system of linear equations A*X = B, where A is an
00026 *  N-by-N Hermitian positive definite tridiagonal matrix and X and B
00027 *  are N-by-NRHS matrices.
00028 *
00029 *  Error bounds on the solution and a condition estimate are also
00030 *  provided.
00031 *
00032 *  Description
00033 *  ===========
00034 *
00035 *  The following steps are performed:
00036 *
00037 *  1. If FACT = 'N', the matrix A is factored as A = L*D*L**H, where L
00038 *     is a unit lower bidiagonal matrix and D is diagonal.  The
00039 *     factorization can also be regarded as having the form
00040 *     A = U**H*D*U.
00041 *
00042 *  2. If the leading i-by-i principal minor is not positive definite,
00043 *     then the routine returns with INFO = i. Otherwise, the factored
00044 *     form of A is used to estimate the condition number of the matrix
00045 *     A.  If the reciprocal of the condition number is less than machine
00046 *     precision, INFO = N+1 is returned as a warning, but the routine
00047 *     still goes on to solve for X and compute error bounds as
00048 *     described below.
00049 *
00050 *  3. The system of equations is solved for X using the factored form
00051 *     of A.
00052 *
00053 *  4. Iterative refinement is applied to improve the computed solution
00054 *     matrix and calculate error bounds and backward error estimates
00055 *     for it.
00056 *
00057 *  Arguments
00058 *  =========
00059 *
00060 *  FACT    (input) CHARACTER*1
00061 *          Specifies whether or not the factored form of the matrix
00062 *          A is supplied on entry.
00063 *          = 'F':  On entry, DF and EF contain the factored form of A.
00064 *                  D, E, DF, and EF will not be modified.
00065 *          = 'N':  The matrix A will be copied to DF and EF and
00066 *                  factored.
00067 *
00068 *  N       (input) INTEGER
00069 *          The order of the matrix A.  N >= 0.
00070 *
00071 *  NRHS    (input) INTEGER
00072 *          The number of right hand sides, i.e., the number of columns
00073 *          of the matrices B and X.  NRHS >= 0.
00074 *
00075 *  D       (input) REAL array, dimension (N)
00076 *          The n diagonal elements of the tridiagonal matrix A.
00077 *
00078 *  E       (input) COMPLEX array, dimension (N-1)
00079 *          The (n-1) subdiagonal elements of the tridiagonal matrix A.
00080 *
00081 *  DF      (input or output) REAL array, dimension (N)
00082 *          If FACT = 'F', then DF is an input argument and on entry
00083 *          contains the n diagonal elements of the diagonal matrix D
00084 *          from the L*D*L**H factorization of A.
00085 *          If FACT = 'N', then DF is an output argument and on exit
00086 *          contains the n diagonal elements of the diagonal matrix D
00087 *          from the L*D*L**H factorization of A.
00088 *
00089 *  EF      (input or output) COMPLEX array, dimension (N-1)
00090 *          If FACT = 'F', then EF is an input argument and on entry
00091 *          contains the (n-1) subdiagonal elements of the unit
00092 *          bidiagonal factor L from the L*D*L**H factorization of A.
00093 *          If FACT = 'N', then EF is an output argument and on exit
00094 *          contains the (n-1) subdiagonal elements of the unit
00095 *          bidiagonal factor L from the L*D*L**H factorization of A.
00096 *
00097 *  B       (input) COMPLEX array, dimension (LDB,NRHS)
00098 *          The N-by-NRHS right hand side matrix B.
00099 *
00100 *  LDB     (input) INTEGER
00101 *          The leading dimension of the array B.  LDB >= max(1,N).
00102 *
00103 *  X       (output) COMPLEX array, dimension (LDX,NRHS)
00104 *          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X.
00105 *
00106 *  LDX     (input) INTEGER
00107 *          The leading dimension of the array X.  LDX >= max(1,N).
00108 *
00109 *  RCOND   (output) REAL
00110 *          The reciprocal condition number of the matrix A.  If RCOND
00111 *          is less than the machine precision (in particular, if
00112 *          RCOND = 0), the matrix is singular to working precision.
00113 *          This condition is indicated by a return code of INFO > 0.
00114 *
00115 *  FERR    (output) REAL array, dimension (NRHS)
00116 *          The forward error bound for each solution vector
00117 *          X(j) (the j-th column of the solution matrix X).
00118 *          If XTRUE is the true solution corresponding to X(j), FERR(j)
00119 *          is an estimated upper bound for the magnitude of the largest
00120 *          element in (X(j) - XTRUE) divided by the magnitude of the
00121 *          largest element in X(j).
00122 *
00123 *  BERR    (output) REAL array, dimension (NRHS)
00124 *          The componentwise relative backward error of each solution
00125 *          vector X(j) (i.e., the smallest relative change in any
00126 *          element of A or B that makes X(j) an exact solution).
00127 *
00128 *  WORK    (workspace) COMPLEX array, dimension (N)
00129 *
00130 *  RWORK   (workspace) REAL array, dimension (N)
00131 *
00132 *  INFO    (output) INTEGER
00133 *          = 0:  successful exit
00134 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00135 *          > 0:  if INFO = i, and i is
00136 *                <= N:  the leading minor of order i of A is
00137 *                       not positive definite, so the factorization
00138 *                       could not be completed, and the solution has not
00139 *                       been computed. RCOND = 0 is returned.
00140 *                = N+1: U is nonsingular, but RCOND is less than machine
00141 *                       precision, meaning that the matrix is singular
00142 *                       to working precision.  Nevertheless, the
00143 *                       solution and error bounds are computed because
00144 *                       there are a number of situations where the
00145 *                       computed solution can be more accurate than the
00146 *                       value of RCOND would suggest.
00147 *
00148 *  =====================================================================
00149 *
00150 *     .. Parameters ..
00151       REAL               ZERO
00152       PARAMETER          ( ZERO = 0.0E+0 )
00153 *     ..
00154 *     .. Local Scalars ..
00155       LOGICAL            NOFACT
00156       REAL               ANORM
00157 *     ..
00158 *     .. External Functions ..
00159       LOGICAL            LSAME
00160       REAL               CLANHT, SLAMCH
00161       EXTERNAL           LSAME, CLANHT, SLAMCH
00162 *     ..
00163 *     .. External Subroutines ..
00164       EXTERNAL           CCOPY, CLACPY, CPTCON, CPTRFS, CPTTRF, CPTTRS,
00165      $                   SCOPY, XERBLA
00166 *     ..
00167 *     .. Intrinsic Functions ..
00168       INTRINSIC          MAX
00169 *     ..
00170 *     .. Executable Statements ..
00171 *
00172 *     Test the input parameters.
00173 *
00174       INFO = 0
00175       NOFACT = LSAME( FACT, 'N' )
00176       IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
00177          INFO = -1
00178       ELSE IF( N.LT.0 ) THEN
00179          INFO = -2
00180       ELSE IF( NRHS.LT.0 ) THEN
00181          INFO = -3
00182       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00183          INFO = -9
00184       ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
00185          INFO = -11
00186       END IF
00187       IF( INFO.NE.0 ) THEN
00188          CALL XERBLA( 'CPTSVX', -INFO )
00189          RETURN
00190       END IF
00191 *
00192       IF( NOFACT ) THEN
00193 *
00194 *        Compute the L*D*L' (or U'*D*U) factorization of A.
00195 *
00196          CALL SCOPY( N, D, 1, DF, 1 )
00197          IF( N.GT.1 )
00198      $      CALL CCOPY( N-1, E, 1, EF, 1 )
00199          CALL CPTTRF( N, DF, EF, INFO )
00200 *
00201 *        Return if INFO is non-zero.
00202 *
00203          IF( INFO.GT.0 )THEN
00204             RCOND = ZERO
00205             RETURN
00206          END IF
00207       END IF
00208 *
00209 *     Compute the norm of the matrix A.
00210 *
00211       ANORM = CLANHT( '1', N, D, E )
00212 *
00213 *     Compute the reciprocal of the condition number of A.
00214 *
00215       CALL CPTCON( N, DF, EF, ANORM, RCOND, RWORK, INFO )
00216 *
00217 *     Compute the solution vectors X.
00218 *
00219       CALL CLACPY( 'Full', N, NRHS, B, LDB, X, LDX )
00220       CALL CPTTRS( 'Lower', N, NRHS, DF, EF, X, LDX, INFO )
00221 *
00222 *     Use iterative refinement to improve the computed solutions and
00223 *     compute error bounds and backward error estimates for them.
00224 *
00225       CALL CPTRFS( 'Lower', N, NRHS, D, E, DF, EF, B, LDB, X, LDX, FERR,
00226      $             BERR, WORK, RWORK, INFO )
00227 *
00228 *     Set INFO = N+1 if the matrix is singular to working precision.
00229 *
00230       IF( RCOND.LT.SLAMCH( 'Epsilon' ) )
00231      $   INFO = N + 1
00232 *
00233       RETURN
00234 *
00235 *     End of CPTSVX
00236 *
00237       END
 All Files Functions