sgeevx.f.html

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

<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Generate orthogonal matrix in VR
</span><span class="comment">*</span><span class="comment">        (Workspace: need 2*N-1, prefer N+(N-1)*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SORGHR.421"></a><a href="sorghr.f.html#SORGHR.1">SORGHR</a>( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
     $                LWORK-IWRK+1, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Perform QR iteration, accumulating Schur vectors in VR
</span><span class="comment">*</span><span class="comment">        (Workspace: need 1, prefer HSWORK (see comments) )
</span><span class="comment">*</span><span class="comment">
</span>         IWRK = ITAU
         CALL <a name="SHSEQR.428"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>( <span class="string">'S'</span>, <span class="string">'V'</span>, N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
     $                WORK( IWRK ), LWORK-IWRK+1, INFO )
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute eigenvalues only
</span><span class="comment">*</span><span class="comment">        If condition numbers desired, compute Schur form
</span><span class="comment">*</span><span class="comment">
</span>         IF( WNTSNN ) THEN
            JOB = <span class="string">'E'</span>
         ELSE
            JOB = <span class="string">'S'</span>
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        (Workspace: need 1, prefer HSWORK (see comments) )
</span><span class="comment">*</span><span class="comment">
</span>         IWRK = ITAU
         CALL <a name="SHSEQR.445"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>( JOB, <span class="string">'N'</span>, N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
     $                WORK( IWRK ), LWORK-IWRK+1, INFO )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If INFO &gt; 0 from <a name="SHSEQR.449"></a><a href="shseqr.f.html#SHSEQR.1">SHSEQR</a>, then quit
</span><span class="comment">*</span><span class="comment">
</span>      IF( INFO.GT.0 )
     $   GO TO 50
<span class="comment">*</span><span class="comment">
</span>      IF( WANTVL .OR. WANTVR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute left and/or right eigenvectors
</span><span class="comment">*</span><span class="comment">        (Workspace: need 3*N)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="STREVC.459"></a><a href="strevc.f.html#STREVC.1">STREVC</a>( SIDE, <span class="string">'B'</span>, SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
     $                N, NOUT, WORK( IWRK ), IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute condition numbers if desired
</span><span class="comment">*</span><span class="comment">     (Workspace: need N*N+6*N unless SENSE = 'E')
</span><span class="comment">*</span><span class="comment">
</span>      IF( .NOT.WNTSNN ) THEN
         CALL <a name="STRSNA.467"></a><a href="strsna.f.html#STRSNA.1">STRSNA</a>( SENSE, <span class="string">'A'</span>, SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
     $                RCONDE, RCONDV, N, NOUT, WORK( IWRK ), N, IWORK,
     $                ICOND )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( WANTVL ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Undo balancing of left eigenvectors
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SGEBAK.476"></a><a href="sgebak.f.html#SGEBAK.1">SGEBAK</a>( BALANC, <span class="string">'L'</span>, N, ILO, IHI, SCALE, N, VL, LDVL,
     $                IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Normalize left eigenvectors and make largest component real
</span><span class="comment">*</span><span class="comment">
</span>         DO 20 I = 1, N
            IF( WI( I ).EQ.ZERO ) THEN
               SCL = ONE / SNRM2( N, VL( 1, I ), 1 )
               CALL SSCAL( N, SCL, VL( 1, I ), 1 )
            ELSE IF( WI( I ).GT.ZERO ) THEN
               SCL = ONE / <a name="SLAPY2.486"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( SNRM2( N, VL( 1, I ), 1 ),
     $               SNRM2( N, VL( 1, I+1 ), 1 ) )
               CALL SSCAL( N, SCL, VL( 1, I ), 1 )
               CALL SSCAL( N, SCL, VL( 1, I+1 ), 1 )
               DO 10 K = 1, N
                  WORK( K ) = VL( K, I )**2 + VL( K, I+1 )**2
   10          CONTINUE
               K = ISAMAX( N, WORK, 1 )
               CALL <a name="SLARTG.494"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( VL( K, I ), VL( K, I+1 ), CS, SN, R )
               CALL SROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
               VL( K, I+1 ) = ZERO
            END IF
   20    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( WANTVR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Undo balancing of right eigenvectors
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SGEBAK.505"></a><a href="sgebak.f.html#SGEBAK.1">SGEBAK</a>( BALANC, <span class="string">'R'</span>, N, ILO, IHI, SCALE, N, VR, LDVR,
     $                IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Normalize right eigenvectors and make largest component real
</span><span class="comment">*</span><span class="comment">
</span>         DO 40 I = 1, N
            IF( WI( I ).EQ.ZERO ) THEN
               SCL = ONE / SNRM2( N, VR( 1, I ), 1 )
               CALL SSCAL( N, SCL, VR( 1, I ), 1 )
            ELSE IF( WI( I ).GT.ZERO ) THEN
               SCL = ONE / <a name="SLAPY2.515"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( SNRM2( N, VR( 1, I ), 1 ),
     $               SNRM2( N, VR( 1, I+1 ), 1 ) )
               CALL SSCAL( N, SCL, VR( 1, I ), 1 )
               CALL SSCAL( N, SCL, VR( 1, I+1 ), 1 )
               DO 30 K = 1, N
                  WORK( K ) = VR( K, I )**2 + VR( K, I+1 )**2
   30          CONTINUE
               K = ISAMAX( N, WORK, 1 )
               CALL <a name="SLARTG.523"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( VR( K, I ), VR( K, I+1 ), CS, SN, R )
               CALL SROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
               VR( K, I+1 ) = ZERO
            END IF
   40    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Undo scaling if necessary
</span><span class="comment">*</span><span class="comment">
</span>   50 CONTINUE
      IF( SCALEA ) THEN
         CALL <a name="SLASCL.534"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
     $                MAX( N-INFO, 1 ), IERR )
         CALL <a name="SLASCL.536"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
     $                MAX( N-INFO, 1 ), IERR )
         IF( INFO.EQ.0 ) THEN
            IF( ( WNTSNV .OR. WNTSNB ) .AND. ICOND.EQ.0 )
     $         CALL <a name="SLASCL.540"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, N, 1, RCONDV, N,
     $                      IERR )
         ELSE
            CALL <a name="SLASCL.543"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
     $                   IERR )
            CALL <a name="SLASCL.545"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
     $                   IERR )
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = MAXWRK
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SGEEVX.553"></a><a href="sgeevx.f.html#SGEEVX.1">SGEEVX</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>