!***********************************************************************
! 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 w_module_svpack

      IMPLICIT NONE
      INTEGER,PARAMETER :: NTR=50
      INTEGER,PARAMETER :: JM=2**10
!      INTEGER,PARAMETER :: JM=2**11
      INTEGER,PARAMETER :: IM=2*JM,MM=JM-1
      INTEGER,PARAMETER :: NM=MM,NN=MM
      INTEGER :: ISEED,L,N,M,LR,LI,LAS,ITR,MAXTD,ICPU
      REAL(8) :: RC,RAN,SL,SLMAX,SLAMAX,GFLOPS,TIM0,TIM1
      real(8), allocatable :: w_data(:)
      real(8), allocatable :: w_data0(:)
      real(8), allocatable :: xy_data(:,:)
!$    INTEGER :: omp_get_max_threads

      ALLOCATE(w_data((mm+1)*(mm+1)))
      ALLOCATE(w_data0((mm+1)*(mm+1)))
      ALLOCATE(xy_data(0:im-1,jm))

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

      call w_initial( mm, im, jm )

      ISEED=1
      DO L=1,(MM+1)*(MM+1)
         CALL MLRAND(ISEED,RAN)
         w_data(l)=2*RAN-1
      END DO
      w_data0 = w_data

      PRINT '(A,I5,A,I5,A,I5,A,I4)','MM=',MM,', IM=',IM,' JM=',JM,', 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

      CALL MLTIME(TIM0)
      DO ITR=1,NTR
         xy_data = xy_w(w_data)
      END DO
      CALL MLTIME(TIM1)
      GFLOPS=RC*NTR/(TIM1-TIM0)/1D9
      PRINT '(A,ES9.2,A,F6.1,A)','S2G: ',(TIM1-TIM0)/NTR,' sec  (', &
     & GFLOPS,' GFlops)'

      CALL MLTIME(TIM0)
      DO ITR=1,NTR
         w_data = w_xy(xy_data)
      END DO
      CALL MLTIME(TIM1)
      GFLOPS=RC*NTR/(TIM1-TIM0)/1D9
      PRINT '(A,ES9.2,A,F6.1,A)','G2S: ',(TIM1-TIM0)/NTR,' sec  (', &
     & GFLOPS,' GFlops)'

      SLMAX=0
      SLAMAX=0
      M=0
      DO N=0,MM
         CALL SVNM2L(MM,N,M,L)
         SL=ABS(w_data0(L)-w_data(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=(w_data0(LR)-w_data(LR))**2+(w_data0(LI)-w_data(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(w_data)
      DEALLOCATE(w_data0)
      DEALLOCATE(xy_data)
        
      END