Mercurial > octave
comparison libcruft/lapack/spbcon.f @ 7789:82be108cc558
First attempt at single precision tyeps
* * *
corrections to qrupdate single precision routines
* * *
prefer demotion to single over promotion to double
* * *
Add single precision support to log2 function
* * *
Trivial PROJECT file update
* * *
Cache optimized hermitian/transpose methods
* * *
Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author | David Bateman <dbateman@free.fr> |
---|---|
date | Sun, 27 Apr 2008 22:34:17 +0200 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
7788:45f5faba05a2 | 7789:82be108cc558 |
---|---|
1 SUBROUTINE SPBCON( UPLO, N, KD, AB, LDAB, ANORM, RCOND, WORK, | |
2 $ IWORK, INFO ) | |
3 * | |
4 * -- LAPACK routine (version 3.1) -- | |
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. | |
6 * November 2006 | |
7 * | |
8 * Modified to call SLACN2 in place of SLACON, 7 Feb 03, SJH. | |
9 * | |
10 * .. Scalar Arguments .. | |
11 CHARACTER UPLO | |
12 INTEGER INFO, KD, LDAB, N | |
13 REAL ANORM, RCOND | |
14 * .. | |
15 * .. Array Arguments .. | |
16 INTEGER IWORK( * ) | |
17 REAL AB( LDAB, * ), WORK( * ) | |
18 * .. | |
19 * | |
20 * Purpose | |
21 * ======= | |
22 * | |
23 * SPBCON estimates the reciprocal of the condition number (in the | |
24 * 1-norm) of a real symmetric positive definite band matrix using the | |
25 * Cholesky factorization A = U**T*U or A = L*L**T computed by SPBTRF. | |
26 * | |
27 * An estimate is obtained for norm(inv(A)), and the reciprocal of the | |
28 * condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))). | |
29 * | |
30 * Arguments | |
31 * ========= | |
32 * | |
33 * UPLO (input) CHARACTER*1 | |
34 * = 'U': Upper triangular factor stored in AB; | |
35 * = 'L': Lower triangular factor stored in AB. | |
36 * | |
37 * N (input) INTEGER | |
38 * The order of the matrix A. N >= 0. | |
39 * | |
40 * KD (input) INTEGER | |
41 * The number of superdiagonals of the matrix A if UPLO = 'U', | |
42 * or the number of subdiagonals if UPLO = 'L'. KD >= 0. | |
43 * | |
44 * AB (input) REAL array, dimension (LDAB,N) | |
45 * The triangular factor U or L from the Cholesky factorization | |
46 * A = U**T*U or A = L*L**T of the band matrix A, stored in the | |
47 * first KD+1 rows of the array. The j-th column of U or L is | |
48 * stored in the j-th column of the array AB as follows: | |
49 * if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j; | |
50 * if UPLO ='L', AB(1+i-j,j) = L(i,j) for j<=i<=min(n,j+kd). | |
51 * | |
52 * LDAB (input) INTEGER | |
53 * The leading dimension of the array AB. LDAB >= KD+1. | |
54 * | |
55 * ANORM (input) REAL | |
56 * The 1-norm (or infinity-norm) of the symmetric band matrix A. | |
57 * | |
58 * RCOND (output) REAL | |
59 * The reciprocal of the condition number of the matrix A, | |
60 * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an | |
61 * estimate of the 1-norm of inv(A) computed in this routine. | |
62 * | |
63 * WORK (workspace) REAL array, dimension (3*N) | |
64 * | |
65 * IWORK (workspace) INTEGER array, dimension (N) | |
66 * | |
67 * INFO (output) INTEGER | |
68 * = 0: successful exit | |
69 * < 0: if INFO = -i, the i-th argument had an illegal value | |
70 * | |
71 * ===================================================================== | |
72 * | |
73 * .. Parameters .. | |
74 REAL ONE, ZERO | |
75 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) | |
76 * .. | |
77 * .. Local Scalars .. | |
78 LOGICAL UPPER | |
79 CHARACTER NORMIN | |
80 INTEGER IX, KASE | |
81 REAL AINVNM, SCALE, SCALEL, SCALEU, SMLNUM | |
82 * .. | |
83 * .. Local Arrays .. | |
84 INTEGER ISAVE( 3 ) | |
85 * .. | |
86 * .. External Functions .. | |
87 LOGICAL LSAME | |
88 INTEGER ISAMAX | |
89 REAL SLAMCH | |
90 EXTERNAL LSAME, ISAMAX, SLAMCH | |
91 * .. | |
92 * .. External Subroutines .. | |
93 EXTERNAL SLACN2, SLATBS, SRSCL, XERBLA | |
94 * .. | |
95 * .. Intrinsic Functions .. | |
96 INTRINSIC ABS | |
97 * .. | |
98 * .. Executable Statements .. | |
99 * | |
100 * Test the input parameters. | |
101 * | |
102 INFO = 0 | |
103 UPPER = LSAME( UPLO, 'U' ) | |
104 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN | |
105 INFO = -1 | |
106 ELSE IF( N.LT.0 ) THEN | |
107 INFO = -2 | |
108 ELSE IF( KD.LT.0 ) THEN | |
109 INFO = -3 | |
110 ELSE IF( LDAB.LT.KD+1 ) THEN | |
111 INFO = -5 | |
112 ELSE IF( ANORM.LT.ZERO ) THEN | |
113 INFO = -6 | |
114 END IF | |
115 IF( INFO.NE.0 ) THEN | |
116 CALL XERBLA( 'SPBCON', -INFO ) | |
117 RETURN | |
118 END IF | |
119 * | |
120 * Quick return if possible | |
121 * | |
122 RCOND = ZERO | |
123 IF( N.EQ.0 ) THEN | |
124 RCOND = ONE | |
125 RETURN | |
126 ELSE IF( ANORM.EQ.ZERO ) THEN | |
127 RETURN | |
128 END IF | |
129 * | |
130 SMLNUM = SLAMCH( 'Safe minimum' ) | |
131 * | |
132 * Estimate the 1-norm of the inverse. | |
133 * | |
134 KASE = 0 | |
135 NORMIN = 'N' | |
136 10 CONTINUE | |
137 CALL SLACN2( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE, ISAVE ) | |
138 IF( KASE.NE.0 ) THEN | |
139 IF( UPPER ) THEN | |
140 * | |
141 * Multiply by inv(U'). | |
142 * | |
143 CALL SLATBS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, | |
144 $ KD, AB, LDAB, WORK, SCALEL, WORK( 2*N+1 ), | |
145 $ INFO ) | |
146 NORMIN = 'Y' | |
147 * | |
148 * Multiply by inv(U). | |
149 * | |
150 CALL SLATBS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N, | |
151 $ KD, AB, LDAB, WORK, SCALEU, WORK( 2*N+1 ), | |
152 $ INFO ) | |
153 ELSE | |
154 * | |
155 * Multiply by inv(L). | |
156 * | |
157 CALL SLATBS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N, | |
158 $ KD, AB, LDAB, WORK, SCALEL, WORK( 2*N+1 ), | |
159 $ INFO ) | |
160 NORMIN = 'Y' | |
161 * | |
162 * Multiply by inv(L'). | |
163 * | |
164 CALL SLATBS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N, | |
165 $ KD, AB, LDAB, WORK, SCALEU, WORK( 2*N+1 ), | |
166 $ INFO ) | |
167 END IF | |
168 * | |
169 * Multiply by 1/SCALE if doing so will not cause overflow. | |
170 * | |
171 SCALE = SCALEL*SCALEU | |
172 IF( SCALE.NE.ONE ) THEN | |
173 IX = ISAMAX( N, WORK, 1 ) | |
174 IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO ) | |
175 $ GO TO 20 | |
176 CALL SRSCL( 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 / AINVNM ) / ANORM | |
185 * | |
186 20 CONTINUE | |
187 * | |
188 RETURN | |
189 * | |
190 * End of SPBCON | |
191 * | |
192 END |