* PACKAGE LEIGEN  !" 固有値解析
*
*"  [HIS] 92/12/14(takepiro) 「行列計算ソフトウェア」使用
*"        92/12/18(takepiro)  Numerical Recipies 使用
*"        93/01/13(takepiro)  4 倍精度計算
*"        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] 
c$$$      COMMON    /COMWRK/ 
*
      DOUBLE PRECISION   WR( MATDMX )             !" 固有値実数部
      DOUBLE PRECISION   WI( MATDMX )             !" 固有値虚数部
      DOUBLE PRECISION   SCALE ( MATDMX )         !" スケーリング情報
      INTEGER            LOW, IGH                 !" 
      DOUBLE PRECISION   INT   ( MATDMX )         !" ヘッセンベルグ変換情報
      DOUBLE PRECISION   A     ( MATDMX, MATDMX ) !" 固有ベクトル
      DOUBLE PRECISION   Z     ( MATDMX, MATDMX ) !" 固有ベクトル
*
      CHARACTER HERR*8                       !" エラー番号出力用文字変数
      INTEGER   NERR                         !" 文字列の長さ
      INTEGER   LENC                         !" 文字の長さを返す関数
      INTEGER            I, J 
*
*"  < 0. 配列をコピーする >
*
      DO 100 J = 1, MATDIM
         DO 100 I = 1, MATDIM
            A( I ,J ) = AMAT( I, J )
 100  CONTINUE 
*
*"  < 1. バランシング >
*
      CALL BALANC
     I      (  MATDMX , MATDIM , 
     M         A      ,
     O         LOW    , IGH    ,SCALE  )
*
*"  < 2. ヘッセンベルグ行列変換 >
*
      CALL ELMHES
     I      (  MATDMX , MATDIM , 
     I         LOW    , IGH    , 
     M         A      , INT      )
*
      CALL ELTRAN
     I      (  MATDMX , MATDIM , 
     I         LOW    , IGH    , 
     I         A      , INT    , 
     O         Z                 )
*
*"  < 3. 固有値, 固有ベクトルの計算 >
*
      CALL HQR2
     I      (  MATDMX , MATDIM , 
     I         LOW    , IGH    , 
     I         A      , 
     O         WR     , WI     , 
     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
     O      (  EIGVLR , EIGVLI , 
     O         EIGVCR , EIGVCI ,
     I         WR     , WI     ,
     I         Z      , 
     I         MATDIM , OEIVEC   )
*
*"  < 5. 終了 > 
*
      RETURN
      END
************************************************************************
      SUBROUTINE EIGSRT      !" 固有値固有ベクトル並べ変え
     O      (  EIGVLR , EIGVLI , 
     O         EIGVCR , EIGVCI , 
     I         WR     , WI     ,
     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]
      DOUBLE PRECISION   WR( MATDMX )             !" 固有値実数部
      DOUBLE PRECISION   WI( MATDMX )             !" 固有値虚数部
      DOUBLE PRECISION   Z ( MATDMX, MATDMX )     !" 固有ベクトル
*
      INTEGER   MATDIM                       !" 計算する次元
      LOGICAL   OEIVEC                       !" 固有関数計算スイッチ
*
*   [INTERNAL WORK]
      INTEGER   J, K
      INTEGER   INDX  ( MATDMX )             !" 順番格納配列
      REAL      WORK  ( MATDMX,MATDMX )      !" 作業領域
*
*" < 0. 固有値の入れ換え >
*
      DO 100 K = 1, MATDIM
         EIGVLR( K ) = WR( K )
         EIGVLI( K ) = WI( K )
 100  CONTINUE 
*
*" < 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( WORK, EIGVLR, MATDIM )             !" 固有値実数部
      DO 2100 J = 1, MATDIM                         !" 固有値実数部の大きい順
         EIGVLR( J ) = WORK( INDX( MATDIM-J+1 ),1  )
 2100 CONTINUE
*
      CALL COPY( WORK, EIGVLI, MATDIM )             !" 固有値虚数部
      DO 2200 J = 1, MATDIM                         !" 固有値実数部の大きい順
         EIGVLI( J ) = WORK( INDX( MATDIM-J+1 ),1 )
 2200 CONTINUE
*
*
*"  < 3. 固有ベクトルのソート >
*
      IF ( OEIVEC )THEN                 
         CALL COPY( WORK, EIGVCR, MATDMX*MATDIM )   !" 固有ベクトル実数部
         DO 2300 K = 1, MATDIM                      !" 固有値実数部の大きい順
            DO 2300 J = 1, MATDIM
               EIGVCR( J,K ) = WORK( J,INDX( MATDIM-K+1 ) )
 2300    CONTINUE
*
         CALL COPY( WORK, EIGVCI, MATDMX*MATDIM )   !" 固有ベクトル虚数部
         DO 2400 K = 1, MATDIM                      !" 固有値実数部の大きい順
            DO 2400 J = 1, MATDIM
               EIGVCI( J,K ) = WORK( J,INDX( MATDIM-K+1 ) )
 2400    CONTINUE
      ENDIF
*
      RETURN
      END
