KTH framework for Nek5000 toolboxes; testing version  0.0.1
zgecon.f
Go to the documentation of this file.
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
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2
subroutine zdrscl(N, SA, SX, INCX)
Definition: zdrscl.f:2
subroutine zgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
Definition: zgecon.f:3
subroutine zlacon(N, V, X, EST, KASE)
Definition: zlacon.f:2
subroutine zlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
Definition: zlatrs.f:3