1 SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, 2 $ INFO ) 3* 4* -- LAPACK routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* March 31, 1993 8* 9* .. Scalar Arguments .. 10 CHARACTER NORM 11 INTEGER INFO, LDA, N 12 DOUBLE PRECISION ANORM, RCOND 13* .. 14* .. Array Arguments .. 15 DOUBLE PRECISION RWORK( * ) 16 COMPLEX*16 A( LDA, * ), WORK( * ) 17* .. 18* 19* Purpose 20* ======= 21* 22* ZGECON estimates the reciprocal of the condition number of a general 23* complex matrix A, in either the 1-norm or the infinity-norm, using 24* the LU factorization computed by ZGETRF. 25* 26* An estimate is obtained for norm(inv(A)), and the reciprocal of the 27* condition number is computed as 28* RCOND = 1 / ( norm(A) * norm(inv(A)) ). 29* 30* Arguments 31* ========= 32* 33* NORM (input) CHARACTER*1 34* Specifies whether the 1-norm condition number or the 35* infinity-norm condition number is required: 36* = '1' or 'O': 1-norm; 37* = 'I': Infinity-norm. 38* 39* N (input) INTEGER 40* The order of the matrix A. N >= 0. 41* 42* A (input) COMPLEX*16 array, dimension (LDA,N) 43* The factors L and U from the factorization A = P*L*U 44* as computed by ZGETRF. 45* 46* LDA (input) INTEGER 47* The leading dimension of the array A. LDA >= max(1,N). 48* 49* ANORM (input) DOUBLE PRECISION 50* If NORM = '1' or 'O', the 1-norm of the original matrix A. 51* If NORM = 'I', the infinity-norm of the original matrix A. 52* 53* RCOND (output) DOUBLE PRECISION 54* The reciprocal of the condition number of the matrix A, 55* computed as RCOND = 1/(norm(A) * norm(inv(A))). 56* 57* WORK (workspace) COMPLEX*16 array, dimension (2*N) 58* 59* RWORK (workspace) DOUBLE PRECISION array, dimension (2*N) 60* 61* INFO (output) INTEGER 62* = 0: successful exit 63* < 0: if INFO = -i, the i-th argument had an illegal value 64* 65* ===================================================================== 66* 67* .. Parameters .. 68 DOUBLE PRECISION ONE, ZERO 69 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 70* .. 71* .. Local Scalars .. 72 LOGICAL ONENRM 73 CHARACTER NORMIN 74 INTEGER IX, KASE, KASE1 75 DOUBLE PRECISION AINVNM, SCALE, SL, SMLNUM, SU 76 COMPLEX*16 ZDUM 77* .. 78* .. External Functions .. 79 LOGICAL LSAME 80 INTEGER IZAMAX 81 DOUBLE PRECISION DLAMCH 82 EXTERNAL LSAME, IZAMAX, DLAMCH 83* .. 84* .. External Subroutines .. 85 EXTERNAL XERBLA, ZDRSCL, ZLACON, ZLATRS 86* .. 87* .. Intrinsic Functions .. 88 INTRINSIC ABS, DBLE, DIMAG, MAX 89* .. 90* .. Statement Functions .. 91 DOUBLE PRECISION CABS1 92* .. 93* .. Statement Function definitions .. 94 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 95* .. 96* .. Executable Statements .. 97* 98* Test the input parameters. 99* 100 INFO = 0 101 ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' ) 102 IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN 103 INFO = -1 104 ELSE IF( N.LT.0 ) THEN 105 INFO = -2 106 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 107 INFO = -4 108 ELSE IF( ANORM.LT.ZERO ) THEN 109 INFO = -5 110 END IF 111 IF( INFO.NE.0 ) THEN 112 CALL XERBLA( 'ZGECON', -INFO ) 113 RETURN 114 END IF 115* 116* Quick return if possible 117* 118 RCOND = ZERO 119 IF( N.EQ.0 ) THEN 120 RCOND = ONE 121 RETURN 122 ELSE IF( ANORM.EQ.ZERO ) THEN 123 RETURN 124 END IF 125* 126 SMLNUM = DLAMCH( 'Safe minimum' ) 127* 128* Estimate the norm of inv(A). 129* 130 AINVNM = ZERO 131 NORMIN = 'N' 132 IF( ONENRM ) THEN 133 KASE1 = 1 134 ELSE 135 KASE1 = 2 136 END IF 137 KASE = 0 138 10 CONTINUE 139 CALL ZLACON( N, WORK( N+1 ), WORK, AINVNM, KASE ) 140 IF( KASE.NE.0 ) THEN 141 IF( KASE.EQ.KASE1 ) THEN 142* 143* Multiply by inv(L). 144* 145 CALL ZLATRS( 'Lower', 'No transpose', 'Unit', NORMIN, N, A, 146 $ LDA, WORK, SL, RWORK, INFO ) 147* 148* Multiply by inv(U). 149* 150 CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, 151 $ A, LDA, WORK, SU, RWORK( N+1 ), INFO ) 152 ELSE 153* 154* Multiply by inv(U'). 155* 156 CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit', 157 $ NORMIN, N, A, LDA, WORK, SU, RWORK( N+1 ), 158 $ INFO ) 159* 160* Multiply by inv(L'). 161* 162 CALL ZLATRS( 'Lower', 'Conjugate transpose', 'Unit', NORMIN, 163 $ N, A, LDA, WORK, SL, RWORK, INFO ) 164 END IF 165* 166* Divide X by 1/(SL*SU) if doing so will not cause overflow. 167* 168 SCALE = SL*SU 169 NORMIN = 'Y' 170 IF( SCALE.NE.ONE ) THEN 171 IX = IZAMAX( N, WORK, 1 ) 172 IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) 173 $ GO TO 20 174 CALL ZDRSCL( N, SCALE, WORK, 1 ) 175 END IF 176 GO TO 10 177 END IF 178* 179* Compute the estimate of the reciprocal condition number. 180* 181 IF( AINVNM.NE.ZERO ) 182 $ RCOND = ( ONE / AINVNM ) / ANORM 183* 184 20 CONTINUE 185 RETURN 186* 187* End of ZGECON 188* 189 END 190