Numérotation des lignes :

C DLALN2    SOURCE    BP208322  18/07/10    21:15:09     9872           *> \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*> &lt;a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaln2.f">*> [TGZ]&lt;/a>*> &lt;a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaln2.f">*> [ZIP]&lt;/a>*> &lt;a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaln2.f">*> [TXT]&lt;/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*       REAL*8   CA, D1, D2, SCALE, SMIN, WI, WR, XNORM*       ..*       .. Array Arguments ..*       REAL*8   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 REAL*8 CA, D1, D2, SCALE, SMIN, WI, WR, XNORM* ..* .. Array Arguments .. REAL*8 A( LDA, * ), B( LDB, * ), X( LDX, * )* ..** =====================================================================** .. Parameters .. REAL*8 ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) REAL*8 TWO PARAMETER ( TWO = 2.0D0 )* ..* .. Local Scalars .. INTEGER ICMAX, J REAL*8 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 ) REAL*8 CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )* ..* .. External Functions .. REAL*8 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 | &lt; 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 | &lt; 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) &lt; 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 &lt; 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) &lt; 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 &lt; 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

© Cast3M 2003 - Tous droits réservés.
Mentions légales