slaqtr.f.html

famous linear algebra library (LAPACK) ports to windows

<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale x if necessary to avoid overflow when adding a
</span><span class="comment">*</span><span class="comment">                 multiple of column j1 of T.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( XJ.GT.ONE ) THEN
                     REC = ONE / XJ
                     IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
                        CALL SSCAL( N2, REC, X, 1 )
                        SCALE = SCALE*REC
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  IF( J1.GT.1 ) THEN
                     CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
                     CALL SAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
     $                           X( N+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span>                     X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 )
                     X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 )
<span class="comment">*</span><span class="comment">
</span>                     XMAX = ZERO
                     DO 50 K = 1, J1 - 1
                        XMAX = MAX( XMAX, ABS( X( K ) )+
     $                         ABS( X( K+N ) ) )
   50                CONTINUE
                  END IF
<span class="comment">*</span><span class="comment">
</span>               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Meet 2 by 2 diagonal block
</span><span class="comment">*</span><span class="comment">
</span>                  D( 1, 1 ) = X( J1 )
                  D( 2, 1 ) = X( J2 )
                  D( 1, 2 ) = X( N+J1 )
                  D( 2, 2 ) = X( N+J2 )
                  CALL <a name="SLALN2.477"></a><a href="slaln2.f.html#SLALN2.1">SLALN2</a>( .FALSE., 2, 2, SMINW, ONE, T( J1, J1 ),
     $                         LDT, ONE, ONE, D, 2, ZERO, -W, V, 2,
     $                         SCALOC, XNORM, IERR )
                  IF( IERR.NE.0 )
     $               INFO = 2
<span class="comment">*</span><span class="comment">
</span>                  IF( SCALOC.NE.ONE ) THEN
                     CALL SSCAL( 2*N, SCALOC, X, 1 )
                     SCALE = SCALOC*SCALE
                  END IF
                  X( J1 ) = V( 1, 1 )
                  X( J2 ) = V( 2, 1 )
                  X( N+J1 ) = V( 1, 2 )
                  X( N+J2 ) = V( 2, 2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale X(J1), .... to avoid overflow in
</span><span class="comment">*</span><span class="comment">                 updating right hand side.
</span><span class="comment">*</span><span class="comment">
</span>                  XJ = MAX( ABS( V( 1, 1 ) )+ABS( V( 1, 2 ) ),
     $                 ABS( V( 2, 1 ) )+ABS( V( 2, 2 ) ) )
                  IF( XJ.GT.ONE ) THEN
                     REC = ONE / XJ
                     IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
     $                   ( BIGNUM-XMAX )*REC ) THEN
                        CALL SSCAL( N2, REC, X, 1 )
                        SCALE = SCALE*REC
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Update the right-hand side.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( J1.GT.1 ) THEN
                     CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
                     CALL SAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
<span class="comment">*</span><span class="comment">
</span>                     CALL SAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
     $                           X( N+1 ), 1 )
                     CALL SAXPY( J1-1, -X( N+J2 ), T( 1, J2 ), 1,
     $                           X( N+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span>                     X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 ) +
     $                        B( J2 )*X( N+J2 )
                     X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 ) -
     $                          B( J2 )*X( J2 )
<span class="comment">*</span><span class="comment">
</span>                     XMAX = ZERO
                     DO 60 K = 1, J1 - 1
                        XMAX = MAX( ABS( X( K ) )+ABS( X( K+N ) ),
     $                         XMAX )
   60                CONTINUE
                  END IF
<span class="comment">*</span><span class="comment">
</span>               END IF
   70       CONTINUE
<span class="comment">*</span><span class="comment">
</span>         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Solve (T + iB)'*(p+iq) = c+id
</span><span class="comment">*</span><span class="comment">
</span>            JNEXT = 1
            DO 80 J = 1, N
               IF( J.LT.JNEXT )
     $            GO TO 80
               J1 = J
               J2 = J
               JNEXT = J + 1
               IF( J.LT.N ) THEN
                  IF( T( J+1, J ).NE.ZERO ) THEN
                     J2 = J + 1
                     JNEXT = J + 2
                  END IF
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( J1.EQ.J2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 1 by 1 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale if necessary to avoid overflow in forming the
</span><span class="comment">*</span><span class="comment">                 right-hand side element by inner product.
</span><span class="comment">*</span><span class="comment">
</span>                  XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
                  IF( XMAX.GT.ONE ) THEN
                     REC = ONE / XMAX
                     IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
                        CALL SSCAL( N2, REC, X, 1 )
                        SCALE = SCALE*REC
                        XMAX = XMAX*REC
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  X( J1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X, 1 )
                  X( N+J1 ) = X( N+J1 ) - SDOT( J1-1, T( 1, J1 ), 1,
     $                        X( N+1 ), 1 )
                  IF( J1.GT.1 ) THEN
                     X( J1 ) = X( J1 ) - B( J1 )*X( N+1 )
                     X( N+J1 ) = X( N+J1 ) + B( J1 )*X( 1 )
                  END IF
                  XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
<span class="comment">*</span><span class="comment">
</span>                  Z = W
                  IF( J1.EQ.1 )
     $               Z = B( 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale if necessary to avoid overflow in
</span><span class="comment">*</span><span class="comment">                 complex division
</span><span class="comment">*</span><span class="comment">
</span>                  TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
                  TMP = T( J1, J1 )
                  IF( TJJ.LT.SMINW ) THEN
                     TMP = SMINW
                     TJJ = SMINW
                     INFO = 1
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  IF( TJJ.LT.ONE ) THEN
                     IF( XJ.GT.BIGNUM*TJJ ) THEN
                        REC = ONE / XJ
                        CALL SSCAL( N2, REC, X, 1 )
                        SCALE = SCALE*REC
                        XMAX = XMAX*REC
                     END IF
                  END IF
                  CALL <a name="SLADIV.599"></a><a href="sladiv.f.html#SLADIV.1">SLADIV</a>( X( J1 ), X( N+J1 ), TMP, -Z, SR, SI )
                  X( J1 ) = SR
                  X( J1+N ) = SI
                  XMAX = MAX( ABS( X( J1 ) )+ABS( X( J1+N ) ), XMAX )
<span class="comment">*</span><span class="comment">
</span>               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 2 by 2 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale if necessary to avoid overflow in forming the
</span><span class="comment">*</span><span class="comment">                 right-hand side element by inner product.
</span><span class="comment">*</span><span class="comment">
</span>                  XJ = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
     $                 ABS( X( J2 ) )+ABS( X( N+J2 ) ) )
                  IF( XMAX.GT.ONE ) THEN
                     REC = ONE / XMAX
                     IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
     $                   ( BIGNUM-XJ ) / XMAX ) THEN
                        CALL SSCAL( N2, REC, X, 1 )
                        SCALE = SCALE*REC
                        XMAX = XMAX*REC
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  D( 1, 1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X,
     $                        1 )
                  D( 2, 1 ) = X( J2 ) - SDOT( J1-1, T( 1, J2 ), 1, X,
     $                        1 )
                  D( 1, 2 ) = X( N+J1 ) - SDOT( J1-1, T( 1, J1 ), 1,
     $                        X( N+1 ), 1 )
                  D( 2, 2 ) = X( N+J2 ) - SDOT( J1-1, T( 1, J2 ), 1,
     $                        X( N+1 ), 1 )
                  D( 1, 1 ) = D( 1, 1 ) - B( J1 )*X( N+1 )
                  D( 2, 1 ) = D( 2, 1 ) - B( J2 )*X( N+1 )
                  D( 1, 2 ) = D( 1, 2 ) + B( J1 )*X( 1 )
                  D( 2, 2 ) = D( 2, 2 ) + B( J2 )*X( 1 )
<span class="comment">*</span><span class="comment">
</span>                  CALL <a name="SLALN2.636"></a><a href="slaln2.f.html#SLALN2.1">SLALN2</a>( .TRUE., 2, 2, SMINW, ONE, T( J1, J1 ),
     $                         LDT, ONE, ONE, D, 2, ZERO, W, V, 2,
     $                         SCALOC, XNORM, IERR )
                  IF( IERR.NE.0 )
     $               INFO = 2
<span class="comment">*</span><span class="comment">
</span>                  IF( SCALOC.NE.ONE ) THEN
                     CALL SSCAL( N2, SCALOC, X, 1 )
                     SCALE = SCALOC*SCALE
                  END IF
                  X( J1 ) = V( 1, 1 )
                  X( J2 ) = V( 2, 1 )
                  X( N+J1 ) = V( 1, 2 )
                  X( N+J2 ) = V( 2, 2 )
                  XMAX = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
     $                   ABS( X( J2 ) )+ABS( X( N+J2 ) ), XMAX )
<span class="comment">*</span><span class="comment">
</span>               END IF
<span class="comment">*</span><span class="comment">
</span>   80       CONTINUE
<span class="comment">*</span><span class="comment">
</span>         END IF
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SLAQTR.663"></a><a href="slaqtr.f.html#SLAQTR.1">SLAQTR</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>