dlaln2()

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

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

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

Purpose:

!>
!> DLALN2 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 () 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.
!>
!>  is a scaling factor (<= 1), computed by DLALN2, 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 is real, 2 if is complex. It may only be 1 !> or 2. !>
[in]	SMIN	!> SMIN is DOUBLE PRECISION !> 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 DOUBLE PRECISION !> The coefficient c, which A is multiplied by. !>
[in]	A	!> A is DOUBLE PRECISION 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 DOUBLE PRECISION !> The 1,1 element in the diagonal matrix D. !>
[in]	D2	!> D2 is DOUBLE PRECISION !> The 2,2 element in the diagonal matrix D. Not used if NA=1. !>
[in]	B	!> B is DOUBLE PRECISION array, dimension (LDB,NW) !> The NA x NW matrix B (right-hand side). If NW=2 ( 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 DOUBLE PRECISION !> The real part of the scalar . !>
[in]	WI	!> WI is DOUBLE PRECISION !> The imaginary part of the scalar . Not used if NW=1. !>
[out]	X	!> X is DOUBLE PRECISION array, dimension (LDX,NW) !> The NA x NW matrix X (unknowns), as computed by DLALN2. !> If NW=2 ( 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 DOUBLE PRECISION !> 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 DOUBLE PRECISION !> 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 214 of file dlaln2.f.

*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      LOGICAL            LTRANS
      INTEGER            INFO, LDA, LDB, LDX, NA, NW
      DOUBLE PRECISION   CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), X( LDX, * )
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d0, one = 1.0d0 )
      DOUBLE PRECISION   TWO
      parameter( two = 2.0d0 )
*     ..
*     .. Local Scalars ..
      INTEGER            ICMAX, J
      DOUBLE PRECISION   BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21,
     $                   CI22, CMAX, CNORM, CR21, CR22, CSI, CSR, LI21,
     $                   LR21, SMINI, SMLNUM, TEMP, U22ABS, UI11, UI11R,
     $                   UI12, UI12S, UI22, UR11, UR11R, UR12, UR12S,
     $                   UR22, XI1, XI2, XR1, XR2
*     ..
*     .. Local Arrays ..
      LOGICAL            RSWAP( 4 ), ZSWAP( 4 )
      INTEGER            IPIVOT( 4, 4 )
      DOUBLE PRECISION   CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           dlamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           dladiv
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max
*     ..
*     .. Equivalences ..
      equivalence( ci( 1, 1 ), civ( 1 ) ),
     $                   ( cr( 1, 1 ), crv( 1 ) )
*     ..
*     .. Data statements ..
      DATA               zswap / .false., .false., .true., .true. /
      DATA               rswap / .false., .true., .false., .true. /
      DATA               ipivot / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4,
     $                   3, 2, 1 /
*     ..
*     .. Executable Statements ..
*
*     Compute BIGNUM
*
      smlnum = two*dlamch( 'Safe minimum' )
      bignum = one / smlnum
      smini = max( smin, smlnum )
*
*     Don't check for input errors
*
      info = 0
*
*     Standard Initializations
*
      scale = one
*
      IF( na.EQ.1 ) THEN
*
*        1 x 1  (i.e., scalar) system   C X = B
*
         IF( nw.EQ.1 ) THEN
*
*           Real 1x1 system.
*
*           C = ca A - w D
*
            csr = ca*a( 1, 1 ) - wr*d1
            cnorm = abs( csr )
*
*           If | C | < SMINI, use C = SMINI
*
            IF( cnorm.LT.smini ) THEN
               csr = smini
               cnorm = smini
               info = 1
            END IF
*
*           Check scaling for  X = B / C
*
            bnorm = abs( b( 1, 1 ) )
            IF( cnorm.LT.one .AND. bnorm.GT.one ) THEN
               IF( bnorm.GT.bignum*cnorm )
     $            scale = one / bnorm
            END IF
*
*           Compute X
*
            x( 1, 1 ) = ( b( 1, 1 )*scale ) / csr
            xnorm = abs( x( 1, 1 ) )
         ELSE
*
*           Complex 1x1 system (w is complex)
*
*           C = ca A - w D
*
            csr = ca*a( 1, 1 ) - wr*d1
            csi = -wi*d1
            cnorm = abs( csr ) + abs( csi )
*
*           If | C | < SMINI, use C = SMINI
*
            IF( cnorm.LT.smini ) THEN
               csr = smini
               csi = zero
               cnorm = smini
               info = 1
            END IF
*
*           Check scaling for  X = B / C
*
            bnorm = abs( b( 1, 1 ) ) + abs( b( 1, 2 ) )
            IF( cnorm.LT.one .AND. bnorm.GT.one ) THEN
               IF( bnorm.GT.bignum*cnorm )
     $            scale = one / bnorm
            END IF
*
*           Compute X
*
            CALL dladiv( scale*b( 1, 1 ), scale*b( 1, 2 ), csr, csi,
     $                   x( 1, 1 ), x( 1, 2 ) )
            xnorm = abs( x( 1, 1 ) ) + abs( x( 1, 2 ) )
         END IF
*
      ELSE
*
*        2x2 System
*
*        Compute the real part of  C = ca A - w D  (or  ca A**T - w D )
*
         cr( 1, 1 ) = ca*a( 1, 1 ) - wr*d1
         cr( 2, 2 ) = ca*a( 2, 2 ) - wr*d2
         IF( ltrans ) THEN
            cr( 1, 2 ) = ca*a( 2, 1 )
            cr( 2, 1 ) = ca*a( 1, 2 )
         ELSE
            cr( 2, 1 ) = ca*a( 2, 1 )
            cr( 1, 2 ) = ca*a( 1, 2 )
         END IF
*
         IF( nw.EQ.1 ) THEN
*
*           Real 2x2 system  (w is real)
*
*           Find the largest element in C
*
            cmax = zero
            icmax = 0
*
            DO 10 j = 1, 4
               IF( abs( crv( j ) ).GT.cmax ) THEN
                  cmax = abs( crv( j ) )
                  icmax = j
               END IF
   10       CONTINUE
*
*           If norm(C) < SMINI, use SMINI*identity.
*
            IF( cmax.LT.smini ) THEN
               bnorm = max( abs( b( 1, 1 ) ), abs( b( 2, 1 ) ) )
               IF( smini.LT.one .AND. bnorm.GT.one ) THEN
                  IF( bnorm.GT.bignum*smini )
     $               scale = one / bnorm
               END IF
               temp = scale / smini
               x( 1, 1 ) = temp*b( 1, 1 )
               x( 2, 1 ) = temp*b( 2, 1 )
               xnorm = temp*bnorm
               info = 1
               RETURN
            END IF
*
*           Gaussian elimination with complete pivoting.
*
            ur11 = crv( icmax )
            cr21 = crv( ipivot( 2, icmax ) )
            ur12 = crv( ipivot( 3, icmax ) )
            cr22 = crv( ipivot( 4, icmax ) )
            ur11r = one / ur11
            lr21 = ur11r*cr21
            ur22 = cr22 - ur12*lr21
*
*           If smaller pivot < SMINI, use SMINI
*
            IF( abs( ur22 ).LT.smini ) THEN
               ur22 = smini
               info = 1
            END IF
            IF( rswap( icmax ) ) THEN
               br1 = b( 2, 1 )
               br2 = b( 1, 1 )
            ELSE
               br1 = b( 1, 1 )
               br2 = b( 2, 1 )
            END IF
            br2 = br2 - lr21*br1
            bbnd = max( abs( br1*( ur22*ur11r ) ), abs( br2 ) )
            IF( bbnd.GT.one .AND. abs( ur22 ).LT.one ) THEN
               IF( bbnd.GE.bignum*abs( ur22 ) )
     $            scale = one / bbnd
            END IF
*
            xr2 = ( br2*scale ) / ur22
            xr1 = ( scale*br1 )*ur11r - xr2*( ur11r*ur12 )
            IF( zswap( icmax ) ) THEN
               x( 1, 1 ) = xr2
               x( 2, 1 ) = xr1
            ELSE
               x( 1, 1 ) = xr1
               x( 2, 1 ) = xr2
            END IF
            xnorm = max( abs( xr1 ), abs( xr2 ) )
*
*           Further scaling if  norm(A) norm(X) > overflow
*
            IF( xnorm.GT.one .AND. cmax.GT.one ) THEN
               IF( xnorm.GT.bignum / cmax ) THEN
                  temp = cmax / bignum
                  x( 1, 1 ) = temp*x( 1, 1 )
                  x( 2, 1 ) = temp*x( 2, 1 )
                  xnorm = temp*xnorm
                  scale = temp*scale
               END IF
            END IF
         ELSE
*
*           Complex 2x2 system  (w is complex)
*
*           Find the largest element in C
*
            ci( 1, 1 ) = -wi*d1
            ci( 2, 1 ) = zero
            ci( 1, 2 ) = zero
            ci( 2, 2 ) = -wi*d2
            cmax = zero
            icmax = 0
*
            DO 20 j = 1, 4
               IF( abs( crv( j ) )+abs( civ( j ) ).GT.cmax ) THEN
                  cmax = abs( crv( j ) ) + abs( civ( j ) )
                  icmax = j
               END IF
   20       CONTINUE
*
*           If norm(C) < SMINI, use SMINI*identity.
*
            IF( cmax.LT.smini ) THEN
               bnorm = max( abs( b( 1, 1 ) )+abs( b( 1, 2 ) ),
     $                 abs( b( 2, 1 ) )+abs( b( 2, 2 ) ) )
               IF( smini.LT.one .AND. bnorm.GT.one ) THEN
                  IF( bnorm.GT.bignum*smini )
     $               scale = one / bnorm
               END IF
               temp = scale / smini
               x( 1, 1 ) = temp*b( 1, 1 )
               x( 2, 1 ) = temp*b( 2, 1 )
               x( 1, 2 ) = temp*b( 1, 2 )
               x( 2, 2 ) = temp*b( 2, 2 )
               xnorm = temp*bnorm
               info = 1
               RETURN
            END IF
*
*           Gaussian elimination with complete pivoting.
*
            ur11 = crv( icmax )
            ui11 = civ( icmax )
            cr21 = crv( ipivot( 2, icmax ) )
            ci21 = civ( ipivot( 2, icmax ) )
            ur12 = crv( ipivot( 3, icmax ) )
            ui12 = civ( ipivot( 3, icmax ) )
            cr22 = crv( ipivot( 4, icmax ) )
            ci22 = civ( ipivot( 4, icmax ) )
            IF( icmax.EQ.1 .OR. icmax.EQ.4 ) THEN
*
*              Code when off-diagonals of pivoted C are real
*
               IF( abs( ur11 ).GT.abs( ui11 ) ) THEN
                  temp = ui11 / ur11
                  ur11r = one / ( ur11*( one+temp**2 ) )
                  ui11r = -temp*ur11r
               ELSE
                  temp = ur11 / ui11
                  ui11r = -one / ( ui11*( one+temp**2 ) )
                  ur11r = -temp*ui11r
               END IF
               lr21 = cr21*ur11r
               li21 = cr21*ui11r
               ur12s = ur12*ur11r
               ui12s = ur12*ui11r
               ur22 = cr22 - ur12*lr21
               ui22 = ci22 - ur12*li21
            ELSE
*
*              Code when diagonals of pivoted C are real
*
               ur11r = one / ur11
               ui11r = zero
               lr21 = cr21*ur11r
               li21 = ci21*ur11r
               ur12s = ur12*ur11r
               ui12s = ui12*ur11r
               ur22 = cr22 - ur12*lr21 + ui12*li21
               ui22 = -ur12*li21 - ui12*lr21
            END IF
            u22abs = abs( ur22 ) + abs( ui22 )
*
*           If smaller pivot < SMINI, use SMINI
*
            IF( u22abs.LT.smini ) THEN
               ur22 = smini
               ui22 = zero
               info = 1
            END IF
            IF( rswap( icmax ) ) THEN
               br2 = b( 1, 1 )
               br1 = b( 2, 1 )
               bi2 = b( 1, 2 )
               bi1 = b( 2, 2 )
            ELSE
               br1 = b( 1, 1 )
               br2 = b( 2, 1 )
               bi1 = b( 1, 2 )
               bi2 = b( 2, 2 )
            END IF
            br2 = br2 - lr21*br1 + li21*bi1
            bi2 = bi2 - li21*br1 - lr21*bi1
            bbnd = max( ( abs( br1 )+abs( bi1 ) )*
     $             ( u22abs*( abs( ur11r )+abs( ui11r ) ) ),
     $             abs( br2 )+abs( bi2 ) )
            IF( bbnd.GT.one .AND. u22abs.LT.one ) THEN
               IF( bbnd.GE.bignum*u22abs ) THEN
                  scale = one / bbnd
                  br1 = scale*br1
                  bi1 = scale*bi1
                  br2 = scale*br2
                  bi2 = scale*bi2
               END IF
            END IF
*
            CALL dladiv( br2, bi2, ur22, ui22, xr2, xi2 )
            xr1 = ur11r*br1 - ui11r*bi1 - ur12s*xr2 + ui12s*xi2
            xi1 = ui11r*br1 + ur11r*bi1 - ui12s*xr2 - ur12s*xi2
            IF( zswap( icmax ) ) THEN
               x( 1, 1 ) = xr2
               x( 2, 1 ) = xr1
               x( 1, 2 ) = xi2
               x( 2, 2 ) = xi1
            ELSE
               x( 1, 1 ) = xr1
               x( 2, 1 ) = xr2
               x( 1, 2 ) = xi1
               x( 2, 2 ) = xi2
            END IF
            xnorm = max( abs( xr1 )+abs( xi1 ), abs( xr2 )+abs( xi2 ) )
*
*           Further scaling if  norm(A) norm(X) > overflow
*
            IF( xnorm.GT.one .AND. cmax.GT.one ) THEN
               IF( xnorm.GT.bignum / cmax ) THEN
                  temp = cmax / bignum
                  x( 1, 1 ) = temp*x( 1, 1 )
                  x( 2, 1 ) = temp*x( 2, 1 )
                  x( 1, 2 ) = temp*x( 1, 2 )
                  x( 2, 2 ) = temp*x( 2, 2 )
                  xnorm = temp*xnorm
                  scale = temp*scale
               END IF
            END IF
         END IF
      END IF
*
      RETURN
*
*     End of DLALN2
*

Here is the call graph for this function:

Here is the caller graph for this function: