KTH framework for Nek5000 toolboxes; testing version  0.0.1
dgecon.f
Go to the documentation of this file.
1  SUBROUTINE dgecon( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK,
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 * February 29, 1992
8 *
9 * .. Scalar Arguments ..
10  CHARACTER NORM
11  INTEGER INFO, LDA, N
12  DOUBLE PRECISION ANORM, RCOND
13 * ..
14 * .. Array Arguments ..
15  INTEGER IWORK( * )
16  DOUBLE PRECISION A( LDA, * ), WORK( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * DGECON estimates the reciprocal of the condition number of a general
23 * real matrix A, in either the 1-norm or the infinity-norm, using
24 * the LU factorization computed by DGETRF.
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) DOUBLE PRECISION array, dimension (LDA,N)
43 * The factors L and U from the factorization A = P*L*U
44 * as computed by DGETRF.
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) DOUBLE PRECISION array, dimension (4*N)
58 *
59 * IWORK (workspace) INTEGER array, dimension (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 * ..
77 * .. External Functions ..
78  LOGICAL LSAME
79  INTEGER IDAMAX
80  DOUBLE PRECISION DLAMCH
81  EXTERNAL lsame, idamax, dlamch
82 * ..
83 * .. External Subroutines ..
84  EXTERNAL dlacon, dlatrs, drscl, xerbla
85 * ..
86 * .. Intrinsic Functions ..
87  INTRINSIC abs, max
88 * ..
89 * .. Executable Statements ..
90 *
91 * Test the input parameters.
92 *
93  info = 0
94  onenrm = norm.EQ.'1' .OR. lsame( norm, 'O' )
95  IF( .NOT.onenrm .AND. .NOT.lsame( norm, 'I' ) ) THEN
96  info = -1
97  ELSE IF( n.LT.0 ) THEN
98  info = -2
99  ELSE IF( lda.LT.max( 1, n ) ) THEN
100  info = -4
101  ELSE IF( anorm.LT.zero ) THEN
102  info = -5
103  END IF
104  IF( info.NE.0 ) THEN
105  CALL xerbla( 'DGECON', -info )
106  RETURN
107  END IF
108 *
109 * Quick return if possible
110 *
111  rcond = zero
112  IF( n.EQ.0 ) THEN
113  rcond = one
114  RETURN
115  ELSE IF( anorm.EQ.zero ) THEN
116  RETURN
117  END IF
118 *
119  smlnum = dlamch( 'Safe minimum' )
120 *
121 * Estimate the norm of inv(A).
122 *
123  ainvnm = zero
124  normin = 'N'
125  IF( onenrm ) THEN
126  kase1 = 1
127  ELSE
128  kase1 = 2
129  END IF
130  kase = 0
131  10 CONTINUE
132  CALL dlacon( n, work( n+1 ), work, iwork, ainvnm, kase )
133  IF( kase.NE.0 ) THEN
134  IF( kase.EQ.kase1 ) THEN
135 *
136 * Multiply by inv(L).
137 *
138  CALL dlatrs( 'Lower', 'No transpose', 'Unit', normin, n, a,
139  $ lda, work, sl, work( 2*n+1 ), info )
140 *
141 * Multiply by inv(U).
142 *
143  CALL dlatrs( 'Upper', 'No transpose', 'Non-unit', normin, n,
144  $ a, lda, work, su, work( 3*n+1 ), info )
145  ELSE
146 *
147 * Multiply by inv(U').
148 *
149  CALL dlatrs( 'Upper', 'Transpose', 'Non-unit', normin, n, a,
150  $ lda, work, su, work( 3*n+1 ), info )
151 *
152 * Multiply by inv(L').
153 *
154  CALL dlatrs( 'Lower', 'Transpose', 'Unit', normin, n, a,
155  $ lda, work, sl, work( 2*n+1 ), info )
156  END IF
157 *
158 * Divide X by 1/(SL*SU) if doing so will not cause overflow.
159 *
160  scale = sl*su
161  normin = 'Y'
162  IF( scale.NE.one ) THEN
163  ix = idamax( n, work, 1 )
164  IF( scale.LT.abs( work( ix ) )*smlnum .OR. scale.EQ.zero )
165  $ GO TO 20
166  CALL drscl( n, scale, work, 1 )
167  END IF
168  GO TO 10
169  END IF
170 *
171 * Compute the estimate of the reciprocal condition number.
172 *
173  IF( ainvnm.NE.zero )
174  $ rcond = ( one / ainvnm ) / anorm
175 *
176  20 CONTINUE
177  RETURN
178 *
179 * End of DGECON
180 *
181  END
subroutine dgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
Definition: dgecon.f:3
subroutine dlacon(N, V, X, ISGN, EST, KASE)
Definition: dlacon.f:2
subroutine dlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
Definition: dlatrs.f:3
subroutine drscl(N, SA, SX, INCX)
Definition: drscl.f:2
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2