1 SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK,
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
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**T).
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**T).
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
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
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**T).
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**T).
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