* PACKAGE LGLU1  !" $B@~7?2r@O(B : $B8GM-%Y%/%H%kJ,2r$K$h$kO"N)J}Dx<0$N2r(B
*
*"  [HIS] 93/07/14(takepiro) $B!V9TNs7W;;%=%U%H%&%'%"!W;HMQ(B
*"                             $B8GM-CM$N@dBPCM$,(B 0 $B$N%b!<%I$O7W;;$7$J$$(B
*"        93/08/11(takepiro) 0 $B8GM-CM%b!<%INe5/$NI>2A$r=$@5(B
*
**********************************************************************
      SUBROUTINE LEIGLU
     M            ( AMAT   , 
     M              BVEC   ,
     O              IERR   , 
     I              MATDIM , OEIGEN  )
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 
      INCLUDE   (ZLDIM)                      !" $B3J;RE@?t(B, $BGH?t(B
#else
#include        "zcdim.F"                    !" 
#include        "zldim.F"                    !" $B78?t9TNs$NBg$-$5(B
#endif
*
*   [MODIFY] 
      REAL      AMAT  ( MATDMX, MATDMX )     !" $B78?t9TNs(B
      REAL      BVEC  ( MATDMX )             !" $B1&JU%Y%/%H%k(B / $B2r(B
*
*   [OUTPUT] 
      INTEGER   IERR                         !" $B%5%V%k!<%A%s%(%i!<%3!<%I(B
*
*   [INPUT] 
      INTEGER   MATDIM                       !" $B7W;;$9$k<!85(B
      LOGICAL   OEIGEN                       !" $B8GM-%Y%/%H%k7W;;%9%$%C%A(B
*
*   [INTERNAL WORK] 
      REAL      EIGVLR( MATDMX )             !" $B8GM-CM<B?tIt(B
      REAL      EIGVLI( MATDMX )             !" $B8GM-CM5u?tIt(B
      REAL      EIGVCR( MATDMX, MATDMX )     !" $B8GM-%Y%/%H%k<B?tIt(B
      REAL      EIGVCI( MATDMX, MATDMX )     !" $B8GM-%Y%/%H%k5u?tIt(B
*
      INTEGER   IP    ( MATDMX )             !" $B%T%\%C%H>pJs(B
      REAL      WORK  ( MATDMX )             !" $B:n6HG[Ns(B
      REAL      EPS                          !" $BFC0[@-$NH=Dj4p=`(B
      DATA      EPS  / 3.52D-15 /
*
      REAL      ABSV2 , ABSB2
      INTEGER   IUNIT, I, J
      LOGICAL   OFIRST 
      DATA      OFIRST  / .TRUE. /
*
      SAVE      EPS, OFIRST
*
      IF ( OFIRST )THEN
         OFIRST = .FALSE.
         CALL GLPGET( 'MSGUNIT', IUNIT )
         WRITE ( IUNIT, * ) ' EXPAND BY EIGENVECTOR : 93/09/09'
      ENDIF
*
*"  < 1. $B8GM-%Y%/%H%k$r5a$a$k(B >
*
      IF ( OEIGEN )THEN
*
         CALL LEIGEN
     M        ( AMAT   , 
     O          EIGVLR , EIGVLI  ,
     O          EIGVCR , EIGVCI  ,
     O          IERR   , 
     I          MATDIM , .TRUE.    )
*
*"  < 2. $B8GM-%Y%/%H%k$G1&JU$rE83+(B >
*
         ABSB2 = 0.
         DO 50 J = 1, MATDIM
            ABSB2 = ABSB2 + BVEC(J)**2
  50     CONTINUE 
*
         DO 100 J = 1, MATDIM
            IF ( EIGVLI(J) .EQ. 0 )THEN
               DO 110 I = 1, MATDIM
                  AMAT ( I,J ) = EIGVCR( I,J )
 110           CONTINUE 
            ELSE IF ( EIGVLI(J) .GT. 0 )THEN
               DO 120 I = 1, MATDIM
                  AMAT ( I,J ) = 2 * EIGVCR( I,J )
 120           CONTINUE 
            ELSE
               DO 130 I = 1, MATDIM
                  AMAT ( I,J ) = - 2 * EIGVCI( I,J )
 130           CONTINUE 
            ENDIF
 100     CONTINUE 
*
         CALL GLU1
     M          ( AMAT     , 
     I            MATDIM   , MATDMX , 
     O            EPS      , 
     W            WORK     , 
     O            IP       , IERR     )
      ENDIF
*
      IF ( IERR .EQ. 0 ) THEN
         CALL GSLV1
     I          ( AMAT     , 
     I            MATDIM   , MATDMX , 
     M            BVEC     , 
     I            IP                )
      ELSE
         CALL MSGDMP
     &         ( 'W', 'LGLU', 'SOL. NOT CALCULATED BECAUSE OF ERRORS' )
      END IF
*
*"  < 3. $B8GM-%Y%/%H%k$N=E$M9g$o$;$G2r$r5a$a$k(B >
*
      CALL RESET( WORK , MATDIM )
*
*
      DO 300 J = 1, MATDIM
         ABSV2 = EIGVLR(J)**2 + EIGVLI(J)**2
         IF( ABSV2 .GT. EPS )THEN
            IF( EIGVLI(J) .EQ. 0.0 )THEN
               DO 310 I = 1, MATDIM
                  WORK(I) = WORK(I) 
     &                 + BVEC(J)*EIGVCR(I,J) / EIGVLR(J)
 310           CONTINUE 
            ELSE IF( EIGVLI(J) .GT. 0.0 )THEN
               DO 320 I = 1, MATDIM
                  WORK(I) = WORK(I) 
     &               + 2 *(   EIGVLR(J) * EIGVCR(I,J)
     &                      + EIGVLI(J) * EIGVCI(I,J) ) * BVEC(J)/ABSV2
 320           CONTINUE 
            ELSE 
               DO 330 I = 1, MATDIM
                  WORK(I) = WORK(I) 
     &               + 2 *( - EIGVLR(J) * EIGVCI(I,J)
     &                      + EIGVLI(J) * EIGVCR(I,J) ) * BVEC(J)/ABSV2
 330           CONTINUE 
            ENDIF
         ELSE
            CALL MSGDMP( 'W','LEIGLU','0-EIGENVALUE MODE NEGLECTED' )
            IF ( BVEC(J) .NE. 0.0 )THEN
               CALL MSGDMP( 'W','LEIGLU','0-EIGENVALUE MODE EXCITED' )
               ABSV2 = 0.
               DO 340 I = 1, MATDIM
                  ABSV2 = ABSV2 + (BVEC(J)*EIGVCR(I,J))**2
 340           CONTINUE 
               WRITE ( 6, * ) SQRT(ABSV2)/SQRT(ABSB2)
            ENDIF
         ENDIF
 300  CONTINUE 
*
      CALL COPY( BVEC , WORK, MATDIM )
*
      RETURN
      END
