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