|
5164
|
1 SUBROUTINE DLATBS( UPLO, TRANS, DIAG, NORMIN, N, KD, AB, LDAB, X, |
|
|
2 $ SCALE, CNORM, INFO ) |
|
|
3 * |
|
|
4 * -- LAPACK auxiliary routine (version 3.0) -- |
|
|
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
6 * Courant Institute, Argonne National Lab, and Rice University |
|
|
7 * June 30, 1992 |
|
|
8 * |
|
|
9 * .. Scalar Arguments .. |
|
|
10 CHARACTER DIAG, NORMIN, TRANS, UPLO |
|
|
11 INTEGER INFO, KD, LDAB, N |
|
|
12 DOUBLE PRECISION SCALE |
|
|
13 * .. |
|
|
14 * .. Array Arguments .. |
|
|
15 DOUBLE PRECISION AB( LDAB, * ), CNORM( * ), X( * ) |
|
|
16 * .. |
|
|
17 * |
|
|
18 * Purpose |
|
|
19 * ======= |
|
|
20 * |
|
|
21 * DLATBS solves one of the triangular systems |
|
|
22 * |
|
|
23 * A *x = s*b or A'*x = s*b |
|
|
24 * |
|
|
25 * with scaling to prevent overflow, where A is an upper or lower |
|
|
26 * triangular band matrix. Here A' denotes the transpose of A, x and b |
|
|
27 * are n-element vectors, and s is a scaling factor, usually less than |
|
|
28 * or equal to 1, chosen so that the components of x will be less than |
|
|
29 * the overflow threshold. If the unscaled problem will not cause |
|
|
30 * overflow, the Level 2 BLAS routine DTBSV is called. If the matrix A |
|
|
31 * is singular (A(j,j) = 0 for some j), then s is set to 0 and a |
|
|
32 * non-trivial solution to A*x = 0 is returned. |
|
|
33 * |
|
|
34 * Arguments |
|
|
35 * ========= |
|
|
36 * |
|
|
37 * UPLO (input) CHARACTER*1 |
|
|
38 * Specifies whether the matrix A is upper or lower triangular. |
|
|
39 * = 'U': Upper triangular |
|
|
40 * = 'L': Lower triangular |
|
|
41 * |
|
|
42 * TRANS (input) CHARACTER*1 |
|
|
43 * Specifies the operation applied to A. |
|
|
44 * = 'N': Solve A * x = s*b (No transpose) |
|
|
45 * = 'T': Solve A'* x = s*b (Transpose) |
|
|
46 * = 'C': Solve A'* x = s*b (Conjugate transpose = Transpose) |
|
|
47 * |
|
|
48 * DIAG (input) CHARACTER*1 |
|
|
49 * Specifies whether or not the matrix A is unit triangular. |
|
|
50 * = 'N': Non-unit triangular |
|
|
51 * = 'U': Unit triangular |
|
|
52 * |
|
|
53 * NORMIN (input) CHARACTER*1 |
|
|
54 * Specifies whether CNORM has been set or not. |
|
|
55 * = 'Y': CNORM contains the column norms on entry |
|
|
56 * = 'N': CNORM is not set on entry. On exit, the norms will |
|
|
57 * be computed and stored in CNORM. |
|
|
58 * |
|
|
59 * N (input) INTEGER |
|
|
60 * The order of the matrix A. N >= 0. |
|
|
61 * |
|
|
62 * KD (input) INTEGER |
|
|
63 * The number of subdiagonals or superdiagonals in the |
|
|
64 * triangular matrix A. KD >= 0. |
|
|
65 * |
|
|
66 * AB (input) DOUBLE PRECISION array, dimension (LDAB,N) |
|
|
67 * The upper or lower triangular band matrix A, stored in the |
|
|
68 * first KD+1 rows of the array. The j-th column of A is stored |
|
|
69 * in the j-th column of the array AB as follows: |
|
|
70 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; |
|
|
71 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). |
|
|
72 * |
|
|
73 * LDAB (input) INTEGER |
|
|
74 * The leading dimension of the array AB. LDAB >= KD+1. |
|
|
75 * |
|
|
76 * X (input/output) DOUBLE PRECISION array, dimension (N) |
|
|
77 * On entry, the right hand side b of the triangular system. |
|
|
78 * On exit, X is overwritten by the solution vector x. |
|
|
79 * |
|
|
80 * SCALE (output) DOUBLE PRECISION |
|
|
81 * The scaling factor s for the triangular system |
|
|
82 * A * x = s*b or A'* x = s*b. |
|
|
83 * If SCALE = 0, the matrix A is singular or badly scaled, and |
|
|
84 * the vector x is an exact or approximate solution to A*x = 0. |
|
|
85 * |
|
|
86 * CNORM (input or output) DOUBLE PRECISION array, dimension (N) |
|
|
87 * |
|
|
88 * If NORMIN = 'Y', CNORM is an input argument and CNORM(j) |
|
|
89 * contains the norm of the off-diagonal part of the j-th column |
|
|
90 * of A. If TRANS = 'N', CNORM(j) must be greater than or equal |
|
|
91 * to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j) |
|
|
92 * must be greater than or equal to the 1-norm. |
|
|
93 * |
|
|
94 * If NORMIN = 'N', CNORM is an output argument and CNORM(j) |
|
|
95 * returns the 1-norm of the offdiagonal part of the j-th column |
|
|
96 * of A. |
|
|
97 * |
|
|
98 * INFO (output) INTEGER |
|
|
99 * = 0: successful exit |
|
|
100 * < 0: if INFO = -k, the k-th argument had an illegal value |
|
|
101 * |
|
|
102 * Further Details |
|
|
103 * ======= ======= |
|
|
104 * |
|
|
105 * A rough bound on x is computed; if that is less than overflow, DTBSV |
|
|
106 * is called, otherwise, specific code is used which checks for possible |
|
|
107 * overflow or divide-by-zero at every operation. |
|
|
108 * |
|
|
109 * A columnwise scheme is used for solving A*x = b. The basic algorithm |
|
|
110 * if A is lower triangular is |
|
|
111 * |
|
|
112 * x[1:n] := b[1:n] |
|
|
113 * for j = 1, ..., n |
|
|
114 * x(j) := x(j) / A(j,j) |
|
|
115 * x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j] |
|
|
116 * end |
|
|
117 * |
|
|
118 * Define bounds on the components of x after j iterations of the loop: |
|
|
119 * M(j) = bound on x[1:j] |
|
|
120 * G(j) = bound on x[j+1:n] |
|
|
121 * Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}. |
|
|
122 * |
|
|
123 * Then for iteration j+1 we have |
|
|
124 * M(j+1) <= G(j) / | A(j+1,j+1) | |
|
|
125 * G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] | |
|
|
126 * <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | ) |
|
|
127 * |
|
|
128 * where CNORM(j+1) is greater than or equal to the infinity-norm of |
|
|
129 * column j+1 of A, not counting the diagonal. Hence |
|
|
130 * |
|
|
131 * G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | ) |
|
|
132 * 1<=i<=j |
|
|
133 * and |
|
|
134 * |
|
|
135 * |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| ) |
|
|
136 * 1<=i< j |
|
|
137 * |
|
|
138 * Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTBSV if the |
|
|
139 * reciprocal of the largest M(j), j=1,..,n, is larger than |
|
|
140 * max(underflow, 1/overflow). |
|
|
141 * |
|
|
142 * The bound on x(j) is also used to determine when a step in the |
|
|
143 * columnwise method can be performed without fear of overflow. If |
|
|
144 * the computed bound is greater than a large constant, x is scaled to |
|
|
145 * prevent overflow, but if the bound overflows, x is set to 0, x(j) to |
|
|
146 * 1, and scale to 0, and a non-trivial solution to A*x = 0 is found. |
|
|
147 * |
|
|
148 * Similarly, a row-wise scheme is used to solve A'*x = b. The basic |
|
|
149 * algorithm for A upper triangular is |
|
|
150 * |
|
|
151 * for j = 1, ..., n |
|
|
152 * x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j) |
|
|
153 * end |
|
|
154 * |
|
|
155 * We simultaneously compute two bounds |
|
|
156 * G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j |
|
|
157 * M(j) = bound on x(i), 1<=i<=j |
|
|
158 * |
|
|
159 * The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we |
|
|
160 * add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1. |
|
|
161 * Then the bound on x(j) is |
|
|
162 * |
|
|
163 * M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) | |
|
|
164 * |
|
|
165 * <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| ) |
|
|
166 * 1<=i<=j |
|
|
167 * |
|
|
168 * and we can safely call DTBSV if 1/M(n) and 1/G(n) are both greater |
|
|
169 * than max(underflow, 1/overflow). |
|
|
170 * |
|
|
171 * ===================================================================== |
|
|
172 * |
|
|
173 * .. Parameters .. |
|
|
174 DOUBLE PRECISION ZERO, HALF, ONE |
|
|
175 PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) |
|
|
176 * .. |
|
|
177 * .. Local Scalars .. |
|
|
178 LOGICAL NOTRAN, NOUNIT, UPPER |
|
|
179 INTEGER I, IMAX, J, JFIRST, JINC, JLAST, JLEN, MAIND |
|
|
180 DOUBLE PRECISION BIGNUM, GROW, REC, SMLNUM, SUMJ, TJJ, TJJS, |
|
|
181 $ TMAX, TSCAL, USCAL, XBND, XJ, XMAX |
|
|
182 * .. |
|
|
183 * .. External Functions .. |
|
|
184 LOGICAL LSAME |
|
|
185 INTEGER IDAMAX |
|
|
186 DOUBLE PRECISION DASUM, DDOT, DLAMCH |
|
|
187 EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH |
|
|
188 * .. |
|
|
189 * .. External Subroutines .. |
|
|
190 EXTERNAL DAXPY, DSCAL, DTBSV, XERBLA |
|
|
191 * .. |
|
|
192 * .. Intrinsic Functions .. |
|
|
193 INTRINSIC ABS, MAX, MIN |
|
|
194 * .. |
|
|
195 * .. Executable Statements .. |
|
|
196 * |
|
|
197 INFO = 0 |
|
|
198 UPPER = LSAME( UPLO, 'U' ) |
|
|
199 NOTRAN = LSAME( TRANS, 'N' ) |
|
|
200 NOUNIT = LSAME( DIAG, 'N' ) |
|
|
201 * |
|
|
202 * Test the input parameters. |
|
|
203 * |
|
|
204 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN |
|
|
205 INFO = -1 |
|
|
206 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT. |
|
|
207 $ LSAME( TRANS, 'C' ) ) THEN |
|
|
208 INFO = -2 |
|
|
209 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN |
|
|
210 INFO = -3 |
|
|
211 ELSE IF( .NOT.LSAME( NORMIN, 'Y' ) .AND. .NOT. |
|
|
212 $ LSAME( NORMIN, 'N' ) ) THEN |
|
|
213 INFO = -4 |
|
|
214 ELSE IF( N.LT.0 ) THEN |
|
|
215 INFO = -5 |
|
|
216 ELSE IF( KD.LT.0 ) THEN |
|
|
217 INFO = -6 |
|
|
218 ELSE IF( LDAB.LT.KD+1 ) THEN |
|
|
219 INFO = -8 |
|
|
220 END IF |
|
|
221 IF( INFO.NE.0 ) THEN |
|
|
222 CALL XERBLA( 'DLATBS', -INFO ) |
|
|
223 RETURN |
|
|
224 END IF |
|
|
225 * |
|
|
226 * Quick return if possible |
|
|
227 * |
|
|
228 IF( N.EQ.0 ) |
|
|
229 $ RETURN |
|
|
230 * |
|
|
231 * Determine machine dependent parameters to control overflow. |
|
|
232 * |
|
|
233 SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' ) |
|
|
234 BIGNUM = ONE / SMLNUM |
|
|
235 SCALE = ONE |
|
|
236 * |
|
|
237 IF( LSAME( NORMIN, 'N' ) ) THEN |
|
|
238 * |
|
|
239 * Compute the 1-norm of each column, not including the diagonal. |
|
|
240 * |
|
|
241 IF( UPPER ) THEN |
|
|
242 * |
|
|
243 * A is upper triangular. |
|
|
244 * |
|
|
245 DO 10 J = 1, N |
|
|
246 JLEN = MIN( KD, J-1 ) |
|
|
247 CNORM( J ) = DASUM( JLEN, AB( KD+1-JLEN, J ), 1 ) |
|
|
248 10 CONTINUE |
|
|
249 ELSE |
|
|
250 * |
|
|
251 * A is lower triangular. |
|
|
252 * |
|
|
253 DO 20 J = 1, N |
|
|
254 JLEN = MIN( KD, N-J ) |
|
|
255 IF( JLEN.GT.0 ) THEN |
|
|
256 CNORM( J ) = DASUM( JLEN, AB( 2, J ), 1 ) |
|
|
257 ELSE |
|
|
258 CNORM( J ) = ZERO |
|
|
259 END IF |
|
|
260 20 CONTINUE |
|
|
261 END IF |
|
|
262 END IF |
|
|
263 * |
|
|
264 * Scale the column norms by TSCAL if the maximum element in CNORM is |
|
|
265 * greater than BIGNUM. |
|
|
266 * |
|
|
267 IMAX = IDAMAX( N, CNORM, 1 ) |
|
|
268 TMAX = CNORM( IMAX ) |
|
|
269 IF( TMAX.LE.BIGNUM ) THEN |
|
|
270 TSCAL = ONE |
|
|
271 ELSE |
|
|
272 TSCAL = ONE / ( SMLNUM*TMAX ) |
|
|
273 CALL DSCAL( N, TSCAL, CNORM, 1 ) |
|
|
274 END IF |
|
|
275 * |
|
|
276 * Compute a bound on the computed solution vector to see if the |
|
|
277 * Level 2 BLAS routine DTBSV can be used. |
|
|
278 * |
|
|
279 J = IDAMAX( N, X, 1 ) |
|
|
280 XMAX = ABS( X( J ) ) |
|
|
281 XBND = XMAX |
|
|
282 IF( NOTRAN ) THEN |
|
|
283 * |
|
|
284 * Compute the growth in A * x = b. |
|
|
285 * |
|
|
286 IF( UPPER ) THEN |
|
|
287 JFIRST = N |
|
|
288 JLAST = 1 |
|
|
289 JINC = -1 |
|
|
290 MAIND = KD + 1 |
|
|
291 ELSE |
|
|
292 JFIRST = 1 |
|
|
293 JLAST = N |
|
|
294 JINC = 1 |
|
|
295 MAIND = 1 |
|
|
296 END IF |
|
|
297 * |
|
|
298 IF( TSCAL.NE.ONE ) THEN |
|
|
299 GROW = ZERO |
|
|
300 GO TO 50 |
|
|
301 END IF |
|
|
302 * |
|
|
303 IF( NOUNIT ) THEN |
|
|
304 * |
|
|
305 * A is non-unit triangular. |
|
|
306 * |
|
|
307 * Compute GROW = 1/G(j) and XBND = 1/M(j). |
|
|
308 * Initially, G(0) = max{x(i), i=1,...,n}. |
|
|
309 * |
|
|
310 GROW = ONE / MAX( XBND, SMLNUM ) |
|
|
311 XBND = GROW |
|
|
312 DO 30 J = JFIRST, JLAST, JINC |
|
|
313 * |
|
|
314 * Exit the loop if the growth factor is too small. |
|
|
315 * |
|
|
316 IF( GROW.LE.SMLNUM ) |
|
|
317 $ GO TO 50 |
|
|
318 * |
|
|
319 * M(j) = G(j-1) / abs(A(j,j)) |
|
|
320 * |
|
|
321 TJJ = ABS( AB( MAIND, J ) ) |
|
|
322 XBND = MIN( XBND, MIN( ONE, TJJ )*GROW ) |
|
|
323 IF( TJJ+CNORM( J ).GE.SMLNUM ) THEN |
|
|
324 * |
|
|
325 * G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) ) |
|
|
326 * |
|
|
327 GROW = GROW*( TJJ / ( TJJ+CNORM( J ) ) ) |
|
|
328 ELSE |
|
|
329 * |
|
|
330 * G(j) could overflow, set GROW to 0. |
|
|
331 * |
|
|
332 GROW = ZERO |
|
|
333 END IF |
|
|
334 30 CONTINUE |
|
|
335 GROW = XBND |
|
|
336 ELSE |
|
|
337 * |
|
|
338 * A is unit triangular. |
|
|
339 * |
|
|
340 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. |
|
|
341 * |
|
|
342 GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) ) |
|
|
343 DO 40 J = JFIRST, JLAST, JINC |
|
|
344 * |
|
|
345 * Exit the loop if the growth factor is too small. |
|
|
346 * |
|
|
347 IF( GROW.LE.SMLNUM ) |
|
|
348 $ GO TO 50 |
|
|
349 * |
|
|
350 * G(j) = G(j-1)*( 1 + CNORM(j) ) |
|
|
351 * |
|
|
352 GROW = GROW*( ONE / ( ONE+CNORM( J ) ) ) |
|
|
353 40 CONTINUE |
|
|
354 END IF |
|
|
355 50 CONTINUE |
|
|
356 * |
|
|
357 ELSE |
|
|
358 * |
|
|
359 * Compute the growth in A' * x = b. |
|
|
360 * |
|
|
361 IF( UPPER ) THEN |
|
|
362 JFIRST = 1 |
|
|
363 JLAST = N |
|
|
364 JINC = 1 |
|
|
365 MAIND = KD + 1 |
|
|
366 ELSE |
|
|
367 JFIRST = N |
|
|
368 JLAST = 1 |
|
|
369 JINC = -1 |
|
|
370 MAIND = 1 |
|
|
371 END IF |
|
|
372 * |
|
|
373 IF( TSCAL.NE.ONE ) THEN |
|
|
374 GROW = ZERO |
|
|
375 GO TO 80 |
|
|
376 END IF |
|
|
377 * |
|
|
378 IF( NOUNIT ) THEN |
|
|
379 * |
|
|
380 * A is non-unit triangular. |
|
|
381 * |
|
|
382 * Compute GROW = 1/G(j) and XBND = 1/M(j). |
|
|
383 * Initially, M(0) = max{x(i), i=1,...,n}. |
|
|
384 * |
|
|
385 GROW = ONE / MAX( XBND, SMLNUM ) |
|
|
386 XBND = GROW |
|
|
387 DO 60 J = JFIRST, JLAST, JINC |
|
|
388 * |
|
|
389 * Exit the loop if the growth factor is too small. |
|
|
390 * |
|
|
391 IF( GROW.LE.SMLNUM ) |
|
|
392 $ GO TO 80 |
|
|
393 * |
|
|
394 * G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) ) |
|
|
395 * |
|
|
396 XJ = ONE + CNORM( J ) |
|
|
397 GROW = MIN( GROW, XBND / XJ ) |
|
|
398 * |
|
|
399 * M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j)) |
|
|
400 * |
|
|
401 TJJ = ABS( AB( MAIND, J ) ) |
|
|
402 IF( XJ.GT.TJJ ) |
|
|
403 $ XBND = XBND*( TJJ / XJ ) |
|
|
404 60 CONTINUE |
|
|
405 GROW = MIN( GROW, XBND ) |
|
|
406 ELSE |
|
|
407 * |
|
|
408 * A is unit triangular. |
|
|
409 * |
|
|
410 * Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}. |
|
|
411 * |
|
|
412 GROW = MIN( ONE, ONE / MAX( XBND, SMLNUM ) ) |
|
|
413 DO 70 J = JFIRST, JLAST, JINC |
|
|
414 * |
|
|
415 * Exit the loop if the growth factor is too small. |
|
|
416 * |
|
|
417 IF( GROW.LE.SMLNUM ) |
|
|
418 $ GO TO 80 |
|
|
419 * |
|
|
420 * G(j) = ( 1 + CNORM(j) )*G(j-1) |
|
|
421 * |
|
|
422 XJ = ONE + CNORM( J ) |
|
|
423 GROW = GROW / XJ |
|
|
424 70 CONTINUE |
|
|
425 END IF |
|
|
426 80 CONTINUE |
|
|
427 END IF |
|
|
428 * |
|
|
429 IF( ( GROW*TSCAL ).GT.SMLNUM ) THEN |
|
|
430 * |
|
|
431 * Use the Level 2 BLAS solve if the reciprocal of the bound on |
|
|
432 * elements of X is not too small. |
|
|
433 * |
|
|
434 CALL DTBSV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, X, 1 ) |
|
|
435 ELSE |
|
|
436 * |
|
|
437 * Use a Level 1 BLAS solve, scaling intermediate results. |
|
|
438 * |
|
|
439 IF( XMAX.GT.BIGNUM ) THEN |
|
|
440 * |
|
|
441 * Scale X so that its components are less than or equal to |
|
|
442 * BIGNUM in absolute value. |
|
|
443 * |
|
|
444 SCALE = BIGNUM / XMAX |
|
|
445 CALL DSCAL( N, SCALE, X, 1 ) |
|
|
446 XMAX = BIGNUM |
|
|
447 END IF |
|
|
448 * |
|
|
449 IF( NOTRAN ) THEN |
|
|
450 * |
|
|
451 * Solve A * x = b |
|
|
452 * |
|
|
453 DO 110 J = JFIRST, JLAST, JINC |
|
|
454 * |
|
|
455 * Compute x(j) = b(j) / A(j,j), scaling x if necessary. |
|
|
456 * |
|
|
457 XJ = ABS( X( J ) ) |
|
|
458 IF( NOUNIT ) THEN |
|
|
459 TJJS = AB( MAIND, J )*TSCAL |
|
|
460 ELSE |
|
|
461 TJJS = TSCAL |
|
|
462 IF( TSCAL.EQ.ONE ) |
|
|
463 $ GO TO 100 |
|
|
464 END IF |
|
|
465 TJJ = ABS( TJJS ) |
|
|
466 IF( TJJ.GT.SMLNUM ) THEN |
|
|
467 * |
|
|
468 * abs(A(j,j)) > SMLNUM: |
|
|
469 * |
|
|
470 IF( TJJ.LT.ONE ) THEN |
|
|
471 IF( XJ.GT.TJJ*BIGNUM ) THEN |
|
|
472 * |
|
|
473 * Scale x by 1/b(j). |
|
|
474 * |
|
|
475 REC = ONE / XJ |
|
|
476 CALL DSCAL( N, REC, X, 1 ) |
|
|
477 SCALE = SCALE*REC |
|
|
478 XMAX = XMAX*REC |
|
|
479 END IF |
|
|
480 END IF |
|
|
481 X( J ) = X( J ) / TJJS |
|
|
482 XJ = ABS( X( J ) ) |
|
|
483 ELSE IF( TJJ.GT.ZERO ) THEN |
|
|
484 * |
|
|
485 * 0 < abs(A(j,j)) <= SMLNUM: |
|
|
486 * |
|
|
487 IF( XJ.GT.TJJ*BIGNUM ) THEN |
|
|
488 * |
|
|
489 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM |
|
|
490 * to avoid overflow when dividing by A(j,j). |
|
|
491 * |
|
|
492 REC = ( TJJ*BIGNUM ) / XJ |
|
|
493 IF( CNORM( J ).GT.ONE ) THEN |
|
|
494 * |
|
|
495 * Scale by 1/CNORM(j) to avoid overflow when |
|
|
496 * multiplying x(j) times column j. |
|
|
497 * |
|
|
498 REC = REC / CNORM( J ) |
|
|
499 END IF |
|
|
500 CALL DSCAL( N, REC, X, 1 ) |
|
|
501 SCALE = SCALE*REC |
|
|
502 XMAX = XMAX*REC |
|
|
503 END IF |
|
|
504 X( J ) = X( J ) / TJJS |
|
|
505 XJ = ABS( X( J ) ) |
|
|
506 ELSE |
|
|
507 * |
|
|
508 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and |
|
|
509 * scale = 0, and compute a solution to A*x = 0. |
|
|
510 * |
|
|
511 DO 90 I = 1, N |
|
|
512 X( I ) = ZERO |
|
|
513 90 CONTINUE |
|
|
514 X( J ) = ONE |
|
|
515 XJ = ONE |
|
|
516 SCALE = ZERO |
|
|
517 XMAX = ZERO |
|
|
518 END IF |
|
|
519 100 CONTINUE |
|
|
520 * |
|
|
521 * Scale x if necessary to avoid overflow when adding a |
|
|
522 * multiple of column j of A. |
|
|
523 * |
|
|
524 IF( XJ.GT.ONE ) THEN |
|
|
525 REC = ONE / XJ |
|
|
526 IF( CNORM( J ).GT.( BIGNUM-XMAX )*REC ) THEN |
|
|
527 * |
|
|
528 * Scale x by 1/(2*abs(x(j))). |
|
|
529 * |
|
|
530 REC = REC*HALF |
|
|
531 CALL DSCAL( N, REC, X, 1 ) |
|
|
532 SCALE = SCALE*REC |
|
|
533 END IF |
|
|
534 ELSE IF( XJ*CNORM( J ).GT.( BIGNUM-XMAX ) ) THEN |
|
|
535 * |
|
|
536 * Scale x by 1/2. |
|
|
537 * |
|
|
538 CALL DSCAL( N, HALF, X, 1 ) |
|
|
539 SCALE = SCALE*HALF |
|
|
540 END IF |
|
|
541 * |
|
|
542 IF( UPPER ) THEN |
|
|
543 IF( J.GT.1 ) THEN |
|
|
544 * |
|
|
545 * Compute the update |
|
|
546 * x(max(1,j-kd):j-1) := x(max(1,j-kd):j-1) - |
|
|
547 * x(j)* A(max(1,j-kd):j-1,j) |
|
|
548 * |
|
|
549 JLEN = MIN( KD, J-1 ) |
|
|
550 CALL DAXPY( JLEN, -X( J )*TSCAL, |
|
|
551 $ AB( KD+1-JLEN, J ), 1, X( J-JLEN ), 1 ) |
|
|
552 I = IDAMAX( J-1, X, 1 ) |
|
|
553 XMAX = ABS( X( I ) ) |
|
|
554 END IF |
|
|
555 ELSE IF( J.LT.N ) THEN |
|
|
556 * |
|
|
557 * Compute the update |
|
|
558 * x(j+1:min(j+kd,n)) := x(j+1:min(j+kd,n)) - |
|
|
559 * x(j) * A(j+1:min(j+kd,n),j) |
|
|
560 * |
|
|
561 JLEN = MIN( KD, N-J ) |
|
|
562 IF( JLEN.GT.0 ) |
|
|
563 $ CALL DAXPY( JLEN, -X( J )*TSCAL, AB( 2, J ), 1, |
|
|
564 $ X( J+1 ), 1 ) |
|
|
565 I = J + IDAMAX( N-J, X( J+1 ), 1 ) |
|
|
566 XMAX = ABS( X( I ) ) |
|
|
567 END IF |
|
|
568 110 CONTINUE |
|
|
569 * |
|
|
570 ELSE |
|
|
571 * |
|
|
572 * Solve A' * x = b |
|
|
573 * |
|
|
574 DO 160 J = JFIRST, JLAST, JINC |
|
|
575 * |
|
|
576 * Compute x(j) = b(j) - sum A(k,j)*x(k). |
|
|
577 * k<>j |
|
|
578 * |
|
|
579 XJ = ABS( X( J ) ) |
|
|
580 USCAL = TSCAL |
|
|
581 REC = ONE / MAX( XMAX, ONE ) |
|
|
582 IF( CNORM( J ).GT.( BIGNUM-XJ )*REC ) THEN |
|
|
583 * |
|
|
584 * If x(j) could overflow, scale x by 1/(2*XMAX). |
|
|
585 * |
|
|
586 REC = REC*HALF |
|
|
587 IF( NOUNIT ) THEN |
|
|
588 TJJS = AB( MAIND, J )*TSCAL |
|
|
589 ELSE |
|
|
590 TJJS = TSCAL |
|
|
591 END IF |
|
|
592 TJJ = ABS( TJJS ) |
|
|
593 IF( TJJ.GT.ONE ) THEN |
|
|
594 * |
|
|
595 * Divide by A(j,j) when scaling x if A(j,j) > 1. |
|
|
596 * |
|
|
597 REC = MIN( ONE, REC*TJJ ) |
|
|
598 USCAL = USCAL / TJJS |
|
|
599 END IF |
|
|
600 IF( REC.LT.ONE ) THEN |
|
|
601 CALL DSCAL( N, REC, X, 1 ) |
|
|
602 SCALE = SCALE*REC |
|
|
603 XMAX = XMAX*REC |
|
|
604 END IF |
|
|
605 END IF |
|
|
606 * |
|
|
607 SUMJ = ZERO |
|
|
608 IF( USCAL.EQ.ONE ) THEN |
|
|
609 * |
|
|
610 * If the scaling needed for A in the dot product is 1, |
|
|
611 * call DDOT to perform the dot product. |
|
|
612 * |
|
|
613 IF( UPPER ) THEN |
|
|
614 JLEN = MIN( KD, J-1 ) |
|
|
615 SUMJ = DDOT( JLEN, AB( KD+1-JLEN, J ), 1, |
|
|
616 $ X( J-JLEN ), 1 ) |
|
|
617 ELSE |
|
|
618 JLEN = MIN( KD, N-J ) |
|
|
619 IF( JLEN.GT.0 ) |
|
|
620 $ SUMJ = DDOT( JLEN, AB( 2, J ), 1, X( J+1 ), 1 ) |
|
|
621 END IF |
|
|
622 ELSE |
|
|
623 * |
|
|
624 * Otherwise, use in-line code for the dot product. |
|
|
625 * |
|
|
626 IF( UPPER ) THEN |
|
|
627 JLEN = MIN( KD, J-1 ) |
|
|
628 DO 120 I = 1, JLEN |
|
|
629 SUMJ = SUMJ + ( AB( KD+I-JLEN, J )*USCAL )* |
|
|
630 $ X( J-JLEN-1+I ) |
|
|
631 120 CONTINUE |
|
|
632 ELSE |
|
|
633 JLEN = MIN( KD, N-J ) |
|
|
634 DO 130 I = 1, JLEN |
|
|
635 SUMJ = SUMJ + ( AB( I+1, J )*USCAL )*X( J+I ) |
|
|
636 130 CONTINUE |
|
|
637 END IF |
|
|
638 END IF |
|
|
639 * |
|
|
640 IF( USCAL.EQ.TSCAL ) THEN |
|
|
641 * |
|
|
642 * Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j) |
|
|
643 * was not used to scale the dotproduct. |
|
|
644 * |
|
|
645 X( J ) = X( J ) - SUMJ |
|
|
646 XJ = ABS( X( J ) ) |
|
|
647 IF( NOUNIT ) THEN |
|
|
648 * |
|
|
649 * Compute x(j) = x(j) / A(j,j), scaling if necessary. |
|
|
650 * |
|
|
651 TJJS = AB( MAIND, J )*TSCAL |
|
|
652 ELSE |
|
|
653 TJJS = TSCAL |
|
|
654 IF( TSCAL.EQ.ONE ) |
|
|
655 $ GO TO 150 |
|
|
656 END IF |
|
|
657 TJJ = ABS( TJJS ) |
|
|
658 IF( TJJ.GT.SMLNUM ) THEN |
|
|
659 * |
|
|
660 * abs(A(j,j)) > SMLNUM: |
|
|
661 * |
|
|
662 IF( TJJ.LT.ONE ) THEN |
|
|
663 IF( XJ.GT.TJJ*BIGNUM ) THEN |
|
|
664 * |
|
|
665 * Scale X by 1/abs(x(j)). |
|
|
666 * |
|
|
667 REC = ONE / XJ |
|
|
668 CALL DSCAL( N, REC, X, 1 ) |
|
|
669 SCALE = SCALE*REC |
|
|
670 XMAX = XMAX*REC |
|
|
671 END IF |
|
|
672 END IF |
|
|
673 X( J ) = X( J ) / TJJS |
|
|
674 ELSE IF( TJJ.GT.ZERO ) THEN |
|
|
675 * |
|
|
676 * 0 < abs(A(j,j)) <= SMLNUM: |
|
|
677 * |
|
|
678 IF( XJ.GT.TJJ*BIGNUM ) THEN |
|
|
679 * |
|
|
680 * Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM. |
|
|
681 * |
|
|
682 REC = ( TJJ*BIGNUM ) / XJ |
|
|
683 CALL DSCAL( N, REC, X, 1 ) |
|
|
684 SCALE = SCALE*REC |
|
|
685 XMAX = XMAX*REC |
|
|
686 END IF |
|
|
687 X( J ) = X( J ) / TJJS |
|
|
688 ELSE |
|
|
689 * |
|
|
690 * A(j,j) = 0: Set x(1:n) = 0, x(j) = 1, and |
|
|
691 * scale = 0, and compute a solution to A'*x = 0. |
|
|
692 * |
|
|
693 DO 140 I = 1, N |
|
|
694 X( I ) = ZERO |
|
|
695 140 CONTINUE |
|
|
696 X( J ) = ONE |
|
|
697 SCALE = ZERO |
|
|
698 XMAX = ZERO |
|
|
699 END IF |
|
|
700 150 CONTINUE |
|
|
701 ELSE |
|
|
702 * |
|
|
703 * Compute x(j) := x(j) / A(j,j) - sumj if the dot |
|
|
704 * product has already been divided by 1/A(j,j). |
|
|
705 * |
|
|
706 X( J ) = X( J ) / TJJS - SUMJ |
|
|
707 END IF |
|
|
708 XMAX = MAX( XMAX, ABS( X( J ) ) ) |
|
|
709 160 CONTINUE |
|
|
710 END IF |
|
|
711 SCALE = SCALE / TSCAL |
|
|
712 END IF |
|
|
713 * |
|
|
714 * Scale the column norms by 1/TSCAL for return. |
|
|
715 * |
|
|
716 IF( TSCAL.NE.ONE ) THEN |
|
|
717 CALL DSCAL( N, ONE / TSCAL, CNORM, 1 ) |
|
|
718 END IF |
|
|
719 * |
|
|
720 RETURN |
|
|
721 * |
|
|
722 * End of DLATBS |
|
|
723 * |
|
|
724 END |
