LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ slaln2()

subroutine slaln2 ( logical  ltrans,
integer  na,
integer  nw,
real  smin,
real  ca,
real, dimension( lda, * )  a,
integer  lda,
real  d1,
real  d2,
real, dimension( ldb, * )  b,
integer  ldb,
real  wr,
real  wi,
real, dimension( ldx, * )  x,
integer  ldx,
real  scale,
real  xnorm,
integer  info 
)

SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.

Download SLALN2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SLALN2 solves a system of the form  (ca A - w D ) X = s B
 or (ca A**T - w D) X = s B   with possible scaling ("s") and
 perturbation of A.  (A**T means A-transpose.)

 A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
 real diagonal matrix, w is a real or complex value, and X and B are
 NA x 1 matrices -- real if w is real, complex if w is complex.  NA
 may be 1 or 2.

 If w is complex, X and B are represented as NA x 2 matrices,
 the first column of each being the real part and the second
 being the imaginary part.

 "s" is a scaling factor (<= 1), computed by SLALN2, which is
 so chosen that X can be computed without overflow.  X is further
 scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
 than overflow.

 If both singular values of (ca A - w D) are less than SMIN,
 SMIN*identity will be used instead of (ca A - w D).  If only one
 singular value is less than SMIN, one element of (ca A - w D) will be
 perturbed enough to make the smallest singular value roughly SMIN.
 If both singular values are at least SMIN, (ca A - w D) will not be
 perturbed.  In any case, the perturbation will be at most some small
 multiple of max( SMIN, ulp*norm(ca A - w D) ).  The singular values
 are computed by infinity-norm approximations, and thus will only be
 correct to a factor of 2 or so.

 Note: all input quantities are assumed to be smaller than overflow
 by a reasonable factor.  (See BIGNUM.)
Parameters
[in]LTRANS
          LTRANS is LOGICAL
          =.TRUE.:  A-transpose will be used.
          =.FALSE.: A will be used (not transposed.)
[in]NA
          NA is INTEGER
          The size of the matrix A.  It may (only) be 1 or 2.
[in]NW
          NW is INTEGER
          1 if "w" is real, 2 if "w" is complex.  It may only be 1
          or 2.
[in]SMIN
          SMIN is REAL
          The desired lower bound on the singular values of A.  This
          should be a safe distance away from underflow or overflow,
          say, between (underflow/machine precision) and  (machine
          precision * overflow ).  (See BIGNUM and ULP.)
[in]CA
          CA is REAL
          The coefficient c, which A is multiplied by.
[in]A
          A is REAL array, dimension (LDA,NA)
          The NA x NA matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  It must be at least NA.
[in]D1
          D1 is REAL
          The 1,1 element in the diagonal matrix D.
[in]D2
          D2 is REAL
          The 2,2 element in the diagonal matrix D.  Not used if NA=1.
[in]B
          B is REAL array, dimension (LDB,NW)
          The NA x NW matrix B (right-hand side).  If NW=2 ("w" is
          complex), column 1 contains the real part of B and column 2
          contains the imaginary part.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  It must be at least NA.
[in]WR
          WR is REAL
          The real part of the scalar "w".
[in]WI
          WI is REAL
          The imaginary part of the scalar "w".  Not used if NW=1.
[out]X
          X is REAL array, dimension (LDX,NW)
          The NA x NW matrix X (unknowns), as computed by SLALN2.
          If NW=2 ("w" is complex), on exit, column 1 will contain
          the real part of X and column 2 will contain the imaginary
          part.
[in]LDX
          LDX is INTEGER
          The leading dimension of X.  It must be at least NA.
[out]SCALE
          SCALE is REAL
          The scale factor that B must be multiplied by to insure
          that overflow does not occur when computing X.  Thus,
          (ca A - w D) X  will be SCALE*B, not B (ignoring
          perturbations of A.)  It will be at most 1.
[out]XNORM
          XNORM is REAL
          The infinity-norm of X, when X is regarded as an NA x NW
          real matrix.
[out]INFO
          INFO is INTEGER
          An error flag.  It will be set to zero if no error occurs,
          a negative number if an argument is in error, or a positive
          number if  ca A - w D  had to be perturbed.
          The possible values are:
          = 0: No error occurred, and (ca A - w D) did not have to be
                 perturbed.
          = 1: (ca A - w D) had to be perturbed to make its smallest
               (or only) singular value greater than SMIN.
          NOTE: In the interests of speed, this routine does not
                check the inputs for errors.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 216 of file slaln2.f.

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 REAL CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
227* ..
228* .. Array Arguments ..
229 REAL A( LDA, * ), B( LDB, * ), X( LDX, * )
230* ..
231*
232* =====================================================================
233*
234* .. Parameters ..
235 REAL ZERO, ONE
236 parameter( zero = 0.0e0, one = 1.0e0 )
237 REAL TWO
238 parameter( two = 2.0e0 )
239* ..
240* .. Local Scalars ..
241 INTEGER ICMAX, J
242 REAL 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 CSWAP( 4 ), RSWAP( 4 )
250 INTEGER IPIVOT( 4, 4 )
251 REAL CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )
252* ..
253* .. External Functions ..
254 REAL SLAMCH
255 EXTERNAL slamch
256* ..
257* .. External Subroutines ..
258 EXTERNAL sladiv
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 cswap / .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*slamch( '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 sladiv( 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( cswap( 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 sladiv( 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( cswap( 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 SLALN2
607*
sladiv
subroutine sladiv(a, b, c, d, p, q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition sladiv.f:91
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: