Mercurial > octave-nkf
comparison libcruft/lapack/zgecon.f @ 4329:d53c33d93440
[project @ 2003-02-18 20:00:48 by jwe]
author | jwe |
---|---|
date | Tue, 18 Feb 2003 20:08:20 +0000 |
parents | |
children | 68db500cb558 |
comparison
equal
deleted
inserted
replaced
4328:f7b63f362168 | 4329:d53c33d93440 |
---|---|
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 |