#include "ESMF_LapackBlas.inc" !> \brief \b DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form. ! ! =========== DOCUMENTATION =========== ! ! Online html documentation available at ! http://www.netlib.org/lapack/explore-html/ ! !> \htmlonly !> Download DLALN2 + dependencies !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaln2.f"> !> [TGZ]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaln2.f"> !> [ZIP]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaln2.f"> !> [TXT]</a> !> \endhtmlonly ! ! Definition: ! =========== ! ! SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, ! LDB, WR, WI, X, LDX, SCALE, XNORM, INFO ) ! ! .. 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, * ) ! .. ! ! !> \par Purpose: ! ============= !> !> \verbatim !> !> 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 ("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 (.LE. 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.) !> \endverbatim ! ! Arguments: ! ========== ! !> \param[in] LTRANS !> \verbatim !> LTRANS is LOGICAL !> =.TRUE.: A-transpose will be used. !> =.FALSE.: A will be used (not transposed.) !> \endverbatim !> !> \param[in] NA !> \verbatim !> NA is INTEGER !> The size of the matrix A. It may (only) be 1 or 2. !> \endverbatim !> !> \param[in] NW !> \verbatim !> NW is INTEGER !> 1 if "w" is real, 2 if "w" is complex. It may only be 1 !> or 2. !> \endverbatim !> !> \param[in] SMIN !> \verbatim !> 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.) !> \endverbatim !> !> \param[in] CA !> \verbatim !> CA is DOUBLE PRECISION !> The coefficient c, which A is multiplied by. !> \endverbatim !> !> \param[in] A !> \verbatim !> A is DOUBLE PRECISION array, dimension (LDA,NA) !> The NA x NA matrix A. !> \endverbatim !> !> \param[in] LDA !> \verbatim !> LDA is INTEGER !> The leading dimension of A. It must be at least NA. !> \endverbatim !> !> \param[in] D1 !> \verbatim !> D1 is DOUBLE PRECISION !> The 1,1 element in the diagonal matrix D. !> \endverbatim !> !> \param[in] D2 !> \verbatim !> D2 is DOUBLE PRECISION !> The 2,2 element in the diagonal matrix D. Not used if NA=1. !> \endverbatim !> !> \param[in] B !> \verbatim !> B is DOUBLE PRECISION 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. !> \endverbatim !> !> \param[in] LDB !> \verbatim !> LDB is INTEGER !> The leading dimension of B. It must be at least NA. !> \endverbatim !> !> \param[in] WR !> \verbatim !> WR is DOUBLE PRECISION !> The real part of the scalar "w". !> \endverbatim !> !> \param[in] WI !> \verbatim !> WI is DOUBLE PRECISION !> The imaginary part of the scalar "w". Not used if NW=1. !> \endverbatim !> !> \param[out] X !> \verbatim !> X is DOUBLE PRECISION array, dimension (LDX,NW) !> The NA x NW matrix X (unknowns), as computed by DLALN2. !> 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. !> \endverbatim !> !> \param[in] LDX !> \verbatim !> LDX is INTEGER !> The leading dimension of X. It must be at least NA. !> \endverbatim !> !> \param[out] SCALE !> \verbatim !> 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. !> \endverbatim !> !> \param[out] XNORM !> \verbatim !> XNORM is DOUBLE PRECISION !> The infinity-norm of X, when X is regarded as an NA x NW !> real matrix. !> \endverbatim !> !> \param[out] INFO !> \verbatim !> 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. !> \endverbatim ! ! Authors: ! ======== ! !> \author Univ. of Tennessee !> \author Univ. of California Berkeley !> \author Univ. of Colorado Denver !> \author NAG Ltd. ! !> \date December 2016 ! !> \ingroup doubleOTHERauxiliary ! ! ===================================================================== SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, & LDB, WR, WI, X, LDX, SCALE, XNORM, INFO ) ! ! -- LAPACK auxiliary routine (version 3.7.0) -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! December 2016 ! ! .. 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 ! END