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