LAPACK: SRC/dlaln2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B,
00002      $                   LDB, WR, WI, X, LDX, SCALE, XNORM, INFO )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       LOGICAL            LTRANS
00011       INTEGER            INFO, LDA, LDB, LDX, NA, NW
00012       DOUBLE PRECISION   CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
00013 *     ..
00014 *     .. Array Arguments ..
00015       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), X( LDX, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DLALN2 solves a system of the form  (ca A - w D ) X = s B
00022 *  or (ca A' - w D) X = s B   with possible scaling ("s") and
00023 *  perturbation of A.  (A' means A-transpose.)
00024 *
00025 *  A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
00026 *  real diagonal matrix, w is a real or complex value, and X and B are
00027 *  NA x 1 matrices -- real if w is real, complex if w is complex.  NA
00028 *  may be 1 or 2.
00029 *
00030 *  If w is complex, X and B are represented as NA x 2 matrices,
00031 *  the first column of each being the real part and the second
00032 *  being the imaginary part.
00033 *
00034 *  "s" is a scaling factor (.LE. 1), computed by DLALN2, which is
00035 *  so chosen that X can be computed without overflow.  X is further
00036 *  scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
00037 *  than overflow.
00038 *
00039 *  If both singular values of (ca A - w D) are less than SMIN,
00040 *  SMIN*identity will be used instead of (ca A - w D).  If only one
00041 *  singular value is less than SMIN, one element of (ca A - w D) will be
00042 *  perturbed enough to make the smallest singular value roughly SMIN.
00043 *  If both singular values are at least SMIN, (ca A - w D) will not be
00044 *  perturbed.  In any case, the perturbation will be at most some small
00045 *  multiple of max( SMIN, ulp*norm(ca A - w D) ).  The singular values
00046 *  are computed by infinity-norm approximations, and thus will only be
00047 *  correct to a factor of 2 or so.
00048 *
00049 *  Note: all input quantities are assumed to be smaller than overflow
00050 *  by a reasonable factor.  (See BIGNUM.)
00051 *
00052 *  Arguments
00053 *  ==========
00054 *
00055 *  LTRANS  (input) LOGICAL
00056 *          =.TRUE.:  A-transpose will be used.
00057 *          =.FALSE.: A will be used (not transposed.)
00058 *
00059 *  NA      (input) INTEGER
00060 *          The size of the matrix A.  It may (only) be 1 or 2.
00061 *
00062 *  NW      (input) INTEGER
00063 *          1 if "w" is real, 2 if "w" is complex.  It may only be 1
00064 *          or 2.
00065 *
00066 *  SMIN    (input) DOUBLE PRECISION
00067 *          The desired lower bound on the singular values of A.  This
00068 *          should be a safe distance away from underflow or overflow,
00069 *          say, between (underflow/machine precision) and  (machine
00070 *          precision * overflow ).  (See BIGNUM and ULP.)
00071 *
00072 *  CA      (input) DOUBLE PRECISION
00073 *          The coefficient c, which A is multiplied by.
00074 *
00075 *  A       (input) DOUBLE PRECISION array, dimension (LDA,NA)
00076 *          The NA x NA matrix A.
00077 *
00078 *  LDA     (input) INTEGER
00079 *          The leading dimension of A.  It must be at least NA.
00080 *
00081 *  D1      (input) DOUBLE PRECISION
00082 *          The 1,1 element in the diagonal matrix D.
00083 *
00084 *  D2      (input) DOUBLE PRECISION
00085 *          The 2,2 element in the diagonal matrix D.  Not used if NW=1.
00086 *
00087 *  B       (input) DOUBLE PRECISION array, dimension (LDB,NW)
00088 *          The NA x NW matrix B (right-hand side).  If NW=2 ("w" is
00089 *          complex), column 1 contains the real part of B and column 2
00090 *          contains the imaginary part.
00091 *
00092 *  LDB     (input) INTEGER
00093 *          The leading dimension of B.  It must be at least NA.
00094 *
00095 *  WR      (input) DOUBLE PRECISION
00096 *          The real part of the scalar "w".
00097 *
00098 *  WI      (input) DOUBLE PRECISION
00099 *          The imaginary part of the scalar "w".  Not used if NW=1.
00100 *
00101 *  X       (output) DOUBLE PRECISION array, dimension (LDX,NW)
00102 *          The NA x NW matrix X (unknowns), as computed by DLALN2.
00103 *          If NW=2 ("w" is complex), on exit, column 1 will contain
00104 *          the real part of X and column 2 will contain the imaginary
00105 *          part.
00106 *
00107 *  LDX     (input) INTEGER
00108 *          The leading dimension of X.  It must be at least NA.
00109 *
00110 *  SCALE   (output) DOUBLE PRECISION
00111 *          The scale factor that B must be multiplied by to insure
00112 *          that overflow does not occur when computing X.  Thus,
00113 *          (ca A - w D) X  will be SCALE*B, not B (ignoring
00114 *          perturbations of A.)  It will be at most 1.
00115 *
00116 *  XNORM   (output) DOUBLE PRECISION
00117 *          The infinity-norm of X, when X is regarded as an NA x NW
00118 *          real matrix.
00119 *
00120 *  INFO    (output) INTEGER
00121 *          An error flag.  It will be set to zero if no error occurs,
00122 *          a negative number if an argument is in error, or a positive
00123 *          number if  ca A - w D  had to be perturbed.
00124 *          The possible values are:
00125 *          = 0: No error occurred, and (ca A - w D) did not have to be
00126 *                 perturbed.
00127 *          = 1: (ca A - w D) had to be perturbed to make its smallest
00128 *               (or only) singular value greater than SMIN.
00129 *          NOTE: In the interests of speed, this routine does not
00130 *                check the inputs for errors.
00131 *
00132 * =====================================================================
00133 *
00134 *     .. Parameters ..
00135       DOUBLE PRECISION   ZERO, ONE
00136       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00137       DOUBLE PRECISION   TWO
00138       PARAMETER          ( TWO = 2.0D0 )
00139 *     ..
00140 *     .. Local Scalars ..
00141       INTEGER            ICMAX, J
00142       DOUBLE PRECISION   BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21,
00143      $                   CI22, CMAX, CNORM, CR21, CR22, CSI, CSR, LI21,
00144      $                   LR21, SMINI, SMLNUM, TEMP, U22ABS, UI11, UI11R,
00145      $                   UI12, UI12S, UI22, UR11, UR11R, UR12, UR12S,
00146      $                   UR22, XI1, XI2, XR1, XR2
00147 *     ..
00148 *     .. Local Arrays ..
00149       LOGICAL            RSWAP( 4 ), ZSWAP( 4 )
00150       INTEGER            IPIVOT( 4, 4 )
00151       DOUBLE PRECISION   CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )
00152 *     ..
00153 *     .. External Functions ..
00154       DOUBLE PRECISION   DLAMCH
00155       EXTERNAL           DLAMCH
00156 *     ..
00157 *     .. External Subroutines ..
00158       EXTERNAL           DLADIV
00159 *     ..
00160 *     .. Intrinsic Functions ..
00161       INTRINSIC          ABS, MAX
00162 *     ..
00163 *     .. Equivalences ..
00164       EQUIVALENCE        ( CI( 1, 1 ), CIV( 1 ) ),
00165      $                   ( CR( 1, 1 ), CRV( 1 ) )
00166 *     ..
00167 *     .. Data statements ..
00168       DATA               ZSWAP / .FALSE., .FALSE., .TRUE., .TRUE. /
00169       DATA               RSWAP / .FALSE., .TRUE., .FALSE., .TRUE. /
00170       DATA               IPIVOT / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4,
00171      $                   3, 2, 1 /
00172 *     ..
00173 *     .. Executable Statements ..
00174 *
00175 *     Compute BIGNUM
00176 *
00177       SMLNUM = TWO*DLAMCH( 'Safe minimum' )
00178       BIGNUM = ONE / SMLNUM
00179       SMINI = MAX( SMIN, SMLNUM )
00180 *
00181 *     Don't check for input errors
00182 *
00183       INFO = 0
00184 *
00185 *     Standard Initializations
00186 *
00187       SCALE = ONE
00188 *
00189       IF( NA.EQ.1 ) THEN
00190 *
00191 *        1 x 1  (i.e., scalar) system   C X = B
00192 *
00193          IF( NW.EQ.1 ) THEN
00194 *
00195 *           Real 1x1 system.
00196 *
00197 *           C = ca A - w D
00198 *
00199             CSR = CA*A( 1, 1 ) - WR*D1
00200             CNORM = ABS( CSR )
00201 *
00202 *           If | C | < SMINI, use C = SMINI
00203 *
00204             IF( CNORM.LT.SMINI ) THEN
00205                CSR = SMINI
00206                CNORM = SMINI
00207                INFO = 1
00208             END IF
00209 *
00210 *           Check scaling for  X = B / C
00211 *
00212             BNORM = ABS( B( 1, 1 ) )
00213             IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN
00214                IF( BNORM.GT.BIGNUM*CNORM )
00215      $            SCALE = ONE / BNORM
00216             END IF
00217 *
00218 *           Compute X
00219 *
00220             X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / CSR
00221             XNORM = ABS( X( 1, 1 ) )
00222          ELSE
00223 *
00224 *           Complex 1x1 system (w is complex)
00225 *
00226 *           C = ca A - w D
00227 *
00228             CSR = CA*A( 1, 1 ) - WR*D1
00229             CSI = -WI*D1
00230             CNORM = ABS( CSR ) + ABS( CSI )
00231 *
00232 *           If | C | < SMINI, use C = SMINI
00233 *
00234             IF( CNORM.LT.SMINI ) THEN
00235                CSR = SMINI
00236                CSI = ZERO
00237                CNORM = SMINI
00238                INFO = 1
00239             END IF
00240 *
00241 *           Check scaling for  X = B / C
00242 *
00243             BNORM = ABS( B( 1, 1 ) ) + ABS( B( 1, 2 ) )
00244             IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN
00245                IF( BNORM.GT.BIGNUM*CNORM )
00246      $            SCALE = ONE / BNORM
00247             END IF
00248 *
00249 *           Compute X
00250 *
00251             CALL DLADIV( SCALE*B( 1, 1 ), SCALE*B( 1, 2 ), CSR, CSI,
00252      $                   X( 1, 1 ), X( 1, 2 ) )
00253             XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
00254          END IF
00255 *
00256       ELSE
00257 *
00258 *        2x2 System
00259 *
00260 *        Compute the real part of  C = ca A - w D  (or  ca A' - w D )
00261 *
00262          CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1
00263          CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2
00264          IF( LTRANS ) THEN
00265             CR( 1, 2 ) = CA*A( 2, 1 )
00266             CR( 2, 1 ) = CA*A( 1, 2 )
00267          ELSE
00268             CR( 2, 1 ) = CA*A( 2, 1 )
00269             CR( 1, 2 ) = CA*A( 1, 2 )
00270          END IF
00271 *
00272          IF( NW.EQ.1 ) THEN
00273 *
00274 *           Real 2x2 system  (w is real)
00275 *
00276 *           Find the largest element in C
00277 *
00278             CMAX = ZERO
00279             ICMAX = 0
00280 *
00281             DO 10 J = 1, 4
00282                IF( ABS( CRV( J ) ).GT.CMAX ) THEN
00283                   CMAX = ABS( CRV( J ) )
00284                   ICMAX = J
00285                END IF
00286    10       CONTINUE
00287 *
00288 *           If norm(C) < SMINI, use SMINI*identity.
00289 *
00290             IF( CMAX.LT.SMINI ) THEN
00291                BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 2, 1 ) ) )
00292                IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN
00293                   IF( BNORM.GT.BIGNUM*SMINI )
00294      $               SCALE = ONE / BNORM
00295                END IF
00296                TEMP = SCALE / SMINI
00297                X( 1, 1 ) = TEMP*B( 1, 1 )
00298                X( 2, 1 ) = TEMP*B( 2, 1 )
00299                XNORM = TEMP*BNORM
00300                INFO = 1
00301                RETURN
00302             END IF
00303 *
00304 *           Gaussian elimination with complete pivoting.
00305 *
00306             UR11 = CRV( ICMAX )
00307             CR21 = CRV( IPIVOT( 2, ICMAX ) )
00308             UR12 = CRV( IPIVOT( 3, ICMAX ) )
00309             CR22 = CRV( IPIVOT( 4, ICMAX ) )
00310             UR11R = ONE / UR11
00311             LR21 = UR11R*CR21
00312             UR22 = CR22 - UR12*LR21
00313 *
00314 *           If smaller pivot < SMINI, use SMINI
00315 *
00316             IF( ABS( UR22 ).LT.SMINI ) THEN
00317                UR22 = SMINI
00318                INFO = 1
00319             END IF
00320             IF( RSWAP( ICMAX ) ) THEN
00321                BR1 = B( 2, 1 )
00322                BR2 = B( 1, 1 )
00323             ELSE
00324                BR1 = B( 1, 1 )
00325                BR2 = B( 2, 1 )
00326             END IF
00327             BR2 = BR2 - LR21*BR1
00328             BBND = MAX( ABS( BR1*( UR22*UR11R ) ), ABS( BR2 ) )
00329             IF( BBND.GT.ONE .AND. ABS( UR22 ).LT.ONE ) THEN
00330                IF( BBND.GE.BIGNUM*ABS( UR22 ) )
00331      $            SCALE = ONE / BBND
00332             END IF
00333 *
00334             XR2 = ( BR2*SCALE ) / UR22
00335             XR1 = ( SCALE*BR1 )*UR11R - XR2*( UR11R*UR12 )
00336             IF( ZSWAP( ICMAX ) ) THEN
00337                X( 1, 1 ) = XR2
00338                X( 2, 1 ) = XR1
00339             ELSE
00340                X( 1, 1 ) = XR1
00341                X( 2, 1 ) = XR2
00342             END IF
00343             XNORM = MAX( ABS( XR1 ), ABS( XR2 ) )
00344 *
00345 *           Further scaling if  norm(A) norm(X) > overflow
00346 *
00347             IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN
00348                IF( XNORM.GT.BIGNUM / CMAX ) THEN
00349                   TEMP = CMAX / BIGNUM
00350                   X( 1, 1 ) = TEMP*X( 1, 1 )
00351                   X( 2, 1 ) = TEMP*X( 2, 1 )
00352                   XNORM = TEMP*XNORM
00353                   SCALE = TEMP*SCALE
00354                END IF
00355             END IF
00356          ELSE
00357 *
00358 *           Complex 2x2 system  (w is complex)
00359 *
00360 *           Find the largest element in C
00361 *
00362             CI( 1, 1 ) = -WI*D1
00363             CI( 2, 1 ) = ZERO
00364             CI( 1, 2 ) = ZERO
00365             CI( 2, 2 ) = -WI*D2
00366             CMAX = ZERO
00367             ICMAX = 0
00368 *
00369             DO 20 J = 1, 4
00370                IF( ABS( CRV( J ) )+ABS( CIV( J ) ).GT.CMAX ) THEN
00371                   CMAX = ABS( CRV( J ) ) + ABS( CIV( J ) )
00372                   ICMAX = J
00373                END IF
00374    20       CONTINUE
00375 *
00376 *           If norm(C) < SMINI, use SMINI*identity.
00377 *
00378             IF( CMAX.LT.SMINI ) THEN
00379                BNORM = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ),
00380      $                 ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) )
00381                IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN
00382                   IF( BNORM.GT.BIGNUM*SMINI )
00383      $               SCALE = ONE / BNORM
00384                END IF
00385                TEMP = SCALE / SMINI
00386                X( 1, 1 ) = TEMP*B( 1, 1 )
00387                X( 2, 1 ) = TEMP*B( 2, 1 )
00388                X( 1, 2 ) = TEMP*B( 1, 2 )
00389                X( 2, 2 ) = TEMP*B( 2, 2 )
00390                XNORM = TEMP*BNORM
00391                INFO = 1
00392                RETURN
00393             END IF
00394 *
00395 *           Gaussian elimination with complete pivoting.
00396 *
00397             UR11 = CRV( ICMAX )
00398             UI11 = CIV( ICMAX )
00399             CR21 = CRV( IPIVOT( 2, ICMAX ) )
00400             CI21 = CIV( IPIVOT( 2, ICMAX ) )
00401             UR12 = CRV( IPIVOT( 3, ICMAX ) )
00402             UI12 = CIV( IPIVOT( 3, ICMAX ) )
00403             CR22 = CRV( IPIVOT( 4, ICMAX ) )
00404             CI22 = CIV( IPIVOT( 4, ICMAX ) )
00405             IF( ICMAX.EQ.1 .OR. ICMAX.EQ.4 ) THEN
00406 *
00407 *              Code when off-diagonals of pivoted C are real
00408 *
00409                IF( ABS( UR11 ).GT.ABS( UI11 ) ) THEN
00410                   TEMP = UI11 / UR11
00411                   UR11R = ONE / ( UR11*( ONE+TEMP**2 ) )
00412                   UI11R = -TEMP*UR11R
00413                ELSE
00414                   TEMP = UR11 / UI11
00415                   UI11R = -ONE / ( UI11*( ONE+TEMP**2 ) )
00416                   UR11R = -TEMP*UI11R
00417                END IF
00418                LR21 = CR21*UR11R
00419                LI21 = CR21*UI11R
00420                UR12S = UR12*UR11R
00421                UI12S = UR12*UI11R
00422                UR22 = CR22 - UR12*LR21
00423                UI22 = CI22 - UR12*LI21
00424             ELSE
00425 *
00426 *              Code when diagonals of pivoted C are real
00427 *
00428                UR11R = ONE / UR11
00429                UI11R = ZERO
00430                LR21 = CR21*UR11R
00431                LI21 = CI21*UR11R
00432                UR12S = UR12*UR11R
00433                UI12S = UI12*UR11R
00434                UR22 = CR22 - UR12*LR21 + UI12*LI21
00435                UI22 = -UR12*LI21 - UI12*LR21
00436             END IF
00437             U22ABS = ABS( UR22 ) + ABS( UI22 )
00438 *
00439 *           If smaller pivot < SMINI, use SMINI
00440 *
00441             IF( U22ABS.LT.SMINI ) THEN
00442                UR22 = SMINI
00443                UI22 = ZERO
00444                INFO = 1
00445             END IF
00446             IF( RSWAP( ICMAX ) ) THEN
00447                BR2 = B( 1, 1 )
00448                BR1 = B( 2, 1 )
00449                BI2 = B( 1, 2 )
00450                BI1 = B( 2, 2 )
00451             ELSE
00452                BR1 = B( 1, 1 )
00453                BR2 = B( 2, 1 )
00454                BI1 = B( 1, 2 )
00455                BI2 = B( 2, 2 )
00456             END IF
00457             BR2 = BR2 - LR21*BR1 + LI21*BI1
00458             BI2 = BI2 - LI21*BR1 - LR21*BI1
00459             BBND = MAX( ( ABS( BR1 )+ABS( BI1 ) )*
00460      $             ( U22ABS*( ABS( UR11R )+ABS( UI11R ) ) ),
00461      $             ABS( BR2 )+ABS( BI2 ) )
00462             IF( BBND.GT.ONE .AND. U22ABS.LT.ONE ) THEN
00463                IF( BBND.GE.BIGNUM*U22ABS ) THEN
00464                   SCALE = ONE / BBND
00465                   BR1 = SCALE*BR1
00466                   BI1 = SCALE*BI1
00467                   BR2 = SCALE*BR2
00468                   BI2 = SCALE*BI2
00469                END IF
00470             END IF
00471 *
00472             CALL DLADIV( BR2, BI2, UR22, UI22, XR2, XI2 )
00473             XR1 = UR11R*BR1 - UI11R*BI1 - UR12S*XR2 + UI12S*XI2
00474             XI1 = UI11R*BR1 + UR11R*BI1 - UI12S*XR2 - UR12S*XI2
00475             IF( ZSWAP( ICMAX ) ) THEN
00476                X( 1, 1 ) = XR2
00477                X( 2, 1 ) = XR1
00478                X( 1, 2 ) = XI2
00479                X( 2, 2 ) = XI1
00480             ELSE
00481                X( 1, 1 ) = XR1
00482                X( 2, 1 ) = XR2
00483                X( 1, 2 ) = XI1
00484                X( 2, 2 ) = XI2
00485             END IF
00486             XNORM = MAX( ABS( XR1 )+ABS( XI1 ), ABS( XR2 )+ABS( XI2 ) )
00487 *
00488 *           Further scaling if  norm(A) norm(X) > overflow
00489 *
00490             IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN
00491                IF( XNORM.GT.BIGNUM / CMAX ) THEN
00492                   TEMP = CMAX / BIGNUM
00493                   X( 1, 1 ) = TEMP*X( 1, 1 )
00494                   X( 2, 1 ) = TEMP*X( 2, 1 )
00495                   X( 1, 2 ) = TEMP*X( 1, 2 )
00496                   X( 2, 2 ) = TEMP*X( 2, 2 )
00497                   XNORM = TEMP*XNORM
00498                   SCALE = TEMP*SCALE
00499                END IF
00500             END IF
00501          END IF
00502       END IF
00503 *
00504       RETURN
00505 *
00506 *     End of DLALN2
00507 *
00508       END