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