|
2329
|
1 SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, |
|
|
2 $ LDB, WR, WI, X, LDX, SCALE, XNORM, INFO ) |
|
|
3 * |
|
7034
|
4 * -- LAPACK auxiliary routine (version 3.1) -- |
|
|
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
6 * November 2006 |
|
2329
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
|
9 LOGICAL LTRANS |
|
|
10 INTEGER INFO, LDA, LDB, LDX, NA, NW |
|
|
11 DOUBLE PRECISION CA, D1, D2, SCALE, SMIN, WI, WR, XNORM |
|
|
12 * .. |
|
|
13 * .. Array Arguments .. |
|
|
14 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * ) |
|
|
15 * .. |
|
|
16 * |
|
|
17 * Purpose |
|
|
18 * ======= |
|
|
19 * |
|
|
20 * DLALN2 solves a system of the form (ca A - w D ) X = s B |
|
|
21 * or (ca A' - w D) X = s B with possible scaling ("s") and |
|
|
22 * perturbation of A. (A' means A-transpose.) |
|
|
23 * |
|
|
24 * A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA |
|
|
25 * real diagonal matrix, w is a real or complex value, and X and B are |
|
|
26 * NA x 1 matrices -- real if w is real, complex if w is complex. NA |
|
|
27 * may be 1 or 2. |
|
|
28 * |
|
|
29 * If w is complex, X and B are represented as NA x 2 matrices, |
|
|
30 * the first column of each being the real part and the second |
|
|
31 * being the imaginary part. |
|
|
32 * |
|
|
33 * "s" is a scaling factor (.LE. 1), computed by DLALN2, which is |
|
|
34 * so chosen that X can be computed without overflow. X is further |
|
|
35 * scaled if necessary to assure that norm(ca A - w D)*norm(X) is less |
|
|
36 * than overflow. |
|
|
37 * |
|
|
38 * If both singular values of (ca A - w D) are less than SMIN, |
|
|
39 * SMIN*identity will be used instead of (ca A - w D). If only one |
|
|
40 * singular value is less than SMIN, one element of (ca A - w D) will be |
|
|
41 * perturbed enough to make the smallest singular value roughly SMIN. |
|
|
42 * If both singular values are at least SMIN, (ca A - w D) will not be |
|
|
43 * perturbed. In any case, the perturbation will be at most some small |
|
|
44 * multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values |
|
|
45 * are computed by infinity-norm approximations, and thus will only be |
|
|
46 * correct to a factor of 2 or so. |
|
|
47 * |
|
|
48 * Note: all input quantities are assumed to be smaller than overflow |
|
|
49 * by a reasonable factor. (See BIGNUM.) |
|
|
50 * |
|
|
51 * Arguments |
|
|
52 * ========== |
|
|
53 * |
|
|
54 * LTRANS (input) LOGICAL |
|
|
55 * =.TRUE.: A-transpose will be used. |
|
|
56 * =.FALSE.: A will be used (not transposed.) |
|
|
57 * |
|
|
58 * NA (input) INTEGER |
|
|
59 * The size of the matrix A. It may (only) be 1 or 2. |
|
|
60 * |
|
|
61 * NW (input) INTEGER |
|
|
62 * 1 if "w" is real, 2 if "w" is complex. It may only be 1 |
|
|
63 * or 2. |
|
|
64 * |
|
|
65 * SMIN (input) DOUBLE PRECISION |
|
|
66 * The desired lower bound on the singular values of A. This |
|
|
67 * should be a safe distance away from underflow or overflow, |
|
|
68 * say, between (underflow/machine precision) and (machine |
|
|
69 * precision * overflow ). (See BIGNUM and ULP.) |
|
|
70 * |
|
|
71 * CA (input) DOUBLE PRECISION |
|
|
72 * The coefficient c, which A is multiplied by. |
|
|
73 * |
|
|
74 * A (input) DOUBLE PRECISION array, dimension (LDA,NA) |
|
|
75 * The NA x NA matrix A. |
|
|
76 * |
|
|
77 * LDA (input) INTEGER |
|
|
78 * The leading dimension of A. It must be at least NA. |
|
|
79 * |
|
|
80 * D1 (input) DOUBLE PRECISION |
|
|
81 * The 1,1 element in the diagonal matrix D. |
|
|
82 * |
|
|
83 * D2 (input) DOUBLE PRECISION |
|
|
84 * The 2,2 element in the diagonal matrix D. Not used if NW=1. |
|
|
85 * |
|
|
86 * B (input) DOUBLE PRECISION array, dimension (LDB,NW) |
|
|
87 * The NA x NW matrix B (right-hand side). If NW=2 ("w" is |
|
|
88 * complex), column 1 contains the real part of B and column 2 |
|
|
89 * contains the imaginary part. |
|
|
90 * |
|
|
91 * LDB (input) INTEGER |
|
|
92 * The leading dimension of B. It must be at least NA. |
|
|
93 * |
|
|
94 * WR (input) DOUBLE PRECISION |
|
|
95 * The real part of the scalar "w". |
|
|
96 * |
|
|
97 * WI (input) DOUBLE PRECISION |
|
|
98 * The imaginary part of the scalar "w". Not used if NW=1. |
|
|
99 * |
|
|
100 * X (output) DOUBLE PRECISION array, dimension (LDX,NW) |
|
|
101 * The NA x NW matrix X (unknowns), as computed by DLALN2. |
|
|
102 * If NW=2 ("w" is complex), on exit, column 1 will contain |
|
|
103 * the real part of X and column 2 will contain the imaginary |
|
|
104 * part. |
|
|
105 * |
|
|
106 * LDX (input) INTEGER |
|
|
107 * The leading dimension of X. It must be at least NA. |
|
|
108 * |
|
|
109 * SCALE (output) DOUBLE PRECISION |
|
|
110 * The scale factor that B must be multiplied by to insure |
|
|
111 * that overflow does not occur when computing X. Thus, |
|
|
112 * (ca A - w D) X will be SCALE*B, not B (ignoring |
|
|
113 * perturbations of A.) It will be at most 1. |
|
|
114 * |
|
|
115 * XNORM (output) DOUBLE PRECISION |
|
|
116 * The infinity-norm of X, when X is regarded as an NA x NW |
|
|
117 * real matrix. |
|
|
118 * |
|
|
119 * INFO (output) INTEGER |
|
|
120 * An error flag. It will be set to zero if no error occurs, |
|
|
121 * a negative number if an argument is in error, or a positive |
|
|
122 * number if ca A - w D had to be perturbed. |
|
|
123 * The possible values are: |
|
|
124 * = 0: No error occurred, and (ca A - w D) did not have to be |
|
|
125 * perturbed. |
|
|
126 * = 1: (ca A - w D) had to be perturbed to make its smallest |
|
|
127 * (or only) singular value greater than SMIN. |
|
|
128 * NOTE: In the interests of speed, this routine does not |
|
|
129 * check the inputs for errors. |
|
|
130 * |
|
|
131 * ===================================================================== |
|
|
132 * |
|
|
133 * .. Parameters .. |
|
|
134 DOUBLE PRECISION ZERO, ONE |
|
|
135 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) |
|
|
136 DOUBLE PRECISION TWO |
|
|
137 PARAMETER ( TWO = 2.0D0 ) |
|
|
138 * .. |
|
|
139 * .. Local Scalars .. |
|
|
140 INTEGER ICMAX, J |
|
|
141 DOUBLE PRECISION BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21, |
|
|
142 $ CI22, CMAX, CNORM, CR21, CR22, CSI, CSR, LI21, |
|
|
143 $ LR21, SMINI, SMLNUM, TEMP, U22ABS, UI11, UI11R, |
|
|
144 $ UI12, UI12S, UI22, UR11, UR11R, UR12, UR12S, |
|
|
145 $ UR22, XI1, XI2, XR1, XR2 |
|
|
146 * .. |
|
|
147 * .. Local Arrays .. |
|
|
148 LOGICAL RSWAP( 4 ), ZSWAP( 4 ) |
|
|
149 INTEGER IPIVOT( 4, 4 ) |
|
|
150 DOUBLE PRECISION CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 ) |
|
|
151 * .. |
|
|
152 * .. External Functions .. |
|
|
153 DOUBLE PRECISION DLAMCH |
|
|
154 EXTERNAL DLAMCH |
|
|
155 * .. |
|
|
156 * .. External Subroutines .. |
|
|
157 EXTERNAL DLADIV |
|
|
158 * .. |
|
|
159 * .. Intrinsic Functions .. |
|
|
160 INTRINSIC ABS, MAX |
|
|
161 * .. |
|
|
162 * .. Equivalences .. |
|
|
163 EQUIVALENCE ( CI( 1, 1 ), CIV( 1 ) ), |
|
|
164 $ ( CR( 1, 1 ), CRV( 1 ) ) |
|
|
165 * .. |
|
|
166 * .. Data statements .. |
|
|
167 DATA ZSWAP / .FALSE., .FALSE., .TRUE., .TRUE. / |
|
|
168 DATA RSWAP / .FALSE., .TRUE., .FALSE., .TRUE. / |
|
|
169 DATA IPIVOT / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4, |
|
|
170 $ 3, 2, 1 / |
|
|
171 * .. |
|
|
172 * .. Executable Statements .. |
|
|
173 * |
|
|
174 * Compute BIGNUM |
|
|
175 * |
|
|
176 SMLNUM = TWO*DLAMCH( 'Safe minimum' ) |
|
|
177 BIGNUM = ONE / SMLNUM |
|
|
178 SMINI = MAX( SMIN, SMLNUM ) |
|
|
179 * |
|
|
180 * Don't check for input errors |
|
|
181 * |
|
|
182 INFO = 0 |
|
|
183 * |
|
|
184 * Standard Initializations |
|
|
185 * |
|
|
186 SCALE = ONE |
|
|
187 * |
|
|
188 IF( NA.EQ.1 ) THEN |
|
|
189 * |
|
|
190 * 1 x 1 (i.e., scalar) system C X = B |
|
|
191 * |
|
|
192 IF( NW.EQ.1 ) THEN |
|
|
193 * |
|
|
194 * Real 1x1 system. |
|
|
195 * |
|
|
196 * C = ca A - w D |
|
|
197 * |
|
|
198 CSR = CA*A( 1, 1 ) - WR*D1 |
|
|
199 CNORM = ABS( CSR ) |
|
|
200 * |
|
|
201 * If | C | < SMINI, use C = SMINI |
|
|
202 * |
|
|
203 IF( CNORM.LT.SMINI ) THEN |
|
|
204 CSR = SMINI |
|
|
205 CNORM = SMINI |
|
|
206 INFO = 1 |
|
|
207 END IF |
|
|
208 * |
|
|
209 * Check scaling for X = B / C |
|
|
210 * |
|
|
211 BNORM = ABS( B( 1, 1 ) ) |
|
|
212 IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN |
|
|
213 IF( BNORM.GT.BIGNUM*CNORM ) |
|
|
214 $ SCALE = ONE / BNORM |
|
|
215 END IF |
|
|
216 * |
|
|
217 * Compute X |
|
|
218 * |
|
|
219 X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / CSR |
|
|
220 XNORM = ABS( X( 1, 1 ) ) |
|
|
221 ELSE |
|
|
222 * |
|
|
223 * Complex 1x1 system (w is complex) |
|
|
224 * |
|
|
225 * C = ca A - w D |
|
|
226 * |
|
|
227 CSR = CA*A( 1, 1 ) - WR*D1 |
|
|
228 CSI = -WI*D1 |
|
|
229 CNORM = ABS( CSR ) + ABS( CSI ) |
|
|
230 * |
|
|
231 * If | C | < SMINI, use C = SMINI |
|
|
232 * |
|
|
233 IF( CNORM.LT.SMINI ) THEN |
|
|
234 CSR = SMINI |
|
|
235 CSI = ZERO |
|
|
236 CNORM = SMINI |
|
|
237 INFO = 1 |
|
|
238 END IF |
|
|
239 * |
|
|
240 * Check scaling for X = B / C |
|
|
241 * |
|
|
242 BNORM = ABS( B( 1, 1 ) ) + ABS( B( 1, 2 ) ) |
|
|
243 IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN |
|
|
244 IF( BNORM.GT.BIGNUM*CNORM ) |
|
|
245 $ SCALE = ONE / BNORM |
|
|
246 END IF |
|
|
247 * |
|
|
248 * Compute X |
|
|
249 * |
|
|
250 CALL DLADIV( SCALE*B( 1, 1 ), SCALE*B( 1, 2 ), CSR, CSI, |
|
|
251 $ X( 1, 1 ), X( 1, 2 ) ) |
|
|
252 XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) ) |
|
|
253 END IF |
|
|
254 * |
|
|
255 ELSE |
|
|
256 * |
|
|
257 * 2x2 System |
|
|
258 * |
|
|
259 * Compute the real part of C = ca A - w D (or ca A' - w D ) |
|
|
260 * |
|
|
261 CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1 |
|
|
262 CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2 |
|
|
263 IF( LTRANS ) THEN |
|
|
264 CR( 1, 2 ) = CA*A( 2, 1 ) |
|
|
265 CR( 2, 1 ) = CA*A( 1, 2 ) |
|
|
266 ELSE |
|
|
267 CR( 2, 1 ) = CA*A( 2, 1 ) |
|
|
268 CR( 1, 2 ) = CA*A( 1, 2 ) |
|
|
269 END IF |
|
|
270 * |
|
|
271 IF( NW.EQ.1 ) THEN |
|
|
272 * |
|
|
273 * Real 2x2 system (w is real) |
|
|
274 * |
|
|
275 * Find the largest element in C |
|
|
276 * |
|
|
277 CMAX = ZERO |
|
|
278 ICMAX = 0 |
|
|
279 * |
|
|
280 DO 10 J = 1, 4 |
|
|
281 IF( ABS( CRV( J ) ).GT.CMAX ) THEN |
|
|
282 CMAX = ABS( CRV( J ) ) |
|
|
283 ICMAX = J |
|
|
284 END IF |
|
|
285 10 CONTINUE |
|
|
286 * |
|
|
287 * If norm(C) < SMINI, use SMINI*identity. |
|
|
288 * |
|
|
289 IF( CMAX.LT.SMINI ) THEN |
|
|
290 BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 2, 1 ) ) ) |
|
|
291 IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN |
|
|
292 IF( BNORM.GT.BIGNUM*SMINI ) |
|
|
293 $ SCALE = ONE / BNORM |
|
|
294 END IF |
|
|
295 TEMP = SCALE / SMINI |
|
|
296 X( 1, 1 ) = TEMP*B( 1, 1 ) |
|
|
297 X( 2, 1 ) = TEMP*B( 2, 1 ) |
|
|
298 XNORM = TEMP*BNORM |
|
|
299 INFO = 1 |
|
|
300 RETURN |
|
|
301 END IF |
|
|
302 * |
|
|
303 * Gaussian elimination with complete pivoting. |
|
|
304 * |
|
|
305 UR11 = CRV( ICMAX ) |
|
|
306 CR21 = CRV( IPIVOT( 2, ICMAX ) ) |
|
|
307 UR12 = CRV( IPIVOT( 3, ICMAX ) ) |
|
|
308 CR22 = CRV( IPIVOT( 4, ICMAX ) ) |
|
|
309 UR11R = ONE / UR11 |
|
|
310 LR21 = UR11R*CR21 |
|
|
311 UR22 = CR22 - UR12*LR21 |
|
|
312 * |
|
|
313 * If smaller pivot < SMINI, use SMINI |
|
|
314 * |
|
|
315 IF( ABS( UR22 ).LT.SMINI ) THEN |
|
|
316 UR22 = SMINI |
|
|
317 INFO = 1 |
|
|
318 END IF |
|
|
319 IF( RSWAP( ICMAX ) ) THEN |
|
|
320 BR1 = B( 2, 1 ) |
|
|
321 BR2 = B( 1, 1 ) |
|
|
322 ELSE |
|
|
323 BR1 = B( 1, 1 ) |
|
|
324 BR2 = B( 2, 1 ) |
|
|
325 END IF |
|
|
326 BR2 = BR2 - LR21*BR1 |
|
|
327 BBND = MAX( ABS( BR1*( UR22*UR11R ) ), ABS( BR2 ) ) |
|
|
328 IF( BBND.GT.ONE .AND. ABS( UR22 ).LT.ONE ) THEN |
|
|
329 IF( BBND.GE.BIGNUM*ABS( UR22 ) ) |
|
|
330 $ SCALE = ONE / BBND |
|
|
331 END IF |
|
|
332 * |
|
|
333 XR2 = ( BR2*SCALE ) / UR22 |
|
|
334 XR1 = ( SCALE*BR1 )*UR11R - XR2*( UR11R*UR12 ) |
|
|
335 IF( ZSWAP( ICMAX ) ) THEN |
|
|
336 X( 1, 1 ) = XR2 |
|
|
337 X( 2, 1 ) = XR1 |
|
|
338 ELSE |
|
|
339 X( 1, 1 ) = XR1 |
|
|
340 X( 2, 1 ) = XR2 |
|
|
341 END IF |
|
|
342 XNORM = MAX( ABS( XR1 ), ABS( XR2 ) ) |
|
|
343 * |
|
|
344 * Further scaling if norm(A) norm(X) > overflow |
|
|
345 * |
|
|
346 IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN |
|
|
347 IF( XNORM.GT.BIGNUM / CMAX ) THEN |
|
|
348 TEMP = CMAX / BIGNUM |
|
|
349 X( 1, 1 ) = TEMP*X( 1, 1 ) |
|
|
350 X( 2, 1 ) = TEMP*X( 2, 1 ) |
|
|
351 XNORM = TEMP*XNORM |
|
|
352 SCALE = TEMP*SCALE |
|
|
353 END IF |
|
|
354 END IF |
|
|
355 ELSE |
|
|
356 * |
|
|
357 * Complex 2x2 system (w is complex) |
|
|
358 * |
|
|
359 * Find the largest element in C |
|
|
360 * |
|
|
361 CI( 1, 1 ) = -WI*D1 |
|
|
362 CI( 2, 1 ) = ZERO |
|
|
363 CI( 1, 2 ) = ZERO |
|
|
364 CI( 2, 2 ) = -WI*D2 |
|
|
365 CMAX = ZERO |
|
|
366 ICMAX = 0 |
|
|
367 * |
|
|
368 DO 20 J = 1, 4 |
|
|
369 IF( ABS( CRV( J ) )+ABS( CIV( J ) ).GT.CMAX ) THEN |
|
|
370 CMAX = ABS( CRV( J ) ) + ABS( CIV( J ) ) |
|
|
371 ICMAX = J |
|
|
372 END IF |
|
|
373 20 CONTINUE |
|
|
374 * |
|
|
375 * If norm(C) < SMINI, use SMINI*identity. |
|
|
376 * |
|
|
377 IF( CMAX.LT.SMINI ) THEN |
|
|
378 BNORM = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ), |
|
|
379 $ ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) ) |
|
|
380 IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN |
|
|
381 IF( BNORM.GT.BIGNUM*SMINI ) |
|
|
382 $ SCALE = ONE / BNORM |
|
|
383 END IF |
|
|
384 TEMP = SCALE / SMINI |
|
|
385 X( 1, 1 ) = TEMP*B( 1, 1 ) |
|
|
386 X( 2, 1 ) = TEMP*B( 2, 1 ) |
|
|
387 X( 1, 2 ) = TEMP*B( 1, 2 ) |
|
|
388 X( 2, 2 ) = TEMP*B( 2, 2 ) |
|
|
389 XNORM = TEMP*BNORM |
|
|
390 INFO = 1 |
|
|
391 RETURN |
|
|
392 END IF |
|
|
393 * |
|
|
394 * Gaussian elimination with complete pivoting. |
|
|
395 * |
|
|
396 UR11 = CRV( ICMAX ) |
|
|
397 UI11 = CIV( ICMAX ) |
|
|
398 CR21 = CRV( IPIVOT( 2, ICMAX ) ) |
|
|
399 CI21 = CIV( IPIVOT( 2, ICMAX ) ) |
|
|
400 UR12 = CRV( IPIVOT( 3, ICMAX ) ) |
|
|
401 UI12 = CIV( IPIVOT( 3, ICMAX ) ) |
|
|
402 CR22 = CRV( IPIVOT( 4, ICMAX ) ) |
|
|
403 CI22 = CIV( IPIVOT( 4, ICMAX ) ) |
|
|
404 IF( ICMAX.EQ.1 .OR. ICMAX.EQ.4 ) THEN |
|
|
405 * |
|
|
406 * Code when off-diagonals of pivoted C are real |
|
|
407 * |
|
|
408 IF( ABS( UR11 ).GT.ABS( UI11 ) ) THEN |
|
|
409 TEMP = UI11 / UR11 |
|
|
410 UR11R = ONE / ( UR11*( ONE+TEMP**2 ) ) |
|
|
411 UI11R = -TEMP*UR11R |
|
|
412 ELSE |
|
|
413 TEMP = UR11 / UI11 |
|
|
414 UI11R = -ONE / ( UI11*( ONE+TEMP**2 ) ) |
|
|
415 UR11R = -TEMP*UI11R |
|
|
416 END IF |
|
|
417 LR21 = CR21*UR11R |
|
|
418 LI21 = CR21*UI11R |
|
|
419 UR12S = UR12*UR11R |
|
|
420 UI12S = UR12*UI11R |
|
|
421 UR22 = CR22 - UR12*LR21 |
|
|
422 UI22 = CI22 - UR12*LI21 |
|
|
423 ELSE |
|
|
424 * |
|
|
425 * Code when diagonals of pivoted C are real |
|
|
426 * |
|
|
427 UR11R = ONE / UR11 |
|
|
428 UI11R = ZERO |
|
|
429 LR21 = CR21*UR11R |
|
|
430 LI21 = CI21*UR11R |
|
|
431 UR12S = UR12*UR11R |
|
|
432 UI12S = UI12*UR11R |
|
|
433 UR22 = CR22 - UR12*LR21 + UI12*LI21 |
|
|
434 UI22 = -UR12*LI21 - UI12*LR21 |
|
|
435 END IF |
|
|
436 U22ABS = ABS( UR22 ) + ABS( UI22 ) |
|
|
437 * |
|
|
438 * If smaller pivot < SMINI, use SMINI |
|
|
439 * |
|
|
440 IF( U22ABS.LT.SMINI ) THEN |
|
|
441 UR22 = SMINI |
|
|
442 UI22 = ZERO |
|
|
443 INFO = 1 |
|
|
444 END IF |
|
|
445 IF( RSWAP( ICMAX ) ) THEN |
|
|
446 BR2 = B( 1, 1 ) |
|
|
447 BR1 = B( 2, 1 ) |
|
|
448 BI2 = B( 1, 2 ) |
|
|
449 BI1 = B( 2, 2 ) |
|
|
450 ELSE |
|
|
451 BR1 = B( 1, 1 ) |
|
|
452 BR2 = B( 2, 1 ) |
|
|
453 BI1 = B( 1, 2 ) |
|
|
454 BI2 = B( 2, 2 ) |
|
|
455 END IF |
|
|
456 BR2 = BR2 - LR21*BR1 + LI21*BI1 |
|
|
457 BI2 = BI2 - LI21*BR1 - LR21*BI1 |
|
|
458 BBND = MAX( ( ABS( BR1 )+ABS( BI1 ) )* |
|
|
459 $ ( U22ABS*( ABS( UR11R )+ABS( UI11R ) ) ), |
|
|
460 $ ABS( BR2 )+ABS( BI2 ) ) |
|
|
461 IF( BBND.GT.ONE .AND. U22ABS.LT.ONE ) THEN |
|
|
462 IF( BBND.GE.BIGNUM*U22ABS ) THEN |
|
|
463 SCALE = ONE / BBND |
|
|
464 BR1 = SCALE*BR1 |
|
|
465 BI1 = SCALE*BI1 |
|
|
466 BR2 = SCALE*BR2 |
|
|
467 BI2 = SCALE*BI2 |
|
|
468 END IF |
|
|
469 END IF |
|
|
470 * |
|
|
471 CALL DLADIV( BR2, BI2, UR22, UI22, XR2, XI2 ) |
|
|
472 XR1 = UR11R*BR1 - UI11R*BI1 - UR12S*XR2 + UI12S*XI2 |
|
|
473 XI1 = UI11R*BR1 + UR11R*BI1 - UI12S*XR2 - UR12S*XI2 |
|
|
474 IF( ZSWAP( ICMAX ) ) THEN |
|
|
475 X( 1, 1 ) = XR2 |
|
|
476 X( 2, 1 ) = XR1 |
|
|
477 X( 1, 2 ) = XI2 |
|
|
478 X( 2, 2 ) = XI1 |
|
|
479 ELSE |
|
|
480 X( 1, 1 ) = XR1 |
|
|
481 X( 2, 1 ) = XR2 |
|
|
482 X( 1, 2 ) = XI1 |
|
|
483 X( 2, 2 ) = XI2 |
|
|
484 END IF |
|
|
485 XNORM = MAX( ABS( XR1 )+ABS( XI1 ), ABS( XR2 )+ABS( XI2 ) ) |
|
|
486 * |
|
|
487 * Further scaling if norm(A) norm(X) > overflow |
|
|
488 * |
|
|
489 IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN |
|
|
490 IF( XNORM.GT.BIGNUM / CMAX ) THEN |
|
|
491 TEMP = CMAX / BIGNUM |
|
|
492 X( 1, 1 ) = TEMP*X( 1, 1 ) |
|
|
493 X( 2, 1 ) = TEMP*X( 2, 1 ) |
|
|
494 X( 1, 2 ) = TEMP*X( 1, 2 ) |
|
|
495 X( 2, 2 ) = TEMP*X( 2, 2 ) |
|
|
496 XNORM = TEMP*XNORM |
|
|
497 SCALE = TEMP*SCALE |
|
|
498 END IF |
|
|
499 END IF |
|
|
500 END IF |
|
|
501 END IF |
|
|
502 * |
|
|
503 RETURN |
|
|
504 * |
|
|
505 * End of DLALN2 |
|
|
506 * |
|
|
507 END |
