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