 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DGBCON( NORM, N, KL, KU, AB, LDAB, IPIV, ANORM, RCOND, 00002 $ WORK, IWORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH. 00010 * 00011 * .. Scalar Arguments .. 00012 CHARACTER NORM 00013 INTEGER INFO, KL, KU, LDAB, N 00014 DOUBLE PRECISION ANORM, RCOND 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER IPIV( * ), IWORK( * ) 00018 DOUBLE PRECISION AB( LDAB, * ), WORK( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * DGBCON estimates the reciprocal of the condition number of a real 00025 * general band matrix A, in either the 1-norm or the infinity-norm, 00026 * using the LU factorization computed by DGBTRF. 00027 * 00028 * An estimate is obtained for norm(inv(A)), and the reciprocal of the 00029 * condition number is computed as 00030 * RCOND = 1 / ( norm(A) * norm(inv(A)) ). 00031 * 00032 * Arguments 00033 * ========= 00034 * 00035 * NORM (input) CHARACTER*1 00036 * Specifies whether the 1-norm condition number or the 00037 * infinity-norm condition number is required: 00038 * = '1' or 'O': 1-norm; 00039 * = 'I': Infinity-norm. 00040 * 00041 * N (input) INTEGER 00042 * The order of the matrix A. N >= 0. 00043 * 00044 * KL (input) INTEGER 00045 * The number of subdiagonals within the band of A. KL >= 0. 00046 * 00047 * KU (input) INTEGER 00048 * The number of superdiagonals within the band of A. KU >= 0. 00049 * 00050 * AB (input) DOUBLE PRECISION array, dimension (LDAB,N) 00051 * Details of the LU factorization of the band matrix A, as 00052 * computed by DGBTRF. U is stored as an upper triangular band 00053 * matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and 00054 * the multipliers used during the factorization are stored in 00055 * rows KL+KU+2 to 2*KL+KU+1. 00056 * 00057 * LDAB (input) INTEGER 00058 * The leading dimension of the array AB. LDAB >= 2*KL+KU+1. 00059 * 00060 * IPIV (input) INTEGER array, dimension (N) 00061 * The pivot indices; for 1 <= i <= N, row i of the matrix was 00062 * interchanged with row IPIV(i). 00063 * 00064 * ANORM (input) DOUBLE PRECISION 00065 * If NORM = '1' or 'O', the 1-norm of the original matrix A. 00066 * If NORM = 'I', the infinity-norm of the original matrix A. 00067 * 00068 * RCOND (output) DOUBLE PRECISION 00069 * The reciprocal of the condition number of the matrix A, 00070 * computed as RCOND = 1/(norm(A) * norm(inv(A))). 00071 * 00072 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N) 00073 * 00074 * IWORK (workspace) INTEGER array, dimension (N) 00075 * 00076 * INFO (output) INTEGER 00077 * = 0: successful exit 00078 * < 0: if INFO = -i, the i-th argument had an illegal value 00079 * 00080 * ===================================================================== 00081 * 00082 * .. Parameters .. 00083 DOUBLE PRECISION ONE, ZERO 00084 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00085 * .. 00086 * .. Local Scalars .. 00087 LOGICAL LNOTI, ONENRM 00088 CHARACTER NORMIN 00089 INTEGER IX, J, JP, KASE, KASE1, KD, LM 00090 DOUBLE PRECISION AINVNM, SCALE, SMLNUM, T 00091 * .. 00092 * .. Local Arrays .. 00093 INTEGER ISAVE( 3 ) 00094 * .. 00095 * .. External Functions .. 00096 LOGICAL LSAME 00097 INTEGER IDAMAX 00098 DOUBLE PRECISION DDOT, DLAMCH 00099 EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH 00100 * .. 00101 * .. External Subroutines .. 00102 EXTERNAL DAXPY, DLACN2, DLATBS, DRSCL, XERBLA 00103 * .. 00104 * .. Intrinsic Functions .. 00105 INTRINSIC ABS, MIN 00106 * .. 00107 * .. Executable Statements .. 00108 * 00109 * Test the input parameters. 00110 * 00111 INFO = 0 00112 ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' ) 00113 IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN 00114 INFO = -1 00115 ELSE IF( N.LT.0 ) THEN 00116 INFO = -2 00117 ELSE IF( KL.LT.0 ) THEN 00118 INFO = -3 00119 ELSE IF( KU.LT.0 ) THEN 00120 INFO = -4 00121 ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN 00122 INFO = -6 00123 ELSE IF( ANORM.LT.ZERO ) THEN 00124 INFO = -8 00125 END IF 00126 IF( INFO.NE.0 ) THEN 00127 CALL XERBLA( 'DGBCON', -INFO ) 00128 RETURN 00129 END IF 00130 * 00131 * Quick return if possible 00132 * 00133 RCOND = ZERO 00134 IF( N.EQ.0 ) THEN 00135 RCOND = ONE 00136 RETURN 00137 ELSE IF( ANORM.EQ.ZERO ) THEN 00138 RETURN 00139 END IF 00140 * 00141 SMLNUM = DLAMCH( 'Safe minimum' ) 00142 * 00143 * Estimate the norm of inv(A). 00144 * 00145 AINVNM = ZERO 00146 NORMIN = 'N' 00147 IF( ONENRM ) THEN 00148 KASE1 = 1 00149 ELSE 00150 KASE1 = 2 00151 END IF 00152 KD = KL + KU + 1 00153 LNOTI = KL.GT.0 00154 KASE = 0 00155 10 CONTINUE 00156 CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) 00157 IF( KASE.NE.0 ) THEN 00158 IF( KASE.EQ.KASE1 ) THEN 00159 * 00160 * Multiply by inv(L). 00161 * 00162 IF( LNOTI ) THEN 00163 DO 20 J = 1, N - 1 00164 LM = MIN( KL, N-J ) 00165 JP = IPIV( J ) 00166 T = WORK( JP ) 00167 IF( JP.NE.J ) THEN 00168 WORK( JP ) = WORK( J ) 00169 WORK( J ) = T 00170 END IF 00171 CALL DAXPY( LM, -T, AB( KD+1, J ), 1, WORK( J+1 ), 1 ) 00172 20 CONTINUE 00173 END IF 00174 * 00175 * Multiply by inv(U). 00176 * 00177 CALL DLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 00178 $ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ), 00179 $ INFO ) 00180 ELSE 00181 * 00182 * Multiply by inv(U**T). 00183 * 00184 CALL DLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, 00185 $ KL+KU, AB, LDAB, WORK, SCALE, WORK( 2*N+1 ), 00186 $ INFO ) 00187 * 00188 * Multiply by inv(L**T). 00189 * 00190 IF( LNOTI ) THEN 00191 DO 30 J = N - 1, 1, -1 00192 LM = MIN( KL, N-J ) 00193 WORK( J ) = WORK( J ) - DDOT( LM, AB( KD+1, J ), 1, 00194 $ WORK( J+1 ), 1 ) 00195 JP = IPIV( J ) 00196 IF( JP.NE.J ) THEN 00197 T = WORK( JP ) 00198 WORK( JP ) = WORK( J ) 00199 WORK( J ) = T 00200 END IF 00201 30 CONTINUE 00202 END IF 00203 END IF 00204 * 00205 * Divide X by 1/SCALE if doing so will not cause overflow. 00206 * 00207 NORMIN = 'Y' 00208 IF( SCALE.NE.ONE ) THEN 00209 IX = IDAMAX( N, WORK, 1 ) 00210 IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 00211 $ GO TO 40 00212 CALL DRSCL( N, SCALE, WORK, 1 ) 00213 END IF 00214 GO TO 10 00215 END IF 00216 * 00217 * Compute the estimate of the reciprocal condition number. 00218 * 00219 IF( AINVNM.NE.ZERO ) 00220 $ RCOND = ( ONE / AINVNM ) / ANORM 00221 * 00222 40 CONTINUE 00223 RETURN 00224 * 00225 * End of DGBCON 00226 * 00227 END
1.7.3 
