LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
dlaln2.f
Go to the documentation of this file.
1*> \brief \b DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download DLALN2 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaln2.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaln2.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaln2.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B,
20* LDB, WR, WI, X, LDX, SCALE, XNORM, INFO )
21*
22* .. Scalar Arguments ..
23* LOGICAL LTRANS
24* INTEGER INFO, LDA, LDB, LDX, NA, NW
25* DOUBLE PRECISION CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
26* ..
27* .. Array Arguments ..
28* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> DLALN2 solves a system of the form (ca A - w D ) X = s B
38*> or (ca A**T - w D) X = s B with possible scaling ("s") and
39*> perturbation of A. (A**T means A-transpose.)
40*>
41*> A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
42*> real diagonal matrix, w is a real or complex value, and X and B are
43*> NA x 1 matrices -- real if w is real, complex if w is complex. NA
44*> may be 1 or 2.
45*>
46*> If w is complex, X and B are represented as NA x 2 matrices,
47*> the first column of each being the real part and the second
48*> being the imaginary part.
49*>
50*> "s" is a scaling factor (<= 1), computed by DLALN2, which is
51*> so chosen that X can be computed without overflow. X is further
52*> scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
53*> than overflow.
54*>
55*> If both singular values of (ca A - w D) are less than SMIN,
56*> SMIN*identity will be used instead of (ca A - w D). If only one
57*> singular value is less than SMIN, one element of (ca A - w D) will be
58*> perturbed enough to make the smallest singular value roughly SMIN.
59*> If both singular values are at least SMIN, (ca A - w D) will not be
60*> perturbed. In any case, the perturbation will be at most some small
61*> multiple of max( SMIN, ulp*norm(ca A - w D) ). The singular values
62*> are computed by infinity-norm approximations, and thus will only be
63*> correct to a factor of 2 or so.
64*>
65*> Note: all input quantities are assumed to be smaller than overflow
66*> by a reasonable factor. (See BIGNUM.)
67*> \endverbatim
68*
69* Arguments:
70* ==========
71*
72*> \param[in] LTRANS
73*> \verbatim
74*> LTRANS is LOGICAL
75*> =.TRUE.: A-transpose will be used.
76*> =.FALSE.: A will be used (not transposed.)
77*> \endverbatim
78*>
79*> \param[in] NA
80*> \verbatim
81*> NA is INTEGER
82*> The size of the matrix A. It may (only) be 1 or 2.
83*> \endverbatim
84*>
85*> \param[in] NW
86*> \verbatim
87*> NW is INTEGER
88*> 1 if "w" is real, 2 if "w" is complex. It may only be 1
89*> or 2.
90*> \endverbatim
91*>
92*> \param[in] SMIN
93*> \verbatim
94*> SMIN is DOUBLE PRECISION
95*> The desired lower bound on the singular values of A. This
96*> should be a safe distance away from underflow or overflow,
97*> say, between (underflow/machine precision) and (machine
98*> precision * overflow ). (See BIGNUM and ULP.)
99*> \endverbatim
100*>
101*> \param[in] CA
102*> \verbatim
103*> CA is DOUBLE PRECISION
104*> The coefficient c, which A is multiplied by.
105*> \endverbatim
106*>
107*> \param[in] A
108*> \verbatim
109*> A is DOUBLE PRECISION array, dimension (LDA,NA)
110*> The NA x NA matrix A.
111*> \endverbatim
112*>
113*> \param[in] LDA
114*> \verbatim
115*> LDA is INTEGER
116*> The leading dimension of A. It must be at least NA.
117*> \endverbatim
118*>
119*> \param[in] D1
120*> \verbatim
121*> D1 is DOUBLE PRECISION
122*> The 1,1 element in the diagonal matrix D.
123*> \endverbatim
124*>
125*> \param[in] D2
126*> \verbatim
127*> D2 is DOUBLE PRECISION
128*> The 2,2 element in the diagonal matrix D. Not used if NA=1.
129*> \endverbatim
130*>
131*> \param[in] B
132*> \verbatim
133*> B is DOUBLE PRECISION array, dimension (LDB,NW)
134*> The NA x NW matrix B (right-hand side). If NW=2 ("w" is
135*> complex), column 1 contains the real part of B and column 2
136*> contains the imaginary part.
137*> \endverbatim
138*>
139*> \param[in] LDB
140*> \verbatim
141*> LDB is INTEGER
142*> The leading dimension of B. It must be at least NA.
143*> \endverbatim
144*>
145*> \param[in] WR
146*> \verbatim
147*> WR is DOUBLE PRECISION
148*> The real part of the scalar "w".
149*> \endverbatim
150*>
151*> \param[in] WI
152*> \verbatim
153*> WI is DOUBLE PRECISION
154*> The imaginary part of the scalar "w". Not used if NW=1.
155*> \endverbatim
156*>
157*> \param[out] X
158*> \verbatim
159*> X is DOUBLE PRECISION array, dimension (LDX,NW)
160*> The NA x NW matrix X (unknowns), as computed by DLALN2.
161*> If NW=2 ("w" is complex), on exit, column 1 will contain
162*> the real part of X and column 2 will contain the imaginary
163*> part.
164*> \endverbatim
165*>
166*> \param[in] LDX
167*> \verbatim
168*> LDX is INTEGER
169*> The leading dimension of X. It must be at least NA.
170*> \endverbatim
171*>
172*> \param[out] SCALE
173*> \verbatim
174*> SCALE is DOUBLE PRECISION
175*> The scale factor that B must be multiplied by to insure
176*> that overflow does not occur when computing X. Thus,
177*> (ca A - w D) X will be SCALE*B, not B (ignoring
178*> perturbations of A.) It will be at most 1.
179*> \endverbatim
180*>
181*> \param[out] XNORM
182*> \verbatim
183*> XNORM is DOUBLE PRECISION
184*> The infinity-norm of X, when X is regarded as an NA x NW
185*> real matrix.
186*> \endverbatim
187*>
188*> \param[out] INFO
189*> \verbatim
190*> INFO is INTEGER
191*> An error flag. It will be set to zero if no error occurs,
192*> a negative number if an argument is in error, or a positive
193*> number if ca A - w D had to be perturbed.
194*> The possible values are:
195*> = 0: No error occurred, and (ca A - w D) did not have to be
196*> perturbed.
197*> = 1: (ca A - w D) had to be perturbed to make its smallest
198*> (or only) singular value greater than SMIN.
199*> NOTE: In the interests of speed, this routine does not
200*> check the inputs for errors.
201*> \endverbatim
202*
203* Authors:
204* ========
205*
206*> \author Univ. of Tennessee
207*> \author Univ. of California Berkeley
208*> \author Univ. of Colorado Denver
209*> \author NAG Ltd.
210*
211*> \ingroup laln2
212*
213* =====================================================================
214 SUBROUTINE dlaln2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B,
215 $ LDB, WR, WI, X, LDX, SCALE, XNORM, INFO )
216*
217* -- LAPACK auxiliary routine --
218* -- LAPACK is a software package provided by Univ. of Tennessee, --
219* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220*
221* .. Scalar Arguments ..
222 LOGICAL LTRANS
223 INTEGER INFO, LDA, LDB, LDX, NA, NW
224 DOUBLE PRECISION CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
225* ..
226* .. Array Arguments ..
227 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * )
228* ..
229*
230* =====================================================================
231*
232* .. Parameters ..
233 DOUBLE PRECISION ZERO, ONE
234 parameter( zero = 0.0d0, one = 1.0d0 )
235 DOUBLE PRECISION TWO
236 parameter( two = 2.0d0 )
237* ..
238* .. Local Scalars ..
239 INTEGER ICMAX, J
240 DOUBLE PRECISION BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21,
241 $ ci22, cmax, cnorm, cr21, cr22, csi, csr, li21,
242 $ lr21, smini, smlnum, temp, u22abs, ui11, ui11r,
243 $ ui12, ui12s, ui22, ur11, ur11r, ur12, ur12s,
244 $ ur22, xi1, xi2, xr1, xr2
245* ..
246* .. Local Arrays ..
247 LOGICAL RSWAP( 4 ), ZSWAP( 4 )
248 INTEGER IPIVOT( 4, 4 )
249 DOUBLE PRECISION CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )
250* ..
251* .. External Functions ..
252 DOUBLE PRECISION DLAMCH
253 EXTERNAL dlamch
254* ..
255* .. External Subroutines ..
256 EXTERNAL dladiv
257* ..
258* .. Intrinsic Functions ..
259 INTRINSIC abs, max
260* ..
261* .. Equivalences ..
262 equivalence( ci( 1, 1 ), civ( 1 ) ),
263 $ ( cr( 1, 1 ), crv( 1 ) )
264* ..
265* .. Data statements ..
266 DATA zswap / .false., .false., .true., .true. /
267 DATA rswap / .false., .true., .false., .true. /
268 DATA ipivot / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4,
269 $ 3, 2, 1 /
270* ..
271* .. Executable Statements ..
272*
273* Compute BIGNUM
274*
275 smlnum = two*dlamch( 'Safe minimum' )
276 bignum = one / smlnum
277 smini = max( smin, smlnum )
278*
279* Don't check for input errors
280*
281 info = 0
282*
283* Standard Initializations
284*
285 scale = one
286*
287 IF( na.EQ.1 ) THEN
288*
289* 1 x 1 (i.e., scalar) system C X = B
290*
291 IF( nw.EQ.1 ) THEN
292*
293* Real 1x1 system.
294*
295* C = ca A - w D
296*
297 csr = ca*a( 1, 1 ) - wr*d1
298 cnorm = abs( csr )
299*
300* If | C | < SMINI, use C = SMINI
301*
302 IF( cnorm.LT.smini ) THEN
303 csr = smini
304 cnorm = smini
305 info = 1
306 END IF
307*
308* Check scaling for X = B / C
309*
310 bnorm = abs( b( 1, 1 ) )
311 IF( cnorm.LT.one .AND. bnorm.GT.one ) THEN
312 IF( bnorm.GT.bignum*cnorm )
313 $ scale = one / bnorm
314 END IF
315*
316* Compute X
317*
318 x( 1, 1 ) = ( b( 1, 1 )*scale ) / csr
319 xnorm = abs( x( 1, 1 ) )
320 ELSE
321*
322* Complex 1x1 system (w is complex)
323*
324* C = ca A - w D
325*
326 csr = ca*a( 1, 1 ) - wr*d1
327 csi = -wi*d1
328 cnorm = abs( csr ) + abs( csi )
329*
330* If | C | < SMINI, use C = SMINI
331*
332 IF( cnorm.LT.smini ) THEN
333 csr = smini
334 csi = zero
335 cnorm = smini
336 info = 1
337 END IF
338*
339* Check scaling for X = B / C
340*
341 bnorm = abs( b( 1, 1 ) ) + abs( b( 1, 2 ) )
342 IF( cnorm.LT.one .AND. bnorm.GT.one ) THEN
343 IF( bnorm.GT.bignum*cnorm )
344 $ scale = one / bnorm
345 END IF
346*
347* Compute X
348*
349 CALL dladiv( scale*b( 1, 1 ), scale*b( 1, 2 ), csr, csi,
350 $ x( 1, 1 ), x( 1, 2 ) )
351 xnorm = abs( x( 1, 1 ) ) + abs( x( 1, 2 ) )
352 END IF
353*
354 ELSE
355*
356* 2x2 System
357*
358* Compute the real part of C = ca A - w D (or ca A**T - w D )
359*
360 cr( 1, 1 ) = ca*a( 1, 1 ) - wr*d1
361 cr( 2, 2 ) = ca*a( 2, 2 ) - wr*d2
362 IF( ltrans ) THEN
363 cr( 1, 2 ) = ca*a( 2, 1 )
364 cr( 2, 1 ) = ca*a( 1, 2 )
365 ELSE
366 cr( 2, 1 ) = ca*a( 2, 1 )
367 cr( 1, 2 ) = ca*a( 1, 2 )
368 END IF
369*
370 IF( nw.EQ.1 ) THEN
371*
372* Real 2x2 system (w is real)
373*
374* Find the largest element in C
375*
376 cmax = zero
377 icmax = 0
378*
379 DO 10 j = 1, 4
380 IF( abs( crv( j ) ).GT.cmax ) THEN
381 cmax = abs( crv( j ) )
382 icmax = j
383 END IF
384 10 CONTINUE
385*
386* If norm(C) < SMINI, use SMINI*identity.
387*
388 IF( cmax.LT.smini ) THEN
389 bnorm = max( abs( b( 1, 1 ) ), abs( b( 2, 1 ) ) )
390 IF( smini.LT.one .AND. bnorm.GT.one ) THEN
391 IF( bnorm.GT.bignum*smini )
392 $ scale = one / bnorm
393 END IF
394 temp = scale / smini
395 x( 1, 1 ) = temp*b( 1, 1 )
396 x( 2, 1 ) = temp*b( 2, 1 )
397 xnorm = temp*bnorm
398 info = 1
399 RETURN
400 END IF
401*
402* Gaussian elimination with complete pivoting.
403*
404 ur11 = crv( icmax )
405 cr21 = crv( ipivot( 2, icmax ) )
406 ur12 = crv( ipivot( 3, icmax ) )
407 cr22 = crv( ipivot( 4, icmax ) )
408 ur11r = one / ur11
409 lr21 = ur11r*cr21
410 ur22 = cr22 - ur12*lr21
411*
412* If smaller pivot < SMINI, use SMINI
413*
414 IF( abs( ur22 ).LT.smini ) THEN
415 ur22 = smini
416 info = 1
417 END IF
418 IF( rswap( icmax ) ) THEN
419 br1 = b( 2, 1 )
420 br2 = b( 1, 1 )
421 ELSE
422 br1 = b( 1, 1 )
423 br2 = b( 2, 1 )
424 END IF
425 br2 = br2 - lr21*br1
426 bbnd = max( abs( br1*( ur22*ur11r ) ), abs( br2 ) )
427 IF( bbnd.GT.one .AND. abs( ur22 ).LT.one ) THEN
428 IF( bbnd.GE.bignum*abs( ur22 ) )
429 $ scale = one / bbnd
430 END IF
431*
432 xr2 = ( br2*scale ) / ur22
433 xr1 = ( scale*br1 )*ur11r - xr2*( ur11r*ur12 )
434 IF( zswap( icmax ) ) THEN
435 x( 1, 1 ) = xr2
436 x( 2, 1 ) = xr1
437 ELSE
438 x( 1, 1 ) = xr1
439 x( 2, 1 ) = xr2
440 END IF
441 xnorm = max( abs( xr1 ), abs( xr2 ) )
442*
443* Further scaling if norm(A) norm(X) > overflow
444*
445 IF( xnorm.GT.one .AND. cmax.GT.one ) THEN
446 IF( xnorm.GT.bignum / cmax ) THEN
447 temp = cmax / bignum
448 x( 1, 1 ) = temp*x( 1, 1 )
449 x( 2, 1 ) = temp*x( 2, 1 )
450 xnorm = temp*xnorm
451 scale = temp*scale
452 END IF
453 END IF
454 ELSE
455*
456* Complex 2x2 system (w is complex)
457*
458* Find the largest element in C
459*
460 ci( 1, 1 ) = -wi*d1
461 ci( 2, 1 ) = zero
462 ci( 1, 2 ) = zero
463 ci( 2, 2 ) = -wi*d2
464 cmax = zero
465 icmax = 0
466*
467 DO 20 j = 1, 4
468 IF( abs( crv( j ) )+abs( civ( j ) ).GT.cmax ) THEN
469 cmax = abs( crv( j ) ) + abs( civ( j ) )
470 icmax = j
471 END IF
472 20 CONTINUE
473*
474* If norm(C) < SMINI, use SMINI*identity.
475*
476 IF( cmax.LT.smini ) THEN
477 bnorm = max( abs( b( 1, 1 ) )+abs( b( 1, 2 ) ),
478 $ abs( b( 2, 1 ) )+abs( b( 2, 2 ) ) )
479 IF( smini.LT.one .AND. bnorm.GT.one ) THEN
480 IF( bnorm.GT.bignum*smini )
481 $ scale = one / bnorm
482 END IF
483 temp = scale / smini
484 x( 1, 1 ) = temp*b( 1, 1 )
485 x( 2, 1 ) = temp*b( 2, 1 )
486 x( 1, 2 ) = temp*b( 1, 2 )
487 x( 2, 2 ) = temp*b( 2, 2 )
488 xnorm = temp*bnorm
489 info = 1
490 RETURN
491 END IF
492*
493* Gaussian elimination with complete pivoting.
494*
495 ur11 = crv( icmax )
496 ui11 = civ( icmax )
497 cr21 = crv( ipivot( 2, icmax ) )
498 ci21 = civ( ipivot( 2, icmax ) )
499 ur12 = crv( ipivot( 3, icmax ) )
500 ui12 = civ( ipivot( 3, icmax ) )
501 cr22 = crv( ipivot( 4, icmax ) )
502 ci22 = civ( ipivot( 4, icmax ) )
503 IF( icmax.EQ.1 .OR. icmax.EQ.4 ) THEN
504*
505* Code when off-diagonals of pivoted C are real
506*
507 IF( abs( ur11 ).GT.abs( ui11 ) ) THEN
508 temp = ui11 / ur11
509 ur11r = one / ( ur11*( one+temp**2 ) )
510 ui11r = -temp*ur11r
511 ELSE
512 temp = ur11 / ui11
513 ui11r = -one / ( ui11*( one+temp**2 ) )
514 ur11r = -temp*ui11r
515 END IF
516 lr21 = cr21*ur11r
517 li21 = cr21*ui11r
518 ur12s = ur12*ur11r
519 ui12s = ur12*ui11r
520 ur22 = cr22 - ur12*lr21
521 ui22 = ci22 - ur12*li21
522 ELSE
523*
524* Code when diagonals of pivoted C are real
525*
526 ur11r = one / ur11
527 ui11r = zero
528 lr21 = cr21*ur11r
529 li21 = ci21*ur11r
530 ur12s = ur12*ur11r
531 ui12s = ui12*ur11r
532 ur22 = cr22 - ur12*lr21 + ui12*li21
533 ui22 = -ur12*li21 - ui12*lr21
534 END IF
535 u22abs = abs( ur22 ) + abs( ui22 )
536*
537* If smaller pivot < SMINI, use SMINI
538*
539 IF( u22abs.LT.smini ) THEN
540 ur22 = smini
541 ui22 = zero
542 info = 1
543 END IF
544 IF( rswap( icmax ) ) THEN
545 br2 = b( 1, 1 )
546 br1 = b( 2, 1 )
547 bi2 = b( 1, 2 )
548 bi1 = b( 2, 2 )
549 ELSE
550 br1 = b( 1, 1 )
551 br2 = b( 2, 1 )
552 bi1 = b( 1, 2 )
553 bi2 = b( 2, 2 )
554 END IF
555 br2 = br2 - lr21*br1 + li21*bi1
556 bi2 = bi2 - li21*br1 - lr21*bi1
557 bbnd = max( ( abs( br1 )+abs( bi1 ) )*
558 $ ( u22abs*( abs( ur11r )+abs( ui11r ) ) ),
559 $ abs( br2 )+abs( bi2 ) )
560 IF( bbnd.GT.one .AND. u22abs.LT.one ) THEN
561 IF( bbnd.GE.bignum*u22abs ) THEN
562 scale = one / bbnd
563 br1 = scale*br1
564 bi1 = scale*bi1
565 br2 = scale*br2
566 bi2 = scale*bi2
567 END IF
568 END IF
569*
570 CALL dladiv( br2, bi2, ur22, ui22, xr2, xi2 )
571 xr1 = ur11r*br1 - ui11r*bi1 - ur12s*xr2 + ui12s*xi2
572 xi1 = ui11r*br1 + ur11r*bi1 - ui12s*xr2 - ur12s*xi2
573 IF( zswap( icmax ) ) THEN
574 x( 1, 1 ) = xr2
575 x( 2, 1 ) = xr1
576 x( 1, 2 ) = xi2
577 x( 2, 2 ) = xi1
578 ELSE
579 x( 1, 1 ) = xr1
580 x( 2, 1 ) = xr2
581 x( 1, 2 ) = xi1
582 x( 2, 2 ) = xi2
583 END IF
584 xnorm = max( abs( xr1 )+abs( xi1 ), abs( xr2 )+abs( xi2 ) )
585*
586* Further scaling if norm(A) norm(X) > overflow
587*
588 IF( xnorm.GT.one .AND. cmax.GT.one ) THEN
589 IF( xnorm.GT.bignum / cmax ) THEN
590 temp = cmax / bignum
591 x( 1, 1 ) = temp*x( 1, 1 )
592 x( 2, 1 ) = temp*x( 2, 1 )
593 x( 1, 2 ) = temp*x( 1, 2 )
594 x( 2, 2 ) = temp*x( 2, 2 )
595 xnorm = temp*xnorm
596 scale = temp*scale
597 END IF
598 END IF
599 END IF
600 END IF
601*
602 RETURN
603*
604* End of DLALN2
605*
606 END
dladiv
subroutine dladiv(a, b, c, d, p, q)
DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition dladiv.f:89
dlaln2
subroutine dlaln2(ltrans, na, nw, smin, ca, a, lda, d1, d2, b, ldb, wr, wi, x, ldx, scale, xnorm, info)
DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
Definition dlaln2.f:216