dggbal.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 494 行 · 第 1/2 页

      DO 150 J = K, L
         DO 110 I = K, LM1
            IP1 = I + 1
            IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
     $         GO TO 120
  110    CONTINUE
         I = L
         GO TO 140
  120    CONTINUE
         DO 130 I = IP1, L
            IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
     $         GO TO 150
  130    CONTINUE
         I = IP1 - 1
  140    CONTINUE
         M = K
         IFLOW = 2
         GO TO 160
  150 CONTINUE
      GO TO 190
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Permute rows M and I
</span><span class="comment">*</span><span class="comment">
</span>  160 CONTINUE
      LSCALE( M ) = I
      IF( I.EQ.M )
     $   GO TO 170
      CALL DSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
      CALL DSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Permute columns M and J
</span><span class="comment">*</span><span class="comment">
</span>  170 CONTINUE
      RSCALE( M ) = J
      IF( J.EQ.M )
     $   GO TO 180
      CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
      CALL DSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
<span class="comment">*</span><span class="comment">
</span>  180 CONTINUE
      GO TO ( 20, 90 )IFLOW
<span class="comment">*</span><span class="comment">
</span>  190 CONTINUE
      ILO = K
      IHI = L
<span class="comment">*</span><span class="comment">
</span>      IF( <a name="LSAME.275"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOB, <span class="string">'P'</span> ) ) THEN
         DO 195 I = ILO, IHI
            LSCALE( I ) = ONE
            RSCALE( I ) = ONE
  195    CONTINUE
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ILO.EQ.IHI )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Balance the submatrix in rows ILO to IHI.
</span><span class="comment">*</span><span class="comment">
</span>      NR = IHI - ILO + 1
      DO 200 I = ILO, IHI
         RSCALE( I ) = ZERO
         LSCALE( I ) = ZERO
<span class="comment">*</span><span class="comment">
</span>         WORK( I ) = ZERO
         WORK( I+N ) = ZERO
         WORK( I+2*N ) = ZERO
         WORK( I+3*N ) = ZERO
         WORK( I+4*N ) = ZERO
         WORK( I+5*N ) = ZERO
  200 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute right side vector in resulting linear equations
</span><span class="comment">*</span><span class="comment">
</span>      BASL = LOG10( SCLFAC )
      DO 240 I = ILO, IHI
         DO 230 J = ILO, IHI
            TB = B( I, J )
            TA = A( I, J )
            IF( TA.EQ.ZERO )
     $         GO TO 210
            TA = LOG10( ABS( TA ) ) / BASL
  210       CONTINUE
            IF( TB.EQ.ZERO )
     $         GO TO 220
            TB = LOG10( ABS( TB ) ) / BASL
  220       CONTINUE
            WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
            WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
  230    CONTINUE
  240 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      COEF = ONE / DBLE( 2*NR )
      COEF2 = COEF*COEF
      COEF5 = HALF*COEF2
      NRP2 = NR + 2
      BETA = ZERO
      IT = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Start generalized conjugate gradient iteration
</span><span class="comment">*</span><span class="comment">
</span>  250 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
     $        DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
<span class="comment">*</span><span class="comment">
</span>      EW = ZERO
      EWC = ZERO
      DO 260 I = ILO, IHI
         EW = EW + WORK( I+4*N )
         EWC = EWC + WORK( I+5*N )
  260 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
      IF( GAMMA.EQ.ZERO )
     $   GO TO 350
      IF( IT.NE.1 )
     $   BETA = GAMMA / PGAMMA
      T = COEF5*( EWC-THREE*EW )
      TC = COEF5*( EW-THREE*EWC )
<span class="comment">*</span><span class="comment">
</span>      CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
      CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
<span class="comment">*</span><span class="comment">
</span>      CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
      CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
<span class="comment">*</span><span class="comment">
</span>      DO 270 I = ILO, IHI
         WORK( I ) = WORK( I ) + TC
         WORK( I+N ) = WORK( I+N ) + T
  270 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Apply matrix to vector
</span><span class="comment">*</span><span class="comment">
</span>      DO 300 I = ILO, IHI
         KOUNT = 0
         SUM = ZERO
         DO 290 J = ILO, IHI
            IF( A( I, J ).EQ.ZERO )
     $         GO TO 280
            KOUNT = KOUNT + 1
            SUM = SUM + WORK( J )
  280       CONTINUE
            IF( B( I, J ).EQ.ZERO )
     $         GO TO 290
            KOUNT = KOUNT + 1
            SUM = SUM + WORK( J )
  290    CONTINUE
         WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
  300 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      DO 330 J = ILO, IHI
         KOUNT = 0
         SUM = ZERO
         DO 320 I = ILO, IHI
            IF( A( I, J ).EQ.ZERO )
     $         GO TO 310
            KOUNT = KOUNT + 1
            SUM = SUM + WORK( I+N )
  310       CONTINUE
            IF( B( I, J ).EQ.ZERO )
     $         GO TO 320
            KOUNT = KOUNT + 1
            SUM = SUM + WORK( I+N )
  320    CONTINUE
         WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
  330 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
     $      DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
      ALPHA = GAMMA / SUM
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Determine correction to current iteration
</span><span class="comment">*</span><span class="comment">
</span>      CMAX = ZERO
      DO 340 I = ILO, IHI
         COR = ALPHA*WORK( I+N )
         IF( ABS( COR ).GT.CMAX )
     $      CMAX = ABS( COR )
         LSCALE( I ) = LSCALE( I ) + COR
         COR = ALPHA*WORK( I )
         IF( ABS( COR ).GT.CMAX )
     $      CMAX = ABS( COR )
         RSCALE( I ) = RSCALE( I ) + COR
  340 CONTINUE
      IF( CMAX.LT.HALF )
     $   GO TO 350
<span class="comment">*</span><span class="comment">
</span>      CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
      CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
<span class="comment">*</span><span class="comment">
</span>      PGAMMA = GAMMA
      IT = IT + 1
      IF( IT.LE.NRP2 )
     $   GO TO 250
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End generalized conjugate gradient iteration
</span><span class="comment">*</span><span class="comment">
</span>  350 CONTINUE
      SFMIN = <a name="DLAMCH.428"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> )
      SFMAX = ONE / SFMIN
      LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
      LSFMAX = INT( LOG10( SFMAX ) / BASL )
      DO 360 I = ILO, IHI
         IRAB = IDAMAX( N-ILO+1, A( I, ILO ), LDA )
         RAB = ABS( A( I, IRAB+ILO-1 ) )
         IRAB = IDAMAX( N-ILO+1, B( I, ILO ), LDB )
         RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
         LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
         IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )
         IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
         LSCALE( I ) = SCLFAC**IR
         ICAB = IDAMAX( IHI, A( 1, I ), 1 )
         CAB = ABS( A( ICAB, I ) )
         ICAB = IDAMAX( IHI, B( 1, I ), 1 )
         CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
         LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
         JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) )
         JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
         RSCALE( I ) = SCLFAC**JC
  360 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Row scaling of matrices A and B
</span><span class="comment">*</span><span class="comment">
</span>      DO 370 I = ILO, IHI
         CALL DSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
         CALL DSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
  370 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Column scaling of matrices A and B
</span><span class="comment">*</span><span class="comment">
</span>      DO 380 J = ILO, IHI
         CALL DSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
         CALL DSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
  380 CONTINUE
<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="DGGBAL.467"></a><a href="dggbal.f.html#DGGBAL.1">DGGBAL</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>