dgegv.f.html

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

</span><span class="comment">*</span><span class="comment">
</span>         IF( ILVL ) THEN
            CALL <a name="DGGBAK.503"></a><a href="dggbak.f.html#DGGBAK.1">DGGBAK</a>( <span class="string">'P'</span>, <span class="string">'L'</span>, N, ILO, IHI, WORK( ILEFT ),
     $                   WORK( IRIGHT ), N, VL, LDVL, IINFO )
            IF( IINFO.NE.0 ) THEN
               INFO = N + 8
               GO TO 120
            END IF
            DO 50 JC = 1, N
               IF( ALPHAI( JC ).LT.ZERO )
     $            GO TO 50
               TEMP = ZERO
               IF( ALPHAI( JC ).EQ.ZERO ) THEN
                  DO 10 JR = 1, N
                     TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
   10             CONTINUE
               ELSE
                  DO 20 JR = 1, N
                     TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
     $                      ABS( VL( JR, JC+1 ) ) )
   20             CONTINUE
               END IF
               IF( TEMP.LT.SAFMIN )
     $            GO TO 50
               TEMP = ONE / TEMP
               IF( ALPHAI( JC ).EQ.ZERO ) THEN
                  DO 30 JR = 1, N
                     VL( JR, JC ) = VL( JR, JC )*TEMP
   30             CONTINUE
               ELSE
                  DO 40 JR = 1, N
                     VL( JR, JC ) = VL( JR, JC )*TEMP
                     VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
   40             CONTINUE
               END IF
   50       CONTINUE
         END IF
         IF( ILVR ) THEN
            CALL <a name="DGGBAK.539"></a><a href="dggbak.f.html#DGGBAK.1">DGGBAK</a>( <span class="string">'P'</span>, <span class="string">'R'</span>, N, ILO, IHI, WORK( ILEFT ),
     $                   WORK( IRIGHT ), N, VR, LDVR, IINFO )
            IF( IINFO.NE.0 ) THEN
               INFO = N + 9
               GO TO 120
            END IF
            DO 100 JC = 1, N
               IF( ALPHAI( JC ).LT.ZERO )
     $            GO TO 100
               TEMP = ZERO
               IF( ALPHAI( JC ).EQ.ZERO ) THEN
                  DO 60 JR = 1, N
                     TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
   60             CONTINUE
               ELSE
                  DO 70 JR = 1, N
                     TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
     $                      ABS( VR( JR, JC+1 ) ) )
   70             CONTINUE
               END IF
               IF( TEMP.LT.SAFMIN )
     $            GO TO 100
               TEMP = ONE / TEMP
               IF( ALPHAI( JC ).EQ.ZERO ) THEN
                  DO 80 JR = 1, N
                     VR( JR, JC ) = VR( JR, JC )*TEMP
   80             CONTINUE
               ELSE
                  DO 90 JR = 1, N
                     VR( JR, JC ) = VR( JR, JC )*TEMP
                     VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
   90             CONTINUE
               END IF
  100       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        End of eigenvector calculation
</span><span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Undo scaling in alpha, beta
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Note: this does not give the alpha and beta for the unscaled
</span><span class="comment">*</span><span class="comment">     problem.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Un-scaling is limited to avoid underflow in alpha and beta
</span><span class="comment">*</span><span class="comment">     if they are significant.
</span><span class="comment">*</span><span class="comment">
</span>      DO 110 JC = 1, N
         ABSAR = ABS( ALPHAR( JC ) )
         ABSAI = ABS( ALPHAI( JC ) )
         ABSB = ABS( BETA( JC ) )
         SALFAR = ANRM*ALPHAR( JC )
         SALFAI = ANRM*ALPHAI( JC )
         SBETA = BNRM*BETA( JC )
         ILIMIT = .FALSE.
         SCALE = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Check for significant underflow in ALPHAI
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE.
     $       MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN
            ILIMIT = .TRUE.
            SCALE = ( ONEPLS*SAFMIN / ANRM1 ) /
     $              MAX( ONEPLS*SAFMIN, ANRM2*ABSAI )
<span class="comment">*</span><span class="comment">
</span>         ELSE IF( SALFAI.EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If insignificant underflow in ALPHAI, then make the
</span><span class="comment">*</span><span class="comment">           conjugate eigenvalue real.
</span><span class="comment">*</span><span class="comment">
</span>            IF( ALPHAI( JC ).LT.ZERO .AND. JC.GT.1 ) THEN
               ALPHAI( JC-1 ) = ZERO
            ELSE IF( ALPHAI( JC ).GT.ZERO .AND. JC.LT.N ) THEN
               ALPHAI( JC+1 ) = ZERO
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Check for significant underflow in ALPHAR
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE.
     $       MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN
            ILIMIT = .TRUE.
            SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / ANRM1 ) /
     $              MAX( ONEPLS*SAFMIN, ANRM2*ABSAR ) )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Check for significant underflow in BETA
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( SBETA ).LT.SAFMIN .AND. ABSB.GE.
     $       MAX( SAFMIN, EPS*ABSAR, EPS*ABSAI ) ) THEN
            ILIMIT = .TRUE.
            SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / BNRM1 ) /
     $              MAX( ONEPLS*SAFMIN, BNRM2*ABSB ) )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Check for possible overflow when limiting scaling
</span><span class="comment">*</span><span class="comment">
</span>         IF( ILIMIT ) THEN
            TEMP = ( SCALE*SAFMIN )*MAX( ABS( SALFAR ), ABS( SALFAI ),
     $             ABS( SBETA ) )
            IF( TEMP.GT.ONE )
     $         SCALE = SCALE / TEMP
            IF( SCALE.LT.ONE )
     $         ILIMIT = .FALSE.
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Recompute un-scaled ALPHAR, ALPHAI, BETA if necessary.
</span><span class="comment">*</span><span class="comment">
</span>         IF( ILIMIT ) THEN
            SALFAR = ( SCALE*ALPHAR( JC ) )*ANRM
            SALFAI = ( SCALE*ALPHAI( JC ) )*ANRM
            SBETA = ( SCALE*BETA( JC ) )*BNRM
         END IF
         ALPHAR( JC ) = SALFAR
         ALPHAI( JC ) = SALFAI
         BETA( JC ) = SBETA
  110 CONTINUE
<span class="comment">*</span><span class="comment">
</span>  120 CONTINUE
      WORK( 1 ) = LWKOPT
<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="DGEGV.663"></a><a href="dgegv.f.html#DGEGV.1">DGEGV</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>