!!  Copyright (C) Stichting Deltares, 2005-2023.
!!
!!  This file is part of iMOD.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License as published by
!!  the Free Software Foundation, either version 3 of the License, or
!!  (at your option) any later version.
!!
!!  This program 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 General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!
!!  Contact: imod.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands.
!!
MODULE MOD_LIS
 USE IMODVAR, ONLY : DP_KIND,SP_KIND     
 USE MOD_UTL, ONLY : ITOS
 
 PUBLIC :: LIS_CREATE,LIS_DESTROY,LIS_INIT,LIS_SOLVE
 INTEGER, PRIVATE :: IERR
 
 CONTAINS

 !###======================================================================
 SUBROUTINE LIS_INIT() ! ISLV,KSP,MAT,RHS_VEC,SOL_VEC,INDX )
 !###======================================================================
 IMPLICIT NONE
 
 CALL LIS_FINITIALIZE_F(IERR)
 
 END SUBROUTINE LIS_INIT
      
 !###======================================================================
 SUBROUTINE LIS_CREATE() ! ISLV,KSP,MAT,RHS_VEC,SOL_VEC,INDX )
 !###======================================================================
 IMPLICIT NONE
!---  Create Jacobian matrix  ---
!
!      CALL lis_matrix_create_f(1,MAT,IERR)
!      CALL lis_matrix_set_size_f(MAT,NLX,NGX,IERR)
!
!---  Create solution vector  ---
!
!      CALL lis_vector_duplicate_f(MAT,SOL_VEC,IERR)
!
!---  Create problem (RHS) vector  ---
!
!      CALL lis_vector_duplicate_f(MAT,RHS_VEC,IERR)
!
!---  Create solver  ---
!
!      CALL lis_solver_create_f(KSP,IERR)
!
!---  Set the options for the solver  ---
!
!      CALL lis_solver_set_option_f('-i bicgstab',KSP,IERR)
!      CALL lis_solver_set_option_f('-p ilu -ilu_fill 0',KSP,IERR)
!      CALL lis_solver_set_option_f('-print none -adds true',KSP,IERR)
!      CALL lis_solver_set_optionc_f(KSP,IERR)
!
 END SUBROUTINE LIS_CREATE

 !###======================================================================
 SUBROUTINE LIS_DESTROY() ! KSP,MAT,RHS_VEC,SOL_VEC )
 !###======================================================================
 IMPLICIT NONE
!!---  Unset Jacobian matrix from DLU array  ---
!!
!      CALL lis_matrix_unset_f(MAT,IERR)
!!
!!---  Destroy Jacobian matrix  ---
!!
!      CALL lis_matrix_destroy_f(MAT,IERR)
!!
!!---  Destroy solution vector  ---
!!
!      CALL lis_vector_destroy_f(SOL_VEC,IERR)
!!
!!---  Destroy problem vector  ---
!!
!      CALL lis_vector_destroy_f(RHS_VEC,IERR)
!!
!!---  Destroy solver  ---
!!
!      CALL lis_solver_destroy_f(KSP,IERR)
!!
  CALL LIS_FINALIZE_F(IERR)
  
 END SUBROUTINE LIS_DESTROY

 !###======================================================================
 SUBROUTINE LIS_SOLVE() ! ISLV,KSP,MAT,RHS_VEC,SOL_VEC,INDX )
 !###======================================================================
 IMPLICIT NONE
!!---  Load Lis problem vector  ---
!!
!      CALL lis_vector_set_values_f(0,NGX,ILU,BLU,RHS_VEC,IERR)
!!
!!---  Load nonzero elements of the Lis matrix for transport  ---
!!
!      IF( ISLV.EQ.-1 ) THEN
!        NNZ = NLUC(NGX+1)
!        CALL lis_matrix_set_csr_f(NNZ,NLUC,MLUC,DLU,MAT,IERR)
!!
!!---  Load nonzero elements of the Lis matrix for geomechanics  ---
!!
!      ELSEIF( INDX.EQ.2 ) THEN
!        NNZ = NLU_GM(NGX+1)
!        CALL lis_matrix_set_csr_f(NNZ,NLU_GM,MLU_GM,DLU,MAT,IERR)
!!
!!---  Load nonzero elements of the Lis matrix for coupled-flow  ---
!!
!      ELSE
!        NNZ = NLU(NGX+1)
!        CALL lis_matrix_set_csr_f(NNZ,NLU,MLU,DLU,MAT,IERR)
!      ENDIF
!      CALL lis_matrix_assemble_f(MAT,IERR)
!!
!!---  Output coupled-flow problem vector  ---
!!
!      IF( ISLC(34).EQ.1 ) THEN
!        CALL lis_output_vector_f(RHS_VEC,2,"rhs_vec.dat",IERR)
!        STOP
!      ENDIF
!!
!!---  Output coupled-flow matrix  ---
!!
!      IF( ISLC(34).EQ.2 ) THEN
!        CALL lis_output_matrix_f(MAT,2,"mat.dat",IERR)
!        STOP
!      ENDIF
!!
!!---  Output transport problem vector  ---
!!
!      IF( ISLC(34).EQ.11 .AND. ISLV.EQ.-1 ) THEN
!        CALL lis_output_vector_f(RHS_VEC,2,"rhs_vec_c.dat",IERR)
!        STOP
!      ENDIF
!!
!!---  Output transport matrix  ---
!!
!      IF( ISLC(34).EQ.12 .AND. ISLV.EQ.-1 ) THEN
!        CALL lis_output_matrix_f(MAT,2,"mat_c.dat",IERR)
!        STOP
!      ENDIF
!!
!!---  Output geomechanics problem vector  ---
!!
!      IF( ISLC(34).EQ.21 .AND. INDX.EQ.2 ) THEN
!        CALL lis_output_vector_f(RHS_VEC,2,"rhs_vec_g.dat",IERR)
!        STOP
!      ENDIF
!!
!!---  Output geomechanics matrix  ---
!!
!      IF( ISLC(34).EQ.22 .AND. INDX.EQ.2 ) THEN
!        CALL lis_output_matrix_f(MAT,2,"mat_g.dat",IERR)
!        STOP
!      ENDIF
!!
!!---  Solve  ---
!!
!      CALL lis_solve_f(MAT,RHS_VEC,SOL_VEC,KSP,IERR)
!!
!!---  Output coupled-flow solution vector  ---
!!
!      IF( ISLC(34).EQ.3 ) THEN
!        CALL lis_output_vector_f(SOL_VEC,2,"fsol_vec.dat",IERR)
!        STOP
!      ENDIF
!!
!!---  Output transport solution vector  ---
!!
!      IF( ISLC(34).EQ.13 .AND. ISLV.EQ.-1 ) THEN
!        CALL lis_output_vector_f(SOL_VEC,2,"fsol_vec_c.dat",IERR)
!        STOP
!      ENDIF
!!
!!---  Output geomechanics solution vector  ---
!!
!      IF( ISLC(34).EQ.23 .AND. INDX.EQ.2 ) THEN
!        CALL lis_output_vector_f(SOL_VEC,2,"fsol_vec_g.dat",IERR)
!        STOP
!      ENDIF
!!
!!---  Get solution  ---
!!
!      ISX = 1
!      ICX = NGX
!      CALL lis_vector_get_values_f(SOL_VEC,ISX,ICX,BLU,IERR)

 END SUBROUTINE LIS_SOLVE

 !###======================================================================
 SUBROUTINE LIS_TEST()
 !###======================================================================
 IMPLICIT NONE
 
 CALL LIS_SOLVE_MATRIX() 
 CALL LIS_SOLVE_EIGENVALUE()
 
 END SUBROUTINE LIS_TEST
 
 !###======================================================================
 SUBROUTINE LIS_SOLVE_MATRIX()
 !###======================================================================
 IMPLICIT NONE
 INTEGER,PARAMETER :: LIS_INS_VALUE=0
 INTEGER,PARAMETER :: LIS_ADD_VALUE=1
 INTEGER,PARAMETER :: LIS_SUB_VALUE=2
 INTEGER,PARAMETER :: LIS_MATRIX_CSR=1
 INTEGER :: I,GN,IS,IE,ITER,IERR,N
 INTEGER(KIND=DP_KIND) :: A
 INTEGER(KIND=DP_KIND) :: B,X,U
 INTEGER(KIND=DP_KIND) :: SOLVER
 REAL(KIND=DP_KIND) :: TIME,RES
 REAL(KIND=DP_KIND),DIMENSION(12) :: RESULTS

 N=12
 CALL LIS_FINITIALIZE_F(IERR)
 CALL LIS_MATRIX_CREATE_F(1,A,IERR)
 CALL LIS_MATRIX_SET_SIZE_F(A,0,N,IERR)
 CALL LIS_MATRIX_GET_SIZE_F(A,N,GN,IERR)
 CALL LIS_MATRIX_GET_RANGE_F(A,IS,IE,IERR)
 DO I=IS,IE-1
 IF(I.GT.1) CALL LIS_MATRIX_SET_VALUE_F(LIS_INS_VALUE,I,I-1,-1.0D0,A,IERR)
 IF(I.LT.GN)CALL LIS_MATRIX_SET_VALUE_F(LIS_INS_VALUE,I,I+1,-1.0D0,A,IERR)
 CALL LIS_MATRIX_SET_VALUE_F(LIS_INS_VALUE,I,I,2.0D0,A,IERR)
 ENDDO
 !## set type of matrix
 CALL LIS_MATRIX_SET_TYPE_F(A,LIS_MATRIX_CSR,IERR)
 !## assemble matrix
 CALL LIS_MATRIX_ASSEMBLE_F(A,IERR)

 !## create vector with dimension of A
 CALL LIS_VECTOR_DUPLICATE_F(A,U,IERR)
 CALL LIS_VECTOR_DUPLICATE_F(A,B,IERR)
 CALL LIS_VECTOR_DUPLICATE_F(A,X,IERR)
 !## set all values to one
 CALL LIS_VECTOR_SET_ALL_F(1.0D0,U,IERR)
 !## multiply Au=b
 CALL LIS_MATVEC_F(A,U,B,IERR)

 CALL LIS_SOLVER_CREATE_F(SOLVER,IERR)
 CALL LIS_SOLVER_SET_OPTIONC_F(SOLVER,IERR)
 !## A=system matrix; B=RHS; x=solution
 CALL LIS_SOLVE_F(A,B,X,SOLVER,IERR)
 !## get solve-information
 CALL LIS_SOLVER_GET_ITER_F(SOLVER,ITER,IERR)
 WRITE(*,*) 'NUMBER OF ITERATIONS = ',ITER
 CALL LIS_SOLVER_GET_TIME_F(SOLVER,TIME,IERR)
 WRITE(*,*) 'SOLUTION TIME (SECS) = ',TIME
 CALL LIS_SOLVER_GET_RESIDUALNORM_F(SOLVER,RES,IERR)
 WRITE(*,*) 'RESIDUAL (SECS) = ',RES
 CALL LIS_VECTOR_PRINT_F(X,IERR)
 CALL LIS_VECTOR_GET_VALUES_F(X,1,N,RESULTS,IERR)

 CALL LIS_MATRIX_DESTROY_F(A,IERR)
 CALL LIS_VECTOR_DESTROY_F(B,IERR)
 CALL LIS_VECTOR_DESTROY_F(X,IERR)
 CALL LIS_VECTOR_DESTROY_F(U,IERR)
 CALL LIS_SOLVER_DESTROY_F(SOLVER,IERR)
 CALL LIS_FINALIZE_F(IERR)

 END SUBROUTINE LIS_SOLVE_MATRIX
 
 !###======================================================================
 SUBROUTINE LIS_SOLVE_EIGENVALUE()
 !###======================================================================
 IMPLICIT NONE
 INTEGER,PARAMETER :: N=3
 REAL(KIND=DP_KIND),DIMENSION(N,N) :: SYSA
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:) :: EIGV
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:) :: EIGW
 REAL(KIND=DP_KIND) :: ESTOP=100.0D0 !99.9
 INTEGER :: I,J
 
 DATA SYSA/2.0D0,0.0D0,0.0D0, &
           0.0D0,3.0D0,4.0D0, &
           0.0D0,4.0D0,9.0D0/
 
 CALL LIS_EIGENVALUE(SYSA,EIGV,EIGW,ESTOP)
 
 DO I=1,N
  WRITE(*,'(A,F17.5)') 'EIGENVALUE (%)'//TRIM(VTOS(I)),EIGW(I)
 ENDDO
 DO J=1,N
  WRITE(*,'(A,99F17.5)') 'EIGENVECTOR '//TRIM(VTOS(I)),(EIGV(I,J),I=1,N)
 ENDDO
 
 END SUBROUTINE LIS_SOLVE_EIGENVALUE

 !###======================================================================
 SUBROUTINE LIS_EIGENVALUE(SYSA,EIGV,EIGW,ESTOP)
 !###======================================================================
 IMPLICIT NONE
 REAL(KIND=DP_KIND),DIMENSION(:,:),INTENT(IN) :: SYSA
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:,:),INTENT(INOUT) :: EIGV
 REAL(KIND=DP_KIND),POINTER,DIMENSION(:),INTENT(INOUT) :: EIGW
 REAL(KIND=DP_KIND),INTENT(IN) :: ESTOP
 REAL(KIND=DP_KIND),DIMENSION(:),ALLOCATABLE :: EIGVALUES,EIGERROR,TMP
 INTEGER,PARAMETER :: LIS_FMT_PLAIN=1
 INTEGER,PARAMETER :: LIS_FMT_MM   =2
 INTEGER,PARAMETER :: LIS_INS_VALUE=0
 INTEGER,PARAMETER :: LIS_MATRIX_COO=10
 INTEGER :: I,J,N,GN,IS,IE,ITER,LN,IERR,NSOL,NN
 INTEGER(KIND=DP_KIND) :: A,X,Y,M,V
 INTEGER(KIND=DP_KIND) :: ESOLVER
 REAL(KIND=DP_KIND) :: EVALUE,TIME,TOT
 CHARACTER(LEN=20) :: ESOLVERNAME
 CHARACTER(LEN=50) :: STRING
  
 CALL LIS_FINITIALIZE_F(IERR)

 N=SIZE(SYSA,1) !3
 CALL LIS_MATRIX_CREATE_F(0,A,IERR)
 !## creates matrix A
 CALL LIS_MATRIX_SET_SIZE_F(A,N,N,IERR)
 !## gets the size of matrix A, ln=local n, gn=global n rows
 CALL LIS_MATRIX_GET_SIZE_F(A,LN,GN,IERR)
 !## get the location of matrix A where partial matrix A starts in global matrix, is=1; ie=n+1
 CALL LIS_MATRIX_GET_RANGE_F(A,IS,IE,IERR)
 !## set type of matrix
 CALL LIS_MATRIX_SET_TYPE_F(A,LIS_MATRIX_COO,IERR)
 NN=0; DO I=IS,IE-1
  DO J=IS,IE-1
   !## insert values in matrix A (only if no-equal to zero)
   IF(SYSA(I,J).NE.0.0D0)THEN
    CALL LIS_MATRIX_SET_VALUE_F(LIS_INS_VALUE,I,J,SYSA(I,J),A,IERR)
    NN=NN+1
   ENDIF
  ENDDO
 ENDDO
 WRITE(*,*) 'inserted ',NN,'elements'
 
 !## assemble matrix
 CALL LIS_MATRIX_ASSEMBLE_F(A,IERR)
 
 CALL LIS_VECTOR_DUPLICATE_F(A,X,IERR) 
 CALL LIS_VECTOR_SET_ALL_F(1.0D0,X,IERR) 
 
! CALL LIS_OUTPUT_MATRIX_F(A,LIS_FMT_MM,'D:\Aa.TXT',IERR);
! !CALL CHKERR(IERR)
! CALL LIS_OUTPUT_VECTOR_F(X,LIS_FMT_MM,'D:\X.TXT',IERR);

! CALL LIS_MATRIX_SET_TYPE_F(A,LIS_MATRIX_CSR,IERR)
! CALL LIS_MATRIX_ASSEMBLE_F(A,IERR) 
 CALL LIS_ESOLVER_CREATE_F(ESOLVER,IERR)
 !## save the residual error
 CALL LIS_ESOLVER_SET_OPTION_F('-EPRINT MEM',ESOLVER,IERR)
 !## retrieve nunber of eigenvalues (3)
 IF(N.GT.10)THEN
  WRITE(STRING,'(A)') '-E SI -SS 10 -P NONE -PRINT ALL'  !## subspace
!  WRITE(STRING,'(A)') '-E LI -SS 10 -P NONE'  !## lanczos
!  WRITE(STRING,'(A)') '-E AI -SS 10 -P NONE'  !## arnoldi
 ELSE
  WRITE(STRING,'(A)') '-E SI -SS '//TRIM(VTOS(N))//' -P NONE -PRINT ALL'
 ENDIF
 CALL LIS_ESOLVER_SET_OPTION_F(TRIM(STRING),ESOLVER,IERR)
! CALL LIS_ESOLVER_SET_OPTION_F('-E SI -SS 3 -P NONE',ESOLVER,IERR)
! CALL LIS_ESOLVER_SET_OPTION_F('-ETOL 1.0E-12',ESOLVER,IERR) ! -TOL 1.0E-12',ESOLVER,IERR) 
! CALL LIS_SOLVER_SET_OPTION_F('-I BICGSTAB -TOL 1.0E-12',SOLVER,IERR) 
 
!      CALL lis_solver_set_option_f('-i bicgstab',KSP,IERR)
!      CALL lis_solver_set_option_f('-p ilu -ilu_fill 0',KSP,IERR)
!      CALL lis_solver_set_option_f('-print none -adds true',KSP,IERR)
 
! CALL LIS_ESOLVER_SET_OPTIONC_F(ESOLVER,IERR) 
 !## solves the eigenvalue problem; a=matrix; x=eigenvalues
 CALL LIS_ESOLVE_F(A,X,EVALUE,ESOLVER,IERR)

! CALL LIS_OUTPUT_MATRIX_F(A,LIS_FMT_MM,'D:\Aa2.TXT',IERR);
! CALL LIS_OUTPUT_VECTOR_F(X,LIS_FMT_MM,'D:\X2.TXT',IERR);

 ALLOCATE(EIGVALUES(N),EIGERROR(N))

 CALL LIS_ESOLVER_GET_ITER_F(ESOLVER,ITER,IERR)
 CALL LIS_ESOLVER_GET_TIME_F(ESOLVER,TIME,IERR)
 CALL LIS_ESOLVER_GET_RESIDUALNORM_F(ESOLVER,EIGERROR,IERR)
 CALL LIS_ESOLVER_GET_ESOLVER_F(ESOLVER,NSOL,IERR)
 CALL LIS_ESOLVER_GET_ESOLVERNAME_F(NSOL,ESOLVERNAME,IERR) 
 
 !## get the eigenvectors
 CALL LIS_MATRIX_CREATE_F(0,M,IERR)
 CALL LIS_MATRIX_SET_SIZE_F(M,N,N,IERR)
 CALL LIS_MATRIX_GET_SIZE_F(M,LN,GN,IERR)
! CALL LIS_MATRIX_SET_ALL_F(0.0D0,M,IERR) 
! CALL LIS_MATRIX_SET_TYPE_F(M,LIS_MATRIX_COO,IERR)
! CALL LIS_MATRIX_ASSEMBLE_F(M,IERR)
 
 CALL LIS_ESOLVER_GET_EVECTORS_F(ESOLVER,M,IERR)
 CALL LIS_OUTPUT_MATRIX_F(M,LIS_FMT_MM,'D:\EIGENVECTORS.TXT',IERR);

 !## get eigenvalues
 CALL LIS_ESOLVER_CREATE_F(V,IERR)
 CALL LIS_VECTOR_DUPLICATE_F(A,V,IERR) 
 CALL LIS_ESOLVER_GET_EVALUES_F(ESOLVER,V,IERR)
 CALL LIS_VECTOR_GET_VALUES_F(V,1,N,EIGVALUES,IERR)
 
! CALL LIS_VECTOR_PRINT_F(EVALUE,IERR)

 CALL LIS_OUTPUT_VECTOR_F(V,LIS_FMT_MM,'D:\EIGENVALUES.TXT',IERR);
 CALL LIS_ESOLVER_OUTPUT_RHISTORY_F(ESOLVER,'D:\EIGENSOLVER.TXT',IERR);
 
 TOT=SUM(EIGVALUES); DO I=1,N
  EIGVALUES(I)=100.0D0*EIGVALUES(I)/TOT
  IF(I.GT.1)EIGVALUES(I)=EIGVALUES(I)+EIGVALUES(I-1)
 ENDDO

 DO I=1,N
  IF(EIGVALUES(I).GT.ESTOP)EXIT
 ENDDO
 ALLOCATE(EIGV(I-1,N),EIGW(I-1),TMP(N))
 
 !## get each eigenvector up to %
 CALL LIS_ESOLVER_CREATE_F(Y,IERR)
 CALL LIS_VECTOR_DUPLICATE_F(A,Y,IERR) 
 DO I=1,N
  IF(EIGVALUES(I).GT.ESTOP)EXIT
  EIGW(I)=EIGVALUES(I)
  CALL LIS_VECTOR_SET_ALL_F(0.0D0,X,IERR) 
!  CALL LIS_VECTOR_SET_ALL_F(0.0D0,Y,IERR) 
  CALL LIS_VECTOR_SET_VALUE_F(LIS_INS_VALUE,I,1.0D0,X,IERR) 
  CALL LIS_MATVEC_F(M,X,Y,IERR)
  CALL LIS_VECTOR_GET_VALUES_F(Y,1,N,TMP,IERR)
!  CALL LIS_VECTOR_GET_VALUES_F(Y,1,N,EIGV(:,I),IERR)
  DO J=1,N
   EIGV(I,J)=TMP(J)
  ENDDO
  CALL LIS_OUTPUT_VECTOR_F(Y,LIS_FMT_MM,'D:\EIGENVECTORS'//TRIM(VTOS(I))//'.TXT',IERR)
 ENDDO
! EIGV=TRANSPOSE(EIGV)
 
 CALL LIS_ESOLVER_DESTROY_F(ESOLVER,IERR) 
 
 DEALLOCATE(EIGVALUES,EIGERROR,TMP)
 CALL LIS_MATRIX_DESTROY_F(A,IERR)
 CALL LIS_MATRIX_DESTROY_F(M,IERR)
 CALL LIS_VECTOR_DESTROY_F(X,IERR)
 CALL LIS_VECTOR_DESTROY_F(Y,IERR)
 CALL LIS_VECTOR_DESTROY_F(V,IERR)

 END SUBROUTINE LIS_EIGENVALUE
 
 END MODULE MOD_LIS