1: SUBROUTINE DTPCON( NORM, UPLO, DIAG, N, AP, 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 DIAG, NORM, UPLO
13: INTEGER INFO, N
14: DOUBLE PRECISION RCOND
15: * ..
16: * .. Array Arguments ..
17: INTEGER IWORK( * )
18: DOUBLE PRECISION AP( * ), WORK( * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * DTPCON estimates the reciprocal of the condition number of a packed
25: * triangular matrix A, in either the 1-norm or the infinity-norm.
26: *
27: * The norm of A is computed and an estimate is obtained for
28: * norm(inv(A)), then the reciprocal of the condition number is
29: * 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: * UPLO (input) CHARACTER*1
42: * = 'U': A is upper triangular;
43: * = 'L': A is lower triangular.
44: *
45: * DIAG (input) CHARACTER*1
46: * = 'N': A is non-unit triangular;
47: * = 'U': A is unit triangular.
48: *
49: * N (input) INTEGER
50: * The order of the matrix A. N >= 0.
51: *
52: * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
53: * The upper or lower triangular matrix A, packed columnwise in
54: * a linear array. The j-th column of A is stored in the array
55: * AP as follows:
56: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
57: * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
58: * If DIAG = 'U', the diagonal elements of A are not referenced
59: * and are assumed to be 1.
60: *
61: * RCOND (output) DOUBLE PRECISION
62: * The reciprocal of the condition number of the matrix A,
63: * computed as RCOND = 1/(norm(A) * norm(inv(A))).
64: *
65: * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
66: *
67: * IWORK (workspace) INTEGER array, dimension (N)
68: *
69: * INFO (output) INTEGER
70: * = 0: successful exit
71: * < 0: if INFO = -i, the i-th argument had an illegal value
72: *
73: * =====================================================================
74: *
75: * .. Parameters ..
76: DOUBLE PRECISION ONE, ZERO
77: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
78: * ..
79: * .. Local Scalars ..
80: LOGICAL NOUNIT, ONENRM, UPPER
81: CHARACTER NORMIN
82: INTEGER IX, KASE, KASE1
83: DOUBLE PRECISION AINVNM, ANORM, SCALE, SMLNUM, XNORM
84: * ..
85: * .. Local Arrays ..
86: INTEGER ISAVE( 3 )
87: * ..
88: * .. External Functions ..
89: LOGICAL LSAME
90: INTEGER IDAMAX
91: DOUBLE PRECISION DLAMCH, DLANTP
92: EXTERNAL LSAME, IDAMAX, DLAMCH, DLANTP
93: * ..
94: * .. External Subroutines ..
95: EXTERNAL DLACN2, DLATPS, DRSCL, XERBLA
96: * ..
97: * .. Intrinsic Functions ..
98: INTRINSIC ABS, DBLE, MAX
99: * ..
100: * .. Executable Statements ..
101: *
102: * Test the input parameters.
103: *
104: INFO = 0
105: UPPER = LSAME( UPLO, 'U' )
106: ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
107: NOUNIT = LSAME( DIAG, 'N' )
108: *
109: IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
110: INFO = -1
111: ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
112: INFO = -2
113: ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
114: INFO = -3
115: ELSE IF( N.LT.0 ) THEN
116: INFO = -4
117: END IF
118: IF( INFO.NE.0 ) THEN
119: CALL XERBLA( 'DTPCON', -INFO )
120: RETURN
121: END IF
122: *
123: * Quick return if possible
124: *
125: IF( N.EQ.0 ) THEN
126: RCOND = ONE
127: RETURN
128: END IF
129: *
130: RCOND = ZERO
131: SMLNUM = DLAMCH( 'Safe minimum' )*DBLE( MAX( 1, N ) )
132: *
133: * Compute the norm of the triangular matrix A.
134: *
135: ANORM = DLANTP( NORM, UPLO, DIAG, N, AP, WORK )
136: *
137: * Continue only if ANORM > 0.
138: *
139: IF( ANORM.GT.ZERO ) THEN
140: *
141: * Estimate the norm of the inverse of A.
142: *
143: AINVNM = ZERO
144: NORMIN = 'N'
145: IF( ONENRM ) THEN
146: KASE1 = 1
147: ELSE
148: KASE1 = 2
149: END IF
150: KASE = 0
151: 10 CONTINUE
152: CALL DLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE )
153: IF( KASE.NE.0 ) THEN
154: IF( KASE.EQ.KASE1 ) THEN
155: *
156: * Multiply by inv(A).
157: *
158: CALL DLATPS( UPLO, 'No transpose', DIAG, NORMIN, N, AP,
159: $ WORK, SCALE, WORK( 2*N+1 ), INFO )
160: ELSE
161: *
162: * Multiply by inv(A').
163: *
164: CALL DLATPS( UPLO, 'Transpose', DIAG, NORMIN, N, AP,
165: $ WORK, SCALE, WORK( 2*N+1 ), INFO )
166: END IF
167: NORMIN = 'Y'
168: *
169: * Multiply by 1/SCALE if doing so will not cause overflow.
170: *
171: IF( SCALE.NE.ONE ) THEN
172: IX = IDAMAX( N, WORK, 1 )
173: XNORM = ABS( WORK( IX ) )
174: IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
175: $ GO TO 20
176: CALL DRSCL( N, SCALE, WORK, 1 )
177: END IF
178: GO TO 10
179: END IF
180: *
181: * Compute the estimate of the reciprocal condition number.
182: *
183: IF( AINVNM.NE.ZERO )
184: $ RCOND = ( ONE / ANORM ) / AINVNM
185: END IF
186: *
187: 20 CONTINUE
188: RETURN
189: *
190: * End of DTPCON
191: *
192: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtpcon.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack