diff options
Diffstat (limited to 'src/lib/lapack/dlaln2.f')
| -rw-r--r-- | src/lib/lapack/dlaln2.f | 507 |
1 files changed, 0 insertions, 507 deletions
diff --git a/src/lib/lapack/dlaln2.f b/src/lib/lapack/dlaln2.f deleted file mode 100644 index 7c99bdbe..00000000 --- a/src/lib/lapack/dlaln2.f +++ /dev/null @@ -1,507 +0,0 @@ - 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.1) -- -* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. -* November 2006 -* -* .. 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, * ) -* .. -* -* Purpose -* ======= -* -* DLALN2 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 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.) -* -* 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) 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.) -* -* CA (input) DOUBLE PRECISION -* The coefficient c, which A is multiplied by. -* -* A (input) DOUBLE PRECISION 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) DOUBLE PRECISION -* The 1,1 element in the diagonal matrix D. -* -* D2 (input) DOUBLE PRECISION -* The 2,2 element in the diagonal matrix D. Not used if NW=1. -* -* B (input) 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. -* -* LDB (input) INTEGER -* The leading dimension of B. It must be at least NA. -* -* WR (input) DOUBLE PRECISION -* The real part of the scalar "w". -* -* WI (input) DOUBLE PRECISION -* The imaginary part of the scalar "w". Not used if NW=1. -* -* X (output) 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. -* -* LDX (input) INTEGER -* The leading dimension of X. It must be at least NA. -* -* SCALE (output) 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. -* -* XNORM (output) DOUBLE PRECISION -* 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 .. - 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' - 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 |