C*********************************************************************
C*********************************************************************
C
      SUBROUTINE PSETUP
     B   (  DENS0,
     G      FDZ  , FDZM,
     M      VMODE, TMODE, VEIGEN, FMODE, TRIGS, IFAC )
C
C*********************************************************************
      IMPLICIT REAL*8 ( A-H, O-Z )
C*********************************************************************
C
      include 'grid_size.f'
C
      dimension trigs(nx+15), ifac(15)
C
      DIMENSION
     \   DENS0(-2:NZB), FDZ(-2:NZB), FDZM(-2:NZB),
     \   VMODE( -2:NZB, NZ ), TMODE( -2:NZB, NZ ),
     \   VEIGEN ( NZ ),
     \   FMODE( -2:NXB, NZ )
C
      DIMENSION
     \   A(NZ, NZ), VR(NZ,NZ), VI(NZ,NZ),
     \   ER(NZ), EI(NZ)
c
      DOUBLE PRECISION   RCONDE( NZ ), RCONDV( NZ ),
     $                   SCALE( NZ ), 
     $                   WRK( NZ*(NZ+6) )
c
      dimension unit(nz,nz)

*     .. Array Arguments ..
      INTEGER            IPIV( NZ ), IWORK( NZ )
      DOUBLE PRECISION   AF( NZ, NZ ), 
     $                   BERR( NZ ), C( NZ ), FERR( NZ ), R( NZ ),
     $                   WORK( 4*NZ ), X( NZ, NZ )
*     ..
C
C=====================================================================
C
      DO 10 I = 1, NZ
      DO 10 J = 1, NZ
         A(I,J) = 0.0D0
   10 CONTINUE
C
      DO 11 IZ = 0, NZ-2
         A(IZ+1,IZ+2) =
     \       + 1.125D0
     \          * ( DENS0(IZ+1) + DENS0(IZ) ) / 2.0D0 / FDZM(IZ+1)
     \          / DENS0(IZ) / FDZ(IZ)
   11 CONTINUE
C
      DO 12 IZ = 1, NZ-1
         A(IZ+1,IZ) =
     \       + 1.125D0
     \          * ( DENS0(IZ-1) + DENS0(IZ) ) / 2.0D0 / FDZM(IZ)
     \          / DENS0(IZ) / FDZ(IZ)
   12 CONTINUE
C
      DO 13 IZ = 0, NZ-2
         A(IZ+1,IZ+1) = A(IZ+1,IZ+1) - A(IZ+1,IZ+2)
   13 CONTINUE
C
      DO 14 IZ = 1, NZ-1
         A(IZ+1,IZ+1) = A(IZ+1,IZ+1) - A(IZ+1,IZ)
   14 CONTINUE
C
      A(1 ,1   ) = 26.0D0 / 27.0D0 * A( 1,1   )
      A(1 ,2   ) = 26.0D0 / 27.0D0 * A( 1,2   )
      A(NZ,NZ  ) = 26.0D0 / 27.0D0 * A(NZ,NZ  )
      A(NZ,NZ-1) = 26.0D0 / 27.0D0 * A(NZ,NZ-1)
C
      DO 21 IZ = 0, NZ-3
         A(IZ+1,IZ+3) =
     \       - 1.0D0 / 24.0D0
     \          * ( DENS0(IZ+2) + DENS0(IZ+1) ) / 2.0D0 / FDZM(IZ+2)
     \          / DENS0(IZ) / FDZ(IZ)
   21 CONTINUE
C
      DO 22 IZ = 2, NZ-1
         A(IZ+1,IZ-1) =
     \       - 1.0D0 / 24.0D0
     \          * ( DENS0(IZ-2) + DENS0(IZ-1) ) / 2.0D0 / FDZM(IZ-1)
     \          / DENS0(IZ) / FDZ(IZ)
   22 CONTINUE
C
      DO 23 IZ = 0, NZ-3
         A(IZ+1,IZ+2) = A(IZ+1,IZ+2) - A(IZ+1,IZ+3)
   23 CONTINUE
C
      DO 24 IZ = 2, NZ-1
         A(IZ+1,IZ  ) = A(IZ+1,IZ  ) - A(IZ+1,IZ-1)
   24 CONTINUE
C
C *** EIGENVALUE PROBLEM 
C
      LWORK = NZ*(NZ+6)

      call DGEEVX( 'B', 'V', 'V', 'B', NZ, A, NZ, ER, EI,
     $                   VI, NZ, VR, NZ, ILO, IHI, SCALE, ABNRM,
     $                   RCONDE, RCONDV, WRK, LWORK, IWORK, INFO )

*      write(0,*) 'eigen info = ', info
*      write(0,*) 'eigen next time you should set LWORK =  ', wrk(1)
C
C
*
*      DO  IM = 1, NZ
*         write(*,*) '----------------------------------------'
*         write(*,*) 'im = ', im
*         write(*,*) 'cond. for eigenvalue = ', 1.0d0 / rconde(im)
*         write(*,*) 'cond. for eigenvector = ', 1.0d0 / rcondv(im)
*         WRITE(6,*)'ER'
*         WRITE(6,*) ER(IM)
*         WRITE(6,*)'VR'
*         WRITE(6,*) (VR(IZ,IM), IZ=1, NZ-1)
*      end do
*
C
C   * ZERO-EIGENVALUE *
C
      IZERO = 0
      EPS = 1.0D-10
C
      DO 101 I = 1, NZ
         IF ( ABS( ER(I) ) .LT. EPS ) IZERO = I
  101 CONTINUE


      DO 121 IM = 1, IZERO - 1
         VEIGEN(IM) = ER(IM)
         DO 122 IZ = 0, NZ-1
            VMODE(IZ,IM) = VR(IZ+1,IM)
  122    CONTINUE
  121 CONTINUE
C
      DO 131 IM = IZERO+1, NZ
         VEIGEN(IM-1) = ER(IM)
         DO 132 IZ = 0, NZ-1
            VMODE(IZ,IM-1) = VR(IZ+1,IM)
  132    CONTINUE
  131 CONTINUE
C
      VEIGEN ( NZ ) = 0.0D0
      DO 141 IZ = 0, NZ-1
         VMODE (IZ,NZ) = SQRT( 1.0D0 / NZ )
  141 CONTINUE
C
C ***  MATRIX INVERSION using DGESVX in netlib *** 95/11/29
C
      DO 151 I = 1, NZ
      DO 151 J = 1, NZ
         A (I,J) = VMODE( I-1, J )
  151 CONTINUE
C
*
      do i = 1, nz
         unit(i,i) = 1.0d0
      end do
*     
      call DGESVX( 'E', 'N', NZ, NZ, A, NZ, AF, NZ, IPIV,
     $             EQUED, R, C, unit, NZ, X, NZ, RCOND, FERR, BERR,
     $             WORK, IWORK, INFO )
*
*      write(0,*) 'info = ', INFO
*      write(0,*) 'detA/detU = ', work(1)
*
*      write(0,*) 'cond = ', 1.0d0 / RCOND
*      write(0,*) 'ferr = ', FERR
*      write(0,*) 'berr = ', BERR
*
*
      DO 201 IM = 1, NZ
      DO 201 IZ = 0, NZ-1
         TMODE (IZ,IM) = X ( IM, IZ+1 )
  201 CONTINUE
C
C *** 圧力方程式用の係数行列　***
C
      DO 301 IM = 1, NZ-1
         FMODE (1 ,IM) = 1.0D0 /   VEIGEN(IM)
         FMODE (NX,IM) = 1.0D0 / ( VEIGEN(IM) - 14.0D0 / 3.0D0 )
  301 CONTINUE
C
      FMODE(1 ,NZ) = 0.0D0
      FMODE(NX,NZ) = - 3.0D0 / 14.0D0
C
      PAI = 4.0D0 * ATAN( 1.0D0 )
C
      DO 302 IM = 1, NZ
      DO 302 KX = 1, NX/2-1
         FMODE ( 2*KX  , IM )
     \     = 1.0D0 / (     VEIGEN(IM)
     \        + (   28.0D0 * COS ( 2.0D0 * PAI * KX / NX )
     \            -          COS ( 4.0D0 * PAI * KX / NX )
     \            - 27.0D0  ) / 12.0D0                        )
         FMODE ( 2*KX+1, IM ) = FMODE ( 2*KX, IM )
  302 CONTINUE
C
C *** 高速フーリエ変換ルーチンの初期設定　*** FTLIB.FORT(FFTEC) ***
C
c      CALL RFFTIM ( NX, TRIGS, IFAC )
c
      call vrffti ( nx, trigs )
C
C
C====================================================================
C
      RETURN
      E N D
