* PACKAGE LEIGEN  !" 力学 固有値解析
*
*"  [HIS] 92/12/14(takepiro) 「行列計算ソフトウェア」使用
*"        92/12/18(takepiro)  Numerical Recipies 使用
*"        93/05/10(takepiro)  エラー処置変更
*"        93/06/11(takepiro)  固有関数計算スイッチ導入
*
**********************************************************************
      SUBROUTINE LEIGEN
     M            ( AMAT   , 
     O              EIGVLR , EIGVLI  ,
     O              EIGVCR , EIGVCI  ,
     O              IERR   , 
     I              MATDIM , OEIVEC    )
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 
      INCLUDE   (ZLDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.F"                    !" 
#include        "zldim.F"                    !" 係数行列の大きさ
#endif
*
*   [MODIFY] 
      REAL      AMAT  ( MATDMX, MATDMX )     !" 係数行列
*
*   [OUTPUT] 
      REAL      EIGVLR( MATDMX )             !" 固有値実数部
      REAL      EIGVLI( MATDMX )             !" 固有値虚数部
      REAL      EIGVCR( MATDMX, MATDMX )     !" 固有ベクトル実数部
      REAL      EIGVCI( MATDMX, MATDMX )     !" 固有ベクトル虚数部
      INTEGER   IERR                         !" サブルーチンエラーコード
*
*   [INPUT] 
      INTEGER   MATDIM                       !" 計算する次元
      LOGICAL   OEIVEC                       !" 固有関数計算スイッチ
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/  Z
*
      REAL      SCALE ( MATDMX )             !" スケーリング情報
      INTEGER   LOW, IGH                     !" 
      REAL      INT   ( MATDMX )             !" ヘッセンベルグ変換情報
      REAL      Z     ( MATDMX, MATDMX )     !" 固有ベクトル
*
      CHARACTER HERR*8                       !" エラー番号出力用文字変数
      INTEGER   NERR                         !" 文字列の長さ
      INTEGER   LENC                         !" 文字の長さを返す関数
*
      SAVE      
*
*"  < 1. バランシング >
*
      CALL BALANC
     I      (  MATDMX , MATDIM , 
     M         AMAT   ,
     O         LOW    , IGH    ,SCALE  )
*
*"  < 2. ヘッセンベルグ行列変換 >
*
      CALL ELMHES
     I      (  MATDMX , MATDIM , 
     I         LOW    , IGH    , 
     M         AMAT   , INT      )
*
      CALL ELTRAN
     I      (  MATDMX , MATDIM , 
     I         LOW    , IGH    , 
     I         AMAT   , INT    , 
     O         Z                 )

*
*"  < 3. 固有値の計算 >
*
      CALL HQR2
     I      (  MATDMX , MATDIM , 
     I         LOW    , IGH    , 
     I         AMAT   , 
     O         EIGVLR , EIGVLI , 
     M         Z      , 
     O         IERR              )
*
      IF ( IERR .NE. 0 ) THEN
         WRITE( HERR , '(I5)' ) IERR
         CALL CLADJ( HERR )
         NERR = LENC( HERR )
         CALL MSGDMP( 'W', 'LEIGEN', 
     &         'IERR = '//HERR(1:NERR)//' : EIGENVALUES NOT CORRECT ' )
      ENDIF
*
*"  < 4. 固有ベクトルの計算 >
*
      IF ( OEIVEC ) THEN
         CALL BALBAK
     I      (  MATDMX , MATDIM , 
     I         LOW    , IGH    , 
     I         SCALE  , MATDIM , 
     M         Z                  )
      ENDIF
*
*"  < 5. 固有ベクトルの並び替え > 
*
      CALL EIGSRT
     M      (  EIGVLR , EIGVLI , 
     O         EIGVCR , EIGVCI ,
     I         Z      , 
     I         MATDIM , OEIVEC   )
*
*"  < 5. 終了 > 
*
      RETURN
      END
************************************************************************
      SUBROUTINE EIGSRT      !" 固有値固有ベクトル並べ変え
     M      (  EIGVLR , EIGVLI , 
     O         EIGVCR , EIGVCI , 
     I         Z      , 
     I         MATDIM , OEIVEC   )
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 
      INCLUDE   (ZLDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.F"                    !" 
#include        "zldim.F"                    !" 係数行列の大きさ
#endif
*
*   [OUTPUT] 
      REAL      EIGVLR( MATDMX )             !" 固有値実数部
      REAL      EIGVLI( MATDMX )             !" 固有値虚数部
      REAL      EIGVCR( MATDMX, MATDMX )     !" 固有ベクトル実数部
      REAL      EIGVCI( MATDMX, MATDMX )     !" 固有ベクトル虚数部
*
*   [INPUT]
      REAL      Z     ( MATDMX, MATDMX )     !" 固有ベクトル
*
      INTEGER   MATDIM                       !" 計算する次元
      LOGICAL   OEIVEC                       !" 固有関数計算スイッチ
*
*   [INTERNAL WORK]
      INTEGER   J, K
      INTEGER   INDX  ( MATDMX )             !" 順番格納配列
*
*" < 1. 固有ベクトルの入れ換え >
*
      IF ( OEIVEC )THEN
         DO 1000 K = 1, MATDIM
            IF ( EIGVLI(K) .EQ. 0.0 ) THEN
               DO 1100 J = 1, MATDIM
                  EIGVCR( J,K ) = Z( J,K )
                  EIGVCI( J,K ) = 0.
 1100          CONTINUE
*
            ELSE IF ( EIGVLI(K) .GE. 0.0 ) THEN
               DO 1200 J = 1, MATDIM
                  EIGVCR( J,K ) =  Z( J,K ) 
                  EIGVCI( J,K ) =  Z( J,K+1 )
 1200          CONTINUE
*
            ELSE
               DO 1300 J = 1, MATDIM
                  EIGVCR(J,K) =    EIGVCR( J,K-1 )
                  EIGVCI(J,K) =  - EIGVCI( J,K-1 )
 1300          CONTINUE
            ENDIF
 1000    CONTINUE
      ENDIF
*
*"  < 2. 固有値のソート >
*
      CALL INDEXX                                   !" 固有値実数部の小さい順
     I       ( MATDIM , EIGVLR ,
     O         INDX              )
*
      CALL COPY( Z, EIGVLR, MATDIM )             !" 固有値実数部
      DO 2100 J = 1, MATDIM                         !" 固有値実数部の大きい順
         EIGVLR( J ) = Z( INDX( MATDIM-J+1 ),1  )
 2100 CONTINUE
*
      CALL COPY( Z, EIGVLI, MATDIM )             !" 固有値虚数部
      DO 2200 J = 1, MATDIM                         !" 固有値実数部の大きい順
         EIGVLI( J ) = Z( INDX( MATDIM-J+1 ),1 )
 2200 CONTINUE
*
*"  < 3. 固有ベクトルのソート >
*
      IF ( OEIVEC )THEN                 
         CALL COPY( Z, EIGVCR, MATDMX*MATDIM )  !" 固有ベクトル実数部
         DO 2300 K = 1, MATDIM                     !" 固有値実数部の大きい順
            DO 2300 J = 1, MATDIM
               EIGVCR( J,K ) = Z( J,INDX( MATDIM-K+1 ) )
 2300    CONTINUE
*
         CALL COPY( Z, EIGVCI, MATDMX*MATDIM )  !" 固有ベクトル虚数部
         DO 2400 K = 1, MATDIM                     !" 固有値実数部の大きい順
            DO 2400 J = 1, MATDIM
               EIGVCI( J,K ) = Z( J,INDX( MATDIM-K+1 ) )
 2400    CONTINUE
      ENDIF
*
      RETURN
      END