sgges.f.html

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

HTML
583 行
字号
123
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize VSL
</span><span class="comment">*</span><span class="comment">     (Workspace: need N, prefer N*NB)
</span><span class="comment">*</span><span class="comment">
</span>      IF( ILVSL ) THEN
         CALL <a name="SLASET.377"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, VSL, LDVSL )
         IF( IROWS.GT.1 ) THEN
            CALL <a name="SLACPY.379"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'L'</span>, IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
     $                   VSL( ILO+1, ILO ), LDVSL )
         END IF
         CALL <a name="SORGQR.382"></a><a href="sorgqr.f.html#SORGQR.1">SORGQR</a>( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
     $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize VSR
</span><span class="comment">*</span><span class="comment">
</span>      IF( ILVSR )
     $   CALL <a name="SLASET.389"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, VSR, LDVSR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Reduce to generalized Hessenberg form
</span><span class="comment">*</span><span class="comment">     (Workspace: none needed)
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="SGGHRD.394"></a><a href="sgghrd.f.html#SGGHRD.1">SGGHRD</a>( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
     $             LDVSL, VSR, LDVSR, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Perform QZ algorithm, computing Schur vectors if desired
</span><span class="comment">*</span><span class="comment">     (Workspace: need N)
</span><span class="comment">*</span><span class="comment">
</span>      IWRK = ITAU
      CALL <a name="SHGEQZ.401"></a><a href="shgeqz.f.html#SHGEQZ.1">SHGEQZ</a>( <span class="string">'S'</span>, JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
     $             ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
     $             WORK( IWRK ), LWORK+1-IWRK, IERR )
      IF( IERR.NE.0 ) THEN
         IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
            INFO = IERR
         ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
            INFO = IERR - N
         ELSE
            INFO = N + 1
         END IF
         GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Sort eigenvalues ALPHA/BETA if desired
</span><span class="comment">*</span><span class="comment">     (Workspace: need 4*N+16 )
</span><span class="comment">*</span><span class="comment">
</span>      SDIM = 0
      IF( WANTST ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Undo scaling on eigenvalues before SELCTGing
</span><span class="comment">*</span><span class="comment">
</span>         IF( ILASCL ) THEN
            CALL <a name="SLASCL.424"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
     $                   IERR )
            CALL <a name="SLASCL.426"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
     $                   IERR )
         END IF
         IF( ILBSCL )
     $      CALL <a name="SLASCL.430"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Select eigenvalues
</span><span class="comment">*</span><span class="comment">
</span>         DO 10 I = 1, N
            BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
   10    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         CALL <a name="STGSEN.438"></a><a href="stgsen.f.html#STGSEN.1">STGSEN</a>( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHAR,
     $                ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL,
     $                PVSR, DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
     $                IERR )
         IF( IERR.EQ.1 )
     $      INFO = N + 3
<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">     Apply back-permutation to VSL and VSR
</span><span class="comment">*</span><span class="comment">     (Workspace: none needed)
</span><span class="comment">*</span><span class="comment">
</span>      IF( ILVSL )
     $   CALL <a name="SGGBAK.451"></a><a href="sggbak.f.html#SGGBAK.1">SGGBAK</a>( <span class="string">'P'</span>, <span class="string">'L'</span>, N, ILO, IHI, WORK( ILEFT ),
     $                WORK( IRIGHT ), N, VSL, LDVSL, IERR )
<span class="comment">*</span><span class="comment">
</span>      IF( ILVSR )
     $   CALL <a name="SGGBAK.455"></a><a href="sggbak.f.html#SGGBAK.1">SGGBAK</a>( <span class="string">'P'</span>, <span class="string">'R'</span>, N, ILO, IHI, WORK( ILEFT ),
     $                WORK( IRIGHT ), N, VSR, LDVSR, IERR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check if unscaling would cause over/underflow, if so, rescale 
</span><span class="comment">*</span><span class="comment">     (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of 
</span><span class="comment">*</span><span class="comment">     B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
</span><span class="comment">*</span><span class="comment">
</span>      IF( ILASCL )THEN
         DO 50 I = 1, N 
            IF( ALPHAI( I ).NE.ZERO ) THEN 
               IF( ( ALPHAR( I )/SAFMAX ).GT.( ANRMTO/ANRM ) .OR.
     $             ( SAFMIN/ALPHAR( I ) ).GT.( ANRM/ANRMTO ) ) THEN
                  WORK( 1 ) = ABS( A( I, I )/ALPHAR( I ) )
                  BETA( I ) = BETA( I )*WORK( 1 )
                  ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
                  ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
               ELSE IF( ( ALPHAI( I )/SAFMAX ).GT.( ANRMTO/ANRM ) .OR.
     $             ( SAFMIN/ALPHAI( I ) ).GT.( ANRM/ANRMTO ) ) THEN
                  WORK( 1 ) = ABS( A( I, I+1 )/ALPHAI( I ) )
                  BETA( I ) = BETA( I )*WORK( 1 )
                  ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
                  ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
               END IF
            END IF
   50    CONTINUE
      END IF 
<span class="comment">*</span><span class="comment">
</span>      IF( ILBSCL )THEN 
         DO 60 I = 1, N
            IF( ALPHAI( I ).NE.ZERO ) THEN
                IF( ( BETA( I )/SAFMAX ).GT.( BNRMTO/BNRM ) .OR.
     $              ( SAFMIN/BETA( I ) ).GT.( BNRM/BNRMTO ) ) THEN
                   WORK( 1 ) = ABS(B( I, I )/BETA( I ))
                   BETA( I ) = BETA( I )*WORK( 1 )
                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
                END IF 
             END IF
   60    CONTINUE 
      END IF 
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Undo scaling
</span><span class="comment">*</span><span class="comment">
</span>      IF( ILASCL ) THEN
         CALL <a name="SLASCL.499"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'H'</span>, 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
         CALL <a name="SLASCL.500"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
         CALL <a name="SLASCL.501"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ILBSCL ) THEN
         CALL <a name="SLASCL.505"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'U'</span>, 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
         CALL <a name="SLASCL.506"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( WANTST ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Check if reordering is correct
</span><span class="comment">*</span><span class="comment">
</span>         LASTSL = .TRUE.
         LST2SL = .TRUE.
         SDIM = 0
         IP = 0
         DO 30 I = 1, N
            CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
            IF( ALPHAI( I ).EQ.ZERO ) THEN
               IF( CURSL )
     $            SDIM = SDIM + 1
               IP = 0
               IF( CURSL .AND. .NOT.LASTSL )
     $            INFO = N + 2
            ELSE
               IF( IP.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Last eigenvalue of conjugate pair
</span><span class="comment">*</span><span class="comment">
</span>                  CURSL = CURSL .OR. LASTSL
                  LASTSL = CURSL
                  IF( CURSL )
     $               SDIM = SDIM + 2
                  IP = -1
                  IF( CURSL .AND. .NOT.LST2SL )
     $               INFO = N + 2
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 First eigenvalue of conjugate pair
</span><span class="comment">*</span><span class="comment">
</span>                  IP = 1
               END IF
            END IF
            LST2SL = LASTSL
            LASTSL = CURSL
   30    CONTINUE
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span>   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = MAXWRK
<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="SGGES.556"></a><a href="sgges.f.html#SGGES.1">SGGES</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
123

sgges.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 sgges.f.html 源码文件，采用 HTML 编程语言编写，共 583 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。

虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?