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