Annotation of rpl/lapack/lapack/dgecon.f, revision 1.5
1.1 bertrand 1: SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK,
2: $ INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
10: *
11: * .. Scalar Arguments ..
12: CHARACTER NORM
13: INTEGER INFO, LDA, N
14: DOUBLE PRECISION ANORM, RCOND
15: * ..
16: * .. Array Arguments ..
17: INTEGER IWORK( * )
18: DOUBLE PRECISION A( LDA, * ), WORK( * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * DGECON estimates the reciprocal of the condition number of a general
25: * real matrix A, in either the 1-norm or the infinity-norm, using
26: * the LU factorization computed by DGETRF.
27: *
28: * An estimate is obtained for norm(inv(A)), and the reciprocal of the
29: * condition number is computed as
30: * RCOND = 1 / ( norm(A) * norm(inv(A)) ).
31: *
32: * Arguments
33: * =========
34: *
35: * NORM (input) CHARACTER*1
36: * Specifies whether the 1-norm condition number or the
37: * infinity-norm condition number is required:
38: * = '1' or 'O': 1-norm;
39: * = 'I': Infinity-norm.
40: *
41: * N (input) INTEGER
42: * The order of the matrix A. N >= 0.
43: *
44: * A (input) DOUBLE PRECISION array, dimension (LDA,N)
45: * The factors L and U from the factorization A = P*L*U
46: * as computed by DGETRF.
47: *
48: * LDA (input) INTEGER
49: * The leading dimension of the array A. LDA >= max(1,N).
50: *
51: * ANORM (input) DOUBLE PRECISION
52: * If NORM = '1' or 'O', the 1-norm of the original matrix A.
53: * If NORM = 'I', the infinity-norm of the original matrix A.
54: *
55: * RCOND (output) DOUBLE PRECISION
56: * The reciprocal of the condition number of the matrix A,
57: * computed as RCOND = 1/(norm(A) * norm(inv(A))).
58: *
59: * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
60: *
61: * IWORK (workspace) INTEGER array, dimension (N)
62: *
63: * INFO (output) INTEGER
64: * = 0: successful exit
65: * < 0: if INFO = -i, the i-th argument had an illegal value
66: *
67: * =====================================================================
68: *
69: * .. Parameters ..
70: DOUBLE PRECISION ONE, ZERO
71: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
72: * ..
73: * .. Local Scalars ..
74: LOGICAL ONENRM
75: CHARACTER NORMIN
76: INTEGER IX, KASE, KASE1
77: DOUBLE PRECISION AINVNM, SCALE, SL, SMLNUM, SU
78: * ..
79: * .. Local Arrays ..
80: INTEGER ISAVE( 3 )
81: * ..
82: * .. External Functions ..
83: LOGICAL LSAME
84: INTEGER IDAMAX
85: DOUBLE PRECISION DLAMCH
86: EXTERNAL LSAME, IDAMAX, DLAMCH
87: * ..
88: * .. External Subroutines ..
89: EXTERNAL DLACN2, DLATRS, DRSCL, XERBLA
90: * ..
91: * .. Intrinsic Functions ..
92: INTRINSIC ABS, MAX
93: * ..
94: * .. Executable Statements ..
95: *
96: * Test the input parameters.
97: *
98: INFO = 0
99: ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
100: IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
101: INFO = -1
102: ELSE IF( N.LT.0 ) THEN
103: INFO = -2
104: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
105: INFO = -4
106: ELSE IF( ANORM.LT.ZERO ) THEN
107: INFO = -5
108: END IF
109: IF( INFO.NE.0 ) THEN
110: CALL XERBLA( 'DGECON', -INFO )
111: RETURN
112: END IF
113: *
114: * Quick return if possible
115: *
116: RCOND = ZERO
117: IF( N.EQ.0 ) THEN
118: RCOND = ONE
119: RETURN
120: ELSE IF( ANORM.EQ.ZERO ) THEN
121: RETURN
122: END IF
123: *
124: SMLNUM = DLAMCH( 'Safe minimum' )
125: *
126: * Estimate the norm of inv(A).
127: *
128: AINVNM = ZERO
129: NORMIN = 'N'
130: IF( ONENRM ) THEN
131: KASE1 = 1
132: ELSE
133: KASE1 = 2
134: END IF
135: KASE = 0
136: 10 CONTINUE
137: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
138: IF( KASE.NE.0 ) THEN
139: IF( KASE.EQ.KASE1 ) THEN
140: *
141: * Multiply by inv(L).
142: *
143: CALL DLATRS( 'Lower', 'No transpose', 'Unit', NORMIN, N, A,
144: $ LDA, WORK, SL, WORK( 2*N+1 ), INFO )
145: *
146: * Multiply by inv(U).
147: *
148: CALL DLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
149: $ A, LDA, WORK, SU, WORK( 3*N+1 ), INFO )
150: ELSE
151: *
152: * Multiply by inv(U').
153: *
154: CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A,
155: $ LDA, WORK, SU, WORK( 3*N+1 ), INFO )
156: *
157: * Multiply by inv(L').
158: *
159: CALL DLATRS( 'Lower', 'Transpose', 'Unit', NORMIN, N, A,
160: $ LDA, WORK, SL, WORK( 2*N+1 ), INFO )
161: END IF
162: *
163: * Divide X by 1/(SL*SU) if doing so will not cause overflow.
164: *
165: SCALE = SL*SU
166: NORMIN = 'Y'
167: IF( SCALE.NE.ONE ) THEN
168: IX = IDAMAX( N, WORK, 1 )
169: IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
170: $ GO TO 20
171: CALL DRSCL( N, SCALE, WORK, 1 )
172: END IF
173: GO TO 10
174: END IF
175: *
176: * Compute the estimate of the reciprocal condition number.
177: *
178: IF( AINVNM.NE.ZERO )
179: $ RCOND = ( ONE / AINVNM ) / ANORM
180: *
181: 20 CONTINUE
182: RETURN
183: *
184: * End of DGECON
185: *
186: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgecon.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack