LAPACK: SRC/dlaqtr.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAQTR( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK,
00002      $                   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            LREAL, LTRAN
00011       INTEGER            INFO, LDT, N
00012       DOUBLE PRECISION   SCALE, W
00013 *     ..
00014 *     .. Array Arguments ..
00015       DOUBLE PRECISION   B( * ), T( LDT, * ), WORK( * ), X( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DLAQTR solves the real quasi-triangular system
00022 *
00023 *               op(T)*p = scale*c,               if LREAL = .TRUE.
00024 *
00025 *  or the complex quasi-triangular systems
00026 *
00027 *             op(T + iB)*(p+iq) = scale*(c+id),  if LREAL = .FALSE.
00028 *
00029 *  in real arithmetic, where T is upper quasi-triangular.
00030 *  If LREAL = .FALSE., then the first diagonal block of T must be
00031 *  1 by 1, B is the specially structured matrix
00032 *
00033 *                 B = [ b(1) b(2) ... b(n) ]
00034 *                     [       w            ]
00035 *                     [           w        ]
00036 *                     [              .     ]
00037 *                     [                 w  ]
00038 *
00039 *  op(A) = A or A', A' denotes the conjugate transpose of
00040 *  matrix A.
00041 *
00042 *  On input, X = [ c ].  On output, X = [ p ].
00043 *                [ d ]                  [ q ]
00044 *
00045 *  This subroutine is designed for the condition number estimation
00046 *  in routine DTRSNA.
00047 *
00048 *  Arguments
00049 *  =========
00050 *
00051 *  LTRAN   (input) LOGICAL
00052 *          On entry, LTRAN specifies the option of conjugate transpose:
00053 *             = .FALSE.,    op(T+i*B) = T+i*B,
00054 *             = .TRUE.,     op(T+i*B) = (T+i*B)'.
00055 *
00056 *  LREAL   (input) LOGICAL
00057 *          On entry, LREAL specifies the input matrix structure:
00058 *             = .FALSE.,    the input is complex
00059 *             = .TRUE.,     the input is real
00060 *
00061 *  N       (input) INTEGER
00062 *          On entry, N specifies the order of T+i*B. N >= 0.
00063 *
00064 *  T       (input) DOUBLE PRECISION array, dimension (LDT,N)
00065 *          On entry, T contains a matrix in Schur canonical form.
00066 *          If LREAL = .FALSE., then the first diagonal block of T mu
00067 *          be 1 by 1.
00068 *
00069 *  LDT     (input) INTEGER
00070 *          The leading dimension of the matrix T. LDT >= max(1,N).
00071 *
00072 *  B       (input) DOUBLE PRECISION array, dimension (N)
00073 *          On entry, B contains the elements to form the matrix
00074 *          B as described above.
00075 *          If LREAL = .TRUE., B is not referenced.
00076 *
00077 *  W       (input) DOUBLE PRECISION
00078 *          On entry, W is the diagonal element of the matrix B.
00079 *          If LREAL = .TRUE., W is not referenced.
00080 *
00081 *  SCALE   (output) DOUBLE PRECISION
00082 *          On exit, SCALE is the scale factor.
00083 *
00084 *  X       (input/output) DOUBLE PRECISION array, dimension (2*N)
00085 *          On entry, X contains the right hand side of the system.
00086 *          On exit, X is overwritten by the solution.
00087 *
00088 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
00089 *
00090 *  INFO    (output) INTEGER
00091 *          On exit, INFO is set to
00092 *             0: successful exit.
00093 *               1: the some diagonal 1 by 1 block has been perturbed by
00094 *                  a small number SMIN to keep nonsingularity.
00095 *               2: the some diagonal 2 by 2 block has been perturbed by
00096 *                  a small number in DLALN2 to keep nonsingularity.
00097 *          NOTE: In the interests of speed, this routine does not
00098 *                check the inputs for errors.
00099 *
00100 * =====================================================================
00101 *
00102 *     .. Parameters ..
00103       DOUBLE PRECISION   ZERO, ONE
00104       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00105 *     ..
00106 *     .. Local Scalars ..
00107       LOGICAL            NOTRAN
00108       INTEGER            I, IERR, J, J1, J2, JNEXT, K, N1, N2
00109       DOUBLE PRECISION   BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW,
00110      $                   SMLNUM, SR, TJJ, TMP, XJ, XMAX, XNORM, Z
00111 *     ..
00112 *     .. Local Arrays ..
00113       DOUBLE PRECISION   D( 2, 2 ), V( 2, 2 )
00114 *     ..
00115 *     .. External Functions ..
00116       INTEGER            IDAMAX
00117       DOUBLE PRECISION   DASUM, DDOT, DLAMCH, DLANGE
00118       EXTERNAL           IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
00119 *     ..
00120 *     .. External Subroutines ..
00121       EXTERNAL           DAXPY, DLADIV, DLALN2, DSCAL
00122 *     ..
00123 *     .. Intrinsic Functions ..
00124       INTRINSIC          ABS, MAX
00125 *     ..
00126 *     .. Executable Statements ..
00127 *
00128 *     Do not test the input parameters for errors
00129 *
00130       NOTRAN = .NOT.LTRAN
00131       INFO = 0
00132 *
00133 *     Quick return if possible
00134 *
00135       IF( N.EQ.0 )
00136      $   RETURN
00137 *
00138 *     Set constants to control overflow
00139 *
00140       EPS = DLAMCH( 'P' )
00141       SMLNUM = DLAMCH( 'S' ) / EPS
00142       BIGNUM = ONE / SMLNUM
00143 *
00144       XNORM = DLANGE( 'M', N, N, T, LDT, D )
00145       IF( .NOT.LREAL )
00146      $   XNORM = MAX( XNORM, ABS( W ), DLANGE( 'M', N, 1, B, N, D ) )
00147       SMIN = MAX( SMLNUM, EPS*XNORM )
00148 *
00149 *     Compute 1-norm of each column of strictly upper triangular
00150 *     part of T to control overflow in triangular solver.
00151 *
00152       WORK( 1 ) = ZERO
00153       DO 10 J = 2, N
00154          WORK( J ) = DASUM( J-1, T( 1, J ), 1 )
00155    10 CONTINUE
00156 *
00157       IF( .NOT.LREAL ) THEN
00158          DO 20 I = 2, N
00159             WORK( I ) = WORK( I ) + ABS( B( I ) )
00160    20    CONTINUE
00161       END IF
00162 *
00163       N2 = 2*N
00164       N1 = N
00165       IF( .NOT.LREAL )
00166      $   N1 = N2
00167       K = IDAMAX( N1, X, 1 )
00168       XMAX = ABS( X( K ) )
00169       SCALE = ONE
00170 *
00171       IF( XMAX.GT.BIGNUM ) THEN
00172          SCALE = BIGNUM / XMAX
00173          CALL DSCAL( N1, SCALE, X, 1 )
00174          XMAX = BIGNUM
00175       END IF
00176 *
00177       IF( LREAL ) THEN
00178 *
00179          IF( NOTRAN ) THEN
00180 *
00181 *           Solve T*p = scale*c
00182 *
00183             JNEXT = N
00184             DO 30 J = N, 1, -1
00185                IF( J.GT.JNEXT )
00186      $            GO TO 30
00187                J1 = J
00188                J2 = J
00189                JNEXT = J - 1
00190                IF( J.GT.1 ) THEN
00191                   IF( T( J, J-1 ).NE.ZERO ) THEN
00192                      J1 = J - 1
00193                      JNEXT = J - 2
00194                   END IF
00195                END IF
00196 *
00197                IF( J1.EQ.J2 ) THEN
00198 *
00199 *                 Meet 1 by 1 diagonal block
00200 *
00201 *                 Scale to avoid overflow when computing
00202 *                     x(j) = b(j)/T(j,j)
00203 *
00204                   XJ = ABS( X( J1 ) )
00205                   TJJ = ABS( T( J1, J1 ) )
00206                   TMP = T( J1, J1 )
00207                   IF( TJJ.LT.SMIN ) THEN
00208                      TMP = SMIN
00209                      TJJ = SMIN
00210                      INFO = 1
00211                   END IF
00212 *
00213                   IF( XJ.EQ.ZERO )
00214      $               GO TO 30
00215 *
00216                   IF( TJJ.LT.ONE ) THEN
00217                      IF( XJ.GT.BIGNUM*TJJ ) THEN
00218                         REC = ONE / XJ
00219                         CALL DSCAL( N, REC, X, 1 )
00220                         SCALE = SCALE*REC
00221                         XMAX = XMAX*REC
00222                      END IF
00223                   END IF
00224                   X( J1 ) = X( J1 ) / TMP
00225                   XJ = ABS( X( J1 ) )
00226 *
00227 *                 Scale x if necessary to avoid overflow when adding a
00228 *                 multiple of column j1 of T.
00229 *
00230                   IF( XJ.GT.ONE ) THEN
00231                      REC = ONE / XJ
00232                      IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
00233                         CALL DSCAL( N, REC, X, 1 )
00234                         SCALE = SCALE*REC
00235                      END IF
00236                   END IF
00237                   IF( J1.GT.1 ) THEN
00238                      CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
00239                      K = IDAMAX( J1-1, X, 1 )
00240                      XMAX = ABS( X( K ) )
00241                   END IF
00242 *
00243                ELSE
00244 *
00245 *                 Meet 2 by 2 diagonal block
00246 *
00247 *                 Call 2 by 2 linear system solve, to take
00248 *                 care of possible overflow by scaling factor.
00249 *
00250                   D( 1, 1 ) = X( J1 )
00251                   D( 2, 1 ) = X( J2 )
00252                   CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, T( J1, J1 ),
00253      $                         LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
00254      $                         SCALOC, XNORM, IERR )
00255                   IF( IERR.NE.0 )
00256      $               INFO = 2
00257 *
00258                   IF( SCALOC.NE.ONE ) THEN
00259                      CALL DSCAL( N, SCALOC, X, 1 )
00260                      SCALE = SCALE*SCALOC
00261                   END IF
00262                   X( J1 ) = V( 1, 1 )
00263                   X( J2 ) = V( 2, 1 )
00264 *
00265 *                 Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2))
00266 *                 to avoid overflow in updating right-hand side.
00267 *
00268                   XJ = MAX( ABS( V( 1, 1 ) ), ABS( V( 2, 1 ) ) )
00269                   IF( XJ.GT.ONE ) THEN
00270                      REC = ONE / XJ
00271                      IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
00272      $                   ( BIGNUM-XMAX )*REC ) THEN
00273                         CALL DSCAL( N, REC, X, 1 )
00274                         SCALE = SCALE*REC
00275                      END IF
00276                   END IF
00277 *
00278 *                 Update right-hand side
00279 *
00280                   IF( J1.GT.1 ) THEN
00281                      CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
00282                      CALL DAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
00283                      K = IDAMAX( J1-1, X, 1 )
00284                      XMAX = ABS( X( K ) )
00285                   END IF
00286 *
00287                END IF
00288 *
00289    30       CONTINUE
00290 *
00291          ELSE
00292 *
00293 *           Solve T'*p = scale*c
00294 *
00295             JNEXT = 1
00296             DO 40 J = 1, N
00297                IF( J.LT.JNEXT )
00298      $            GO TO 40
00299                J1 = J
00300                J2 = J
00301                JNEXT = J + 1
00302                IF( J.LT.N ) THEN
00303                   IF( T( J+1, J ).NE.ZERO ) THEN
00304                      J2 = J + 1
00305                      JNEXT = J + 2
00306                   END IF
00307                END IF
00308 *
00309                IF( J1.EQ.J2 ) THEN
00310 *
00311 *                 1 by 1 diagonal block
00312 *
00313 *                 Scale if necessary to avoid overflow in forming the
00314 *                 right-hand side element by inner product.
00315 *
00316                   XJ = ABS( X( J1 ) )
00317                   IF( XMAX.GT.ONE ) THEN
00318                      REC = ONE / XMAX
00319                      IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
00320                         CALL DSCAL( N, REC, X, 1 )
00321                         SCALE = SCALE*REC
00322                         XMAX = XMAX*REC
00323                      END IF
00324                   END IF
00325 *
00326                   X( J1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 1 )
00327 *
00328                   XJ = ABS( X( J1 ) )
00329                   TJJ = ABS( T( J1, J1 ) )
00330                   TMP = T( J1, J1 )
00331                   IF( TJJ.LT.SMIN ) THEN
00332                      TMP = SMIN
00333                      TJJ = SMIN
00334                      INFO = 1
00335                   END IF
00336 *
00337                   IF( TJJ.LT.ONE ) THEN
00338                      IF( XJ.GT.BIGNUM*TJJ ) THEN
00339                         REC = ONE / XJ
00340                         CALL DSCAL( N, REC, X, 1 )
00341                         SCALE = SCALE*REC
00342                         XMAX = XMAX*REC
00343                      END IF
00344                   END IF
00345                   X( J1 ) = X( J1 ) / TMP
00346                   XMAX = MAX( XMAX, ABS( X( J1 ) ) )
00347 *
00348                ELSE
00349 *
00350 *                 2 by 2 diagonal block
00351 *
00352 *                 Scale if necessary to avoid overflow in forming the
00353 *                 right-hand side elements by inner product.
00354 *
00355                   XJ = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ) )
00356                   IF( XMAX.GT.ONE ) THEN
00357                      REC = ONE / XMAX
00358                      IF( MAX( WORK( J2 ), WORK( J1 ) ).GT.( BIGNUM-XJ )*
00359      $                   REC ) THEN
00360                         CALL DSCAL( N, REC, X, 1 )
00361                         SCALE = SCALE*REC
00362                         XMAX = XMAX*REC
00363                      END IF
00364                   END IF
00365 *
00366                   D( 1, 1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X,
00367      $                        1 )
00368                   D( 2, 1 ) = X( J2 ) - DDOT( J1-1, T( 1, J2 ), 1, X,
00369      $                        1 )
00370 *
00371                   CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J1, J1 ),
00372      $                         LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
00373      $                         SCALOC, XNORM, IERR )
00374                   IF( IERR.NE.0 )
00375      $               INFO = 2
00376 *
00377                   IF( SCALOC.NE.ONE ) THEN
00378                      CALL DSCAL( N, SCALOC, X, 1 )
00379                      SCALE = SCALE*SCALOC
00380                   END IF
00381                   X( J1 ) = V( 1, 1 )
00382                   X( J2 ) = V( 2, 1 )
00383                   XMAX = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ), XMAX )
00384 *
00385                END IF
00386    40       CONTINUE
00387          END IF
00388 *
00389       ELSE
00390 *
00391          SMINW = MAX( EPS*ABS( W ), SMIN )
00392          IF( NOTRAN ) THEN
00393 *
00394 *           Solve (T + iB)*(p+iq) = c+id
00395 *
00396             JNEXT = N
00397             DO 70 J = N, 1, -1
00398                IF( J.GT.JNEXT )
00399      $            GO TO 70
00400                J1 = J
00401                J2 = J
00402                JNEXT = J - 1
00403                IF( J.GT.1 ) THEN
00404                   IF( T( J, J-1 ).NE.ZERO ) THEN
00405                      J1 = J - 1
00406                      JNEXT = J - 2
00407                   END IF
00408                END IF
00409 *
00410                IF( J1.EQ.J2 ) THEN
00411 *
00412 *                 1 by 1 diagonal block
00413 *
00414 *                 Scale if necessary to avoid overflow in division
00415 *
00416                   Z = W
00417                   IF( J1.EQ.1 )
00418      $               Z = B( 1 )
00419                   XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) )
00420                   TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
00421                   TMP = T( J1, J1 )
00422                   IF( TJJ.LT.SMINW ) THEN
00423                      TMP = SMINW
00424                      TJJ = SMINW
00425                      INFO = 1
00426                   END IF
00427 *
00428                   IF( XJ.EQ.ZERO )
00429      $               GO TO 70
00430 *
00431                   IF( TJJ.LT.ONE ) THEN
00432                      IF( XJ.GT.BIGNUM*TJJ ) THEN
00433                         REC = ONE / XJ
00434                         CALL DSCAL( N2, REC, X, 1 )
00435                         SCALE = SCALE*REC
00436                         XMAX = XMAX*REC
00437                      END IF
00438                   END IF
00439                   CALL DLADIV( X( J1 ), X( N+J1 ), TMP, Z, SR, SI )
00440                   X( J1 ) = SR
00441                   X( N+J1 ) = SI
00442                   XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) )
00443 *
00444 *                 Scale x if necessary to avoid overflow when adding a
00445 *                 multiple of column j1 of T.
00446 *
00447                   IF( XJ.GT.ONE ) THEN
00448                      REC = ONE / XJ
00449                      IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
00450                         CALL DSCAL( N2, REC, X, 1 )
00451                         SCALE = SCALE*REC
00452                      END IF
00453                   END IF
00454 *
00455                   IF( J1.GT.1 ) THEN
00456                      CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
00457                      CALL DAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
00458      $                           X( N+1 ), 1 )
00459 *
00460                      X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 )
00461                      X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 )
00462 *
00463                      XMAX = ZERO
00464                      DO 50 K = 1, J1 - 1
00465                         XMAX = MAX( XMAX, ABS( X( K ) )+
00466      $                         ABS( X( K+N ) ) )
00467    50                CONTINUE
00468                   END IF
00469 *
00470                ELSE
00471 *
00472 *                 Meet 2 by 2 diagonal block
00473 *
00474                   D( 1, 1 ) = X( J1 )
00475                   D( 2, 1 ) = X( J2 )
00476                   D( 1, 2 ) = X( N+J1 )
00477                   D( 2, 2 ) = X( N+J2 )
00478                   CALL DLALN2( .FALSE., 2, 2, SMINW, ONE, T( J1, J1 ),
00479      $                         LDT, ONE, ONE, D, 2, ZERO, -W, V, 2,
00480      $                         SCALOC, XNORM, IERR )
00481                   IF( IERR.NE.0 )
00482      $               INFO = 2
00483 *
00484                   IF( SCALOC.NE.ONE ) THEN
00485                      CALL DSCAL( 2*N, SCALOC, X, 1 )
00486                      SCALE = SCALOC*SCALE
00487                   END IF
00488                   X( J1 ) = V( 1, 1 )
00489                   X( J2 ) = V( 2, 1 )
00490                   X( N+J1 ) = V( 1, 2 )
00491                   X( N+J2 ) = V( 2, 2 )
00492 *
00493 *                 Scale X(J1), .... to avoid overflow in
00494 *                 updating right hand side.
00495 *
00496                   XJ = MAX( ABS( V( 1, 1 ) )+ABS( V( 1, 2 ) ),
00497      $                 ABS( V( 2, 1 ) )+ABS( V( 2, 2 ) ) )
00498                   IF( XJ.GT.ONE ) THEN
00499                      REC = ONE / XJ
00500                      IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
00501      $                   ( BIGNUM-XMAX )*REC ) THEN
00502                         CALL DSCAL( N2, REC, X, 1 )
00503                         SCALE = SCALE*REC
00504                      END IF
00505                   END IF
00506 *
00507 *                 Update the right-hand side.
00508 *
00509                   IF( J1.GT.1 ) THEN
00510                      CALL DAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
00511                      CALL DAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
00512 *
00513                      CALL DAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
00514      $                           X( N+1 ), 1 )
00515                      CALL DAXPY( J1-1, -X( N+J2 ), T( 1, J2 ), 1,
00516      $                           X( N+1 ), 1 )
00517 *
00518                      X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 ) +
00519      $                        B( J2 )*X( N+J2 )
00520                      X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 ) -
00521      $                          B( J2 )*X( J2 )
00522 *
00523                      XMAX = ZERO
00524                      DO 60 K = 1, J1 - 1
00525                         XMAX = MAX( ABS( X( K ) )+ABS( X( K+N ) ),
00526      $                         XMAX )
00527    60                CONTINUE
00528                   END IF
00529 *
00530                END IF
00531    70       CONTINUE
00532 *
00533          ELSE
00534 *
00535 *           Solve (T + iB)'*(p+iq) = c+id
00536 *
00537             JNEXT = 1
00538             DO 80 J = 1, N
00539                IF( J.LT.JNEXT )
00540      $            GO TO 80
00541                J1 = J
00542                J2 = J
00543                JNEXT = J + 1
00544                IF( J.LT.N ) THEN
00545                   IF( T( J+1, J ).NE.ZERO ) THEN
00546                      J2 = J + 1
00547                      JNEXT = J + 2
00548                   END IF
00549                END IF
00550 *
00551                IF( J1.EQ.J2 ) THEN
00552 *
00553 *                 1 by 1 diagonal block
00554 *
00555 *                 Scale if necessary to avoid overflow in forming the
00556 *                 right-hand side element by inner product.
00557 *
00558                   XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
00559                   IF( XMAX.GT.ONE ) THEN
00560                      REC = ONE / XMAX
00561                      IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
00562                         CALL DSCAL( N2, REC, X, 1 )
00563                         SCALE = SCALE*REC
00564                         XMAX = XMAX*REC
00565                      END IF
00566                   END IF
00567 *
00568                   X( J1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X, 1 )
00569                   X( N+J1 ) = X( N+J1 ) - DDOT( J1-1, T( 1, J1 ), 1,
00570      $                        X( N+1 ), 1 )
00571                   IF( J1.GT.1 ) THEN
00572                      X( J1 ) = X( J1 ) - B( J1 )*X( N+1 )
00573                      X( N+J1 ) = X( N+J1 ) + B( J1 )*X( 1 )
00574                   END IF
00575                   XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
00576 *
00577                   Z = W
00578                   IF( J1.EQ.1 )
00579      $               Z = B( 1 )
00580 *
00581 *                 Scale if necessary to avoid overflow in
00582 *                 complex division
00583 *
00584                   TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
00585                   TMP = T( J1, J1 )
00586                   IF( TJJ.LT.SMINW ) THEN
00587                      TMP = SMINW
00588                      TJJ = SMINW
00589                      INFO = 1
00590                   END IF
00591 *
00592                   IF( TJJ.LT.ONE ) THEN
00593                      IF( XJ.GT.BIGNUM*TJJ ) THEN
00594                         REC = ONE / XJ
00595                         CALL DSCAL( N2, REC, X, 1 )
00596                         SCALE = SCALE*REC
00597                         XMAX = XMAX*REC
00598                      END IF
00599                   END IF
00600                   CALL DLADIV( X( J1 ), X( N+J1 ), TMP, -Z, SR, SI )
00601                   X( J1 ) = SR
00602                   X( J1+N ) = SI
00603                   XMAX = MAX( ABS( X( J1 ) )+ABS( X( J1+N ) ), XMAX )
00604 *
00605                ELSE
00606 *
00607 *                 2 by 2 diagonal block
00608 *
00609 *                 Scale if necessary to avoid overflow in forming the
00610 *                 right-hand side element by inner product.
00611 *
00612                   XJ = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
00613      $                 ABS( X( J2 ) )+ABS( X( N+J2 ) ) )
00614                   IF( XMAX.GT.ONE ) THEN
00615                      REC = ONE / XMAX
00616                      IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
00617      $                   ( BIGNUM-XJ ) / XMAX ) THEN
00618                         CALL DSCAL( N2, REC, X, 1 )
00619                         SCALE = SCALE*REC
00620                         XMAX = XMAX*REC
00621                      END IF
00622                   END IF
00623 *
00624                   D( 1, 1 ) = X( J1 ) - DDOT( J1-1, T( 1, J1 ), 1, X,
00625      $                        1 )
00626                   D( 2, 1 ) = X( J2 ) - DDOT( J1-1, T( 1, J2 ), 1, X,
00627      $                        1 )
00628                   D( 1, 2 ) = X( N+J1 ) - DDOT( J1-1, T( 1, J1 ), 1,
00629      $                        X( N+1 ), 1 )
00630                   D( 2, 2 ) = X( N+J2 ) - DDOT( J1-1, T( 1, J2 ), 1,
00631      $                        X( N+1 ), 1 )
00632                   D( 1, 1 ) = D( 1, 1 ) - B( J1 )*X( N+1 )
00633                   D( 2, 1 ) = D( 2, 1 ) - B( J2 )*X( N+1 )
00634                   D( 1, 2 ) = D( 1, 2 ) + B( J1 )*X( 1 )
00635                   D( 2, 2 ) = D( 2, 2 ) + B( J2 )*X( 1 )
00636 *
00637                   CALL DLALN2( .TRUE., 2, 2, SMINW, ONE, T( J1, J1 ),
00638      $                         LDT, ONE, ONE, D, 2, ZERO, W, V, 2,
00639      $                         SCALOC, XNORM, IERR )
00640                   IF( IERR.NE.0 )
00641      $               INFO = 2
00642 *
00643                   IF( SCALOC.NE.ONE ) THEN
00644                      CALL DSCAL( N2, SCALOC, X, 1 )
00645                      SCALE = SCALOC*SCALE
00646                   END IF
00647                   X( J1 ) = V( 1, 1 )
00648                   X( J2 ) = V( 2, 1 )
00649                   X( N+J1 ) = V( 1, 2 )
00650                   X( N+J2 ) = V( 2, 2 )
00651                   XMAX = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
00652      $                   ABS( X( J2 ) )+ABS( X( N+J2 ) ), XMAX )
00653 *
00654                END IF
00655 *
00656    80       CONTINUE
00657 *
00658          END IF
00659 *
00660       END IF
00661 *
00662       RETURN
00663 *
00664 *     End of DLAQTR
00665 *
00666       END