![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dla_porcond.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_porcond.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: DOUBLE PRECISION FUNCTION DLA_PORCOND( UPLO, N, A, LDA, AF, LDAF,
2: $ CMODE, C, INFO, WORK,
3: $ IWORK )
4: *
5: * -- LAPACK routine (version 3.2.2) --
6: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
7: * -- Jason Riedy of Univ. of California Berkeley. --
8: * -- June 2010 --
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 UPLO
17: INTEGER N, LDA, LDAF, INFO, CMODE
18: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * ),
19: $ C( * )
20: * ..
21: * .. Array Arguments ..
22: INTEGER IWORK( * )
23: * ..
24: *
25: * Purpose
26: * =======
27: *
28: * DLA_PORCOND 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: * UPLO (input) CHARACTER*1
42: * = 'U': Upper triangle of A is stored;
43: * = 'L': Lower triangle of A is stored.
44: *
45: * N (input) INTEGER
46: * The number of linear equations, i.e., the order of the
47: * matrix A. N >= 0.
48: *
49: * A (input) DOUBLE PRECISION array, dimension (LDA,N)
50: * On entry, the N-by-N matrix A.
51: *
52: * LDA (input) INTEGER
53: * The leading dimension of the array A. LDA >= max(1,N).
54: *
55: * AF (input) DOUBLE PRECISION array, dimension (LDAF,N)
56: * The triangular factor U or L from the Cholesky factorization
57: * A = U**T*U or A = L*L**T, as computed by DPOTRF.
58: *
59: * LDAF (input) INTEGER
60: * The leading dimension of the array AF. LDAF >= max(1,N).
61: *
62: * CMODE (input) INTEGER
63: * Determines op2(C) in the formula op(A) * op2(C) as follows:
64: * CMODE = 1 op2(C) = C
65: * CMODE = 0 op2(C) = I
66: * CMODE = -1 op2(C) = inv(C)
67: *
68: * C (input) DOUBLE PRECISION array, dimension (N)
69: * The vector C in the formula op(A) * op2(C).
70: *
71: * INFO (output) INTEGER
72: * = 0: Successful exit.
73: * i > 0: The ith argument is invalid.
74: *
75: * WORK (input) DOUBLE PRECISION array, dimension (3*N).
76: * Workspace.
77: *
78: * IWORK (input) INTEGER array, dimension (N).
79: * Workspace.
80: *
81: * =====================================================================
82: *
83: * .. Local Scalars ..
84: INTEGER KASE, I, J
85: DOUBLE PRECISION AINVNM, TMP
86: LOGICAL UP
87: * ..
88: * .. Array Arguments ..
89: INTEGER ISAVE( 3 )
90: * ..
91: * .. External Functions ..
92: LOGICAL LSAME
93: INTEGER IDAMAX
94: EXTERNAL LSAME, IDAMAX
95: * ..
96: * .. External Subroutines ..
97: EXTERNAL DLACN2, DPOTRS, XERBLA
98: * ..
99: * .. Intrinsic Functions ..
100: INTRINSIC ABS, MAX
101: * ..
102: * .. Executable Statements ..
103: *
104: DLA_PORCOND = 0.0D+0
105: *
106: INFO = 0
107: IF( N.LT.0 ) THEN
108: INFO = -2
109: END IF
110: IF( INFO.NE.0 ) THEN
111: CALL XERBLA( 'DLA_PORCOND', -INFO )
112: RETURN
113: END IF
114:
115: IF( N.EQ.0 ) THEN
116: DLA_PORCOND = 1.0D+0
117: RETURN
118: END IF
119: UP = .FALSE.
120: IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
121: *
122: * Compute the equilibration matrix R such that
123: * inv(R)*A*C has unit 1-norm.
124: *
125: IF ( UP ) THEN
126: DO I = 1, N
127: TMP = 0.0D+0
128: IF ( CMODE .EQ. 1 ) THEN
129: DO J = 1, I
130: TMP = TMP + ABS( A( J, I ) * C( J ) )
131: END DO
132: DO J = I+1, N
133: TMP = TMP + ABS( A( I, J ) * C( J ) )
134: END DO
135: ELSE IF ( CMODE .EQ. 0 ) THEN
136: DO J = 1, I
137: TMP = TMP + ABS( A( J, I ) )
138: END DO
139: DO J = I+1, N
140: TMP = TMP + ABS( A( I, J ) )
141: END DO
142: ELSE
143: DO J = 1, I
144: TMP = TMP + ABS( A( J ,I ) / C( J ) )
145: END DO
146: DO J = I+1, N
147: TMP = TMP + ABS( A( I, J ) / C( J ) )
148: END DO
149: END IF
150: WORK( 2*N+I ) = TMP
151: END DO
152: ELSE
153: DO I = 1, N
154: TMP = 0.0D+0
155: IF ( CMODE .EQ. 1 ) THEN
156: DO J = 1, I
157: TMP = TMP + ABS( A( I, J ) * C( J ) )
158: END DO
159: DO J = I+1, N
160: TMP = TMP + ABS( A( J, I ) * C( J ) )
161: END DO
162: ELSE IF ( CMODE .EQ. 0 ) THEN
163: DO J = 1, I
164: TMP = TMP + ABS( A( I, J ) )
165: END DO
166: DO J = I+1, N
167: TMP = TMP + ABS( A( J, I ) )
168: END DO
169: ELSE
170: DO J = 1, I
171: TMP = TMP + ABS( A( I, J ) / C( J ) )
172: END DO
173: DO J = I+1, N
174: TMP = TMP + ABS( A( J, I ) / C( J ) )
175: END DO
176: END IF
177: WORK( 2*N+I ) = TMP
178: END DO
179: ENDIF
180: *
181: * Estimate the norm of inv(op(A)).
182: *
183: AINVNM = 0.0D+0
184:
185: KASE = 0
186: 10 CONTINUE
187: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
188: IF( KASE.NE.0 ) THEN
189: IF( KASE.EQ.2 ) THEN
190: *
191: * Multiply by R.
192: *
193: DO I = 1, N
194: WORK( I ) = WORK( I ) * WORK( 2*N+I )
195: END DO
196:
197: IF (UP) THEN
198: CALL DPOTRS( 'Upper', N, 1, AF, LDAF, WORK, N, INFO )
199: ELSE
200: CALL DPOTRS( 'Lower', N, 1, AF, LDAF, WORK, N, INFO )
201: ENDIF
202: *
203: * Multiply by inv(C).
204: *
205: IF ( CMODE .EQ. 1 ) THEN
206: DO I = 1, N
207: WORK( I ) = WORK( I ) / C( I )
208: END DO
209: ELSE IF ( CMODE .EQ. -1 ) THEN
210: DO I = 1, N
211: WORK( I ) = WORK( I ) * C( I )
212: END DO
213: END IF
214: ELSE
215: *
216: * Multiply by inv(C').
217: *
218: IF ( CMODE .EQ. 1 ) THEN
219: DO I = 1, N
220: WORK( I ) = WORK( I ) / C( I )
221: END DO
222: ELSE IF ( CMODE .EQ. -1 ) THEN
223: DO I = 1, N
224: WORK( I ) = WORK( I ) * C( I )
225: END DO
226: END IF
227:
228: IF ( UP ) THEN
229: CALL DPOTRS( 'Upper', N, 1, AF, LDAF, WORK, N, INFO )
230: ELSE
231: CALL DPOTRS( 'Lower', N, 1, AF, LDAF, WORK, N, INFO )
232: ENDIF
233: *
234: * Multiply by R.
235: *
236: DO I = 1, N
237: WORK( I ) = WORK( I ) * WORK( 2*N+I )
238: END DO
239: END IF
240: GO TO 10
241: END IF
242: *
243: * Compute the estimate of the reciprocal condition number.
244: *
245: IF( AINVNM .NE. 0.0D+0 )
246: $ DLA_PORCOND = ( 1.0D+0 / AINVNM )
247: *
248: RETURN
249: *
250: END