!***********************************************************************
! ISPACK FORTRAN SUBROUTINE LIBRARY FOR SCIENTIFIC COMPUTING
! Copyright (C) 1998--2016 Keiichi Ishioka <ishioka@gfd-dennou.org>
!
! This library is free software; you can redistribute it and/or
! modify it under the terms of the GNU Lesser General Public
! License as published by the Free Software Foundation; either
! version 2.1 of the License, or (at your option) any later version.
!
! This library is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
! Lesser General Public License for more details.
! 
! You should have received a copy of the GNU Lesser General Public
! License along with this library; if not, write to the Free Software
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
! 02110-1301 USA.
!***********************************************************************
      use wt_mpi_module_supack
      
!      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT NONE
      INCLUDE 'mpif.h'
      INTEGER,PARAMETER :: JM=2**9
!      INTEGER,PARAMETER :: JM=2**10
!      INTEGER,PARAMETER :: MM=JM-1,IM=JM*2,JV=256,NM=MM,NN=NM
      INTEGER,PARAMETER :: MM=JM-1,IM=JM*2,JV=4,NM=MM,NN=NM
      INTEGER,PARAMETER :: KM=64, LM=48
      INTEGER,PARAMETER :: NTR=20
      REAL(8),PARAMETER :: RI=0.5, RO=1.5
      INTEGER :: N,M,L,LR,LI,LAS,ISEED,ICPU,IERR,NP,IP,IPOW,ITR,MAXTD,K
      REAL(8) :: SLMAX,RAN,SLAMAX,SL,GFLOPS,RC,TIM0,TIM1
      
      real(8), allocatable :: wt_data(:,:)
      real(8), allocatable :: wt_data0(:,:)
      real(8), allocatable :: xvz_data(:,:,:)

      REAL(8),DIMENSION(:),ALLOCATABLE :: SALL
      REAL(8),DIMENSION(:),ALLOCATABLE :: SDALL
      REAL(8),DIMENSION(:),ALLOCATABLE :: SBUF
!$    INTEGER :: omp_get_max_threads

      RC=(1D0*5*IM*LOG(1D0*IM)/LOG(2D0)*0.5D0*JM+1D0*(MM+1)*(MM+1)*JM) &
           * ( 1D0*5*KM*LOG(1D0*KM)/LOG(2D0)*0.25D0 )
      ! 逆変換/正変換の1回あたりの演算数概算

      CALL MPI_INIT(IERR)
      CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NP,IERR)
      CALL MPI_COMM_RANK(MPI_COMM_WORLD,IP,IERR)

      call wt_mpi_initial( im, jm, km, mm, lm, ri, ro )

      ALLOCATE(wt_data(nc,0:lm))
      ALLOCATE(wt_data0(nc,0:lm))
      ALLOCATE(xvz_data(0:im-1,jc,0:km))

      IF(IP.EQ.0) THEN
!!$         ALLOCATE(SALL((MM+1)*(MM+1)))
!!$         ALLOCATE(SDALL((MM+1)*(MM+1)))
!!$         ALLOCATE(SBUF((MM/NP+1)*(2*(NN+1)-MM/NP*NP)*NP))
         ISEED=1
         DO K=0,LM
            DO L=1,nc
               CALL MLRAND(ISEED,RAN)
               wt_Data(l,k) = 2*RAN-1
            END DO
         END DO
!!$         DO L=1,(MM+1)*(MM+1)
!!$            CALL MLRAND(ISEED,RAN)
!!$            SALL(L)=2*RAN-1
!!$            SDALL(L)=SALL(L)
!!$         END DO
         PRINT '(A,I5,A,I5,A,I5,A,I5,A,I4)','MM=',MM,', IM=',IM,' JM=',JM,' JV=',JV,', NTR=',NTR
         CALL LVGCPU(ICPU)
         IF(ICPU.EQ.0) THEN
            PRINT '(A)','SSE=fort'
         ELSE IF(ICPU.EQ.10) THEN
            PRINT '(A)','SSE=avx'
         ELSE IF(ICPU.EQ.20) THEN
            PRINT '(A)','SSE=fma'
         ELSE IF(ICPU.EQ.100) THEN
            PRINT '(A)','SSE=fort, for SX'
         END IF

         MAXTD=1
!$       MAXTD=omp_get_max_threads()      
         PRINT '(A,I3)','number of threads =',MAXTD
         PRINT '(A,I3)','number of processes =',NP
      END IF

!!$      CALL SUSS2S(MM,NN,SALL,w_data,SBUF)

      IPOW=0
      CALL MLTIME(TIM0)
      DO ITR=1,NTR
         xvz_data = xvz_wt(wt_data)
      END DO
      CALL MLTIME(TIM1)
      GFLOPS=RC*NTR/(TIM1-TIM0)/1D9
      IF(IP.EQ.0) THEN
        PRINT '(A,ES9.2,A,F6.1,A)','S2G: ',(TIM1-TIM0)/NTR,' sec  (',  &
     &      GFLOPS,' GFlops)'
      END IF

      IPOW=0
      CALL MLTIME(TIM0)
      DO ITR=1,NTR
         wt_data = wt_xvz(xvz_data)
      END DO
      CALL MLTIME(TIM1)
      GFLOPS=RC*NTR/(TIM1-TIM0)/1D9
      IF(IP.EQ.0) THEN
        PRINT '(A,ES9.2,A,F6.1,A)','G2S: ',(TIM1-TIM0)/NTR,' sec  (',  &
     &      GFLOPS,' GFlops)'
      END IF

!!$      CALL SUGS2S(MM,NN,w_data,SDALL,SBUF)
!!$
!!$      IF(IP.EQ.0) THEN
!!$         DEALLOCATE(SBUF)
!!$         SLMAX=0
!!$         SLAMAX=0
!!$         M=0
!!$         DO N=0,MM
!!$            CALL SVNM2L(MM,N,M,L)
!!$            SL=ABS(SDALL(L)-SALL(L))
!!$            IF(SL.GT.SLMAX) THEN
!!$               SLMAX=SL
!!$               LAS=L
!!$            END IF
!!$            SLAMAX=SLAMAX+SL**2
!!$            !        print *,N,M,SL
!!$         END DO
!!$
!!$         DO M=1,MM
!!$            DO N=M,MM
!!$               CALL SVNM2L(MM,N,M,LR)
!!$               CALL SVNM2L(MM,N,-M,LI)
!!$               SL=(SDALL(LR)-SALL(LR))**2+(SDALL(LI)-SALL(LI))**2
!!$               SL=SQRT(SL)
!!$               IF(SL.GT.SLMAX) THEN
!!$                  SLMAX=SL
!!$                  LAS=LR
!!$               END IF
!!$               !          print *,N,M,SL
!!$               SLAMAX=SLAMAX+SL**2
!!$            END DO
!!$         END DO
!!$         CALL SVL2NM(MM,LAS,N,M)
!!$         PRINT '(A,ES9.2,A,I5,A,I5,A)','maxerror =',SLMAX,' (n=',N,', m=',M,')'
!!$         PRINT '(A,ES9.2)','rmserror =',SQRT(SLAMAX/((MM+1)*(MM+2)/2))
!!$         DEALLOCATE(SALL)
!!$         DEALLOCATE(SDALL)
!!$      END IF

      CALL MPI_FINALIZE(IERR)

      END