508 lines
16 KiB
FortranFixed
508 lines
16 KiB
FortranFixed
SUBROUTINE SLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B,
|
|
$ LDB, WR, WI, X, LDX, SCALE, XNORM, INFO )
|
|
*
|
|
* -- LAPACK auxiliary routine (version 3.1) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
|
|
* November 2006
|
|
*
|
|
* .. Scalar Arguments ..
|
|
LOGICAL LTRANS
|
|
INTEGER INFO, LDA, LDB, LDX, NA, NW
|
|
REAL CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
|
|
* ..
|
|
* .. Array Arguments ..
|
|
REAL A( LDA, * ), B( LDB, * ), X( LDX, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* SLALN2 solves a system of the form (ca A - w D ) X = s B
|
|
* or (ca A' - w D) X = s B with possible scaling ("s") and
|
|
* perturbation of A. (A' 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 (.LE. 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.)
|
|
*
|
|
* Arguments
|
|
* ==========
|
|
*
|
|
* LTRANS (input) LOGICAL
|
|
* =.TRUE.: A-transpose will be used.
|
|
* =.FALSE.: A will be used (not transposed.)
|
|
*
|
|
* NA (input) INTEGER
|
|
* The size of the matrix A. It may (only) be 1 or 2.
|
|
*
|
|
* NW (input) INTEGER
|
|
* 1 if "w" is real, 2 if "w" is complex. It may only be 1
|
|
* or 2.
|
|
*
|
|
* SMIN (input) 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.)
|
|
*
|
|
* CA (input) REAL
|
|
* The coefficient c, which A is multiplied by.
|
|
*
|
|
* A (input) REAL array, dimension (LDA,NA)
|
|
* The NA x NA matrix A.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of A. It must be at least NA.
|
|
*
|
|
* D1 (input) REAL
|
|
* The 1,1 element in the diagonal matrix D.
|
|
*
|
|
* D2 (input) REAL
|
|
* The 2,2 element in the diagonal matrix D. Not used if NW=1.
|
|
*
|
|
* B (input) 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.
|
|
*
|
|
* LDB (input) INTEGER
|
|
* The leading dimension of B. It must be at least NA.
|
|
*
|
|
* WR (input) REAL
|
|
* The real part of the scalar "w".
|
|
*
|
|
* WI (input) REAL
|
|
* The imaginary part of the scalar "w". Not used if NW=1.
|
|
*
|
|
* X (output) 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.
|
|
*
|
|
* LDX (input) INTEGER
|
|
* The leading dimension of X. It must be at least NA.
|
|
*
|
|
* SCALE (output) 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.
|
|
*
|
|
* XNORM (output) REAL
|
|
* The infinity-norm of X, when X is regarded as an NA x NW
|
|
* real matrix.
|
|
*
|
|
* INFO (output) 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.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
REAL ZERO, ONE
|
|
PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
|
|
REAL TWO
|
|
PARAMETER ( TWO = 2.0E0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
INTEGER ICMAX, J
|
|
REAL 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 CSWAP( 4 ), RSWAP( 4 )
|
|
INTEGER IPIVOT( 4, 4 )
|
|
REAL CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )
|
|
* ..
|
|
* .. External Functions ..
|
|
REAL SLAMCH
|
|
EXTERNAL SLAMCH
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL SLADIV
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC ABS, MAX
|
|
* ..
|
|
* .. Equivalences ..
|
|
EQUIVALENCE ( CI( 1, 1 ), CIV( 1 ) ),
|
|
$ ( CR( 1, 1 ), CRV( 1 ) )
|
|
* ..
|
|
* .. Data statements ..
|
|
DATA CSWAP / .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*SLAMCH( '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 SLADIV( 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' - 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( CSWAP( 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 SLADIV( 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( CSWAP( 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 SLALN2
|
|
*
|
|
END
|