![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dla_gercond.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dla_gercond.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: DOUBLE PRECISION FUNCTION DLA_GERCOND ( TRANS, N, A, LDA, AF,
2: $ LDAF, IPIV, CMODE, C,
3: $ INFO, WORK, IWORK )
4: *
5: * -- LAPACK routine (version 3.2.1) --
6: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
7: * -- Jason Riedy of Univ. of California Berkeley. --
8: * -- April 2009 --
9: *
10: * -- LAPACK is a software package provided by Univ. of Tennessee, --
11: * -- Univ. of California Berkeley and NAG Ltd. --
12: *
13: IMPLICIT NONE
14: * ..
15: * .. Scalar Arguments ..
16: CHARACTER TRANS
17: INTEGER N, LDA, LDAF, INFO, CMODE
18: * ..
19: * .. Array Arguments ..
20: INTEGER IPIV( * ), IWORK( * )
21: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ),
22: $ C( * )
23: * ..
24: *
25: * Purpose
26: * =======
27: *
28: * DLA_GERCOND estimates the Skeel condition number of op(A) * op2(C)
29: * where op2 is determined by CMODE as follows
30: * CMODE = 1 op2(C) = C
31: * CMODE = 0 op2(C) = I
32: * CMODE = -1 op2(C) = inv(C)
33: * The Skeel condition number cond(A) = norminf( |inv(A)||A| )
34: * is computed by computing scaling factors R such that
35: * diag(R)*A*op2(C) is row equilibrated and computing the standard
36: * infinity-norm condition number.
37: *
38: * Arguments
39: * ==========
40: *
41: * TRANS (input) CHARACTER*1
42: * Specifies the form of the system of equations:
43: * = 'N': A * X = B (No transpose)
44: * = 'T': A**T * X = B (Transpose)
45: * = 'C': A**H * X = B (Conjugate Transpose = Transpose)
46: *
47: * N (input) INTEGER
48: * The number of linear equations, i.e., the order of the
49: * matrix A. N >= 0.
50: *
51: * A (input) DOUBLE PRECISION array, dimension (LDA,N)
52: * On entry, the N-by-N matrix A.
53: *
54: * LDA (input) INTEGER
55: * The leading dimension of the array A. LDA >= max(1,N).
56: *
57: * AF (input) DOUBLE PRECISION array, dimension (LDAF,N)
58: * The factors L and U from the factorization
59: * A = P*L*U as computed by DGETRF.
60: *
61: * LDAF (input) INTEGER
62: * The leading dimension of the array AF. LDAF >= max(1,N).
63: *
64: * IPIV (input) INTEGER array, dimension (N)
65: * The pivot indices from the factorization A = P*L*U
66: * as computed by DGETRF; row i of the matrix was interchanged
67: * with row IPIV(i).
68: *
69: * CMODE (input) INTEGER
70: * Determines op2(C) in the formula op(A) * op2(C) as follows:
71: * CMODE = 1 op2(C) = C
72: * CMODE = 0 op2(C) = I
73: * CMODE = -1 op2(C) = inv(C)
74: *
75: * C (input) DOUBLE PRECISION array, dimension (N)
76: * The vector C in the formula op(A) * op2(C).
77: *
78: * INFO (output) INTEGER
79: * = 0: Successful exit.
80: * i > 0: The ith argument is invalid.
81: *
82: * WORK (input) DOUBLE PRECISION array, dimension (3*N).
83: * Workspace.
84: *
85: * IWORK (input) INTEGER array, dimension (N).
86: * Workspace.
87: *
88: * =====================================================================
89: *
90: * .. Local Scalars ..
91: LOGICAL NOTRANS
92: INTEGER KASE, I, J
93: DOUBLE PRECISION AINVNM, TMP
94: * ..
95: * .. Local Arrays ..
96: INTEGER ISAVE( 3 )
97: * ..
98: * .. External Functions ..
99: LOGICAL LSAME
100: EXTERNAL LSAME
101: * ..
102: * .. External Subroutines ..
103: EXTERNAL DLACN2, DGETRS, XERBLA
104: * ..
105: * .. Intrinsic Functions ..
106: INTRINSIC ABS, MAX
107: * ..
108: * .. Executable Statements ..
109: *
110: DLA_GERCOND = 0.0D+0
111: *
112: INFO = 0
113: NOTRANS = LSAME( TRANS, 'N' )
114: IF ( .NOT. NOTRANS .AND. .NOT. LSAME(TRANS, 'T')
115: $ .AND. .NOT. LSAME(TRANS, 'C') ) THEN
116: INFO = -1
117: ELSE IF( N.LT.0 ) THEN
118: INFO = -2
119: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
120: INFO = -4
121: ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
122: INFO = -6
123: END IF
124: IF( INFO.NE.0 ) THEN
125: CALL XERBLA( 'DLA_GERCOND', -INFO )
126: RETURN
127: END IF
128: IF( N.EQ.0 ) THEN
129: DLA_GERCOND = 1.0D+0
130: RETURN
131: END IF
132: *
133: * Compute the equilibration matrix R such that
134: * inv(R)*A*C has unit 1-norm.
135: *
136: IF (NOTRANS) THEN
137: DO I = 1, N
138: TMP = 0.0D+0
139: IF ( CMODE .EQ. 1 ) THEN
140: DO J = 1, N
141: TMP = TMP + ABS( A( I, J ) * C( J ) )
142: END DO
143: ELSE IF ( CMODE .EQ. 0 ) THEN
144: DO J = 1, N
145: TMP = TMP + ABS( A( I, J ) )
146: END DO
147: ELSE
148: DO J = 1, N
149: TMP = TMP + ABS( A( I, J ) / C( J ) )
150: END DO
151: END IF
152: WORK( 2*N+I ) = TMP
153: END DO
154: ELSE
155: DO I = 1, N
156: TMP = 0.0D+0
157: IF ( CMODE .EQ. 1 ) THEN
158: DO J = 1, N
159: TMP = TMP + ABS( A( J, I ) * C( J ) )
160: END DO
161: ELSE IF ( CMODE .EQ. 0 ) THEN
162: DO J = 1, N
163: TMP = TMP + ABS( A( J, I ) )
164: END DO
165: ELSE
166: DO J = 1, N
167: TMP = TMP + ABS( A( J, I ) / C( J ) )
168: END DO
169: END IF
170: WORK( 2*N+I ) = TMP
171: END DO
172: END IF
173: *
174: * Estimate the norm of inv(op(A)).
175: *
176: AINVNM = 0.0D+0
177:
178: KASE = 0
179: 10 CONTINUE
180: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
181: IF( KASE.NE.0 ) THEN
182: IF( KASE.EQ.2 ) THEN
183: *
184: * Multiply by R.
185: *
186: DO I = 1, N
187: WORK(I) = WORK(I) * WORK(2*N+I)
188: END DO
189:
190: IF (NOTRANS) THEN
191: CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
192: $ WORK, N, INFO )
193: ELSE
194: CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV,
195: $ WORK, N, INFO )
196: END IF
197: *
198: * Multiply by inv(C).
199: *
200: IF ( CMODE .EQ. 1 ) THEN
201: DO I = 1, N
202: WORK( I ) = WORK( I ) / C( I )
203: END DO
204: ELSE IF ( CMODE .EQ. -1 ) THEN
205: DO I = 1, N
206: WORK( I ) = WORK( I ) * C( I )
207: END DO
208: END IF
209: ELSE
210: *
211: * Multiply by inv(C').
212: *
213: IF ( CMODE .EQ. 1 ) THEN
214: DO I = 1, N
215: WORK( I ) = WORK( I ) / C( I )
216: END DO
217: ELSE IF ( CMODE .EQ. -1 ) THEN
218: DO I = 1, N
219: WORK( I ) = WORK( I ) * C( I )
220: END DO
221: END IF
222:
223: IF (NOTRANS) THEN
224: CALL DGETRS( 'Transpose', N, 1, AF, LDAF, IPIV,
225: $ WORK, N, INFO )
226: ELSE
227: CALL DGETRS( 'No transpose', N, 1, AF, LDAF, IPIV,
228: $ WORK, N, INFO )
229: END IF
230: *
231: * Multiply by R.
232: *
233: DO I = 1, N
234: WORK( I ) = WORK( I ) * WORK( 2*N+I )
235: END DO
236: END IF
237: GO TO 10
238: END IF
239: *
240: * Compute the estimate of the reciprocal condition number.
241: *
242: IF( AINVNM .NE. 0.0D+0 )
243: $ DLA_GERCOND = ( 1.0D+0 / AINVNM )
244: *
245: RETURN
246: *
247: END