!>  
!!
!!
!!<B>      SOLVES A LINEAR PROBLEM</B>, using \ref linsysslv_mod::linsys_solve "linSys_solve"  
!!
!!  This is a tutorial example for \ref  linsysslv_mod "linSysSlv_mod".
!!
!!  Shows various methods to solve \f$A x = b\f$
!!
!> @author Charles PIERRE, September 2021
!>  

program linSysSlv_mod_expl

  use cumin  

  implicit none

  !! !!!!!!!!!!!!!!!!! VARAIBLE DEFINITION : BEGIN
  !!
  !!   N         = problem size
  !!
  integer, parameter :: N = 100

  
  !! MATRICES
  !!
  !!   A     = system matrix, csr format
  !!   b, x  = rhs and solution
  !!
  type(csr)      :: A 
  real(RP), dimension(N) :: b, x

  !! SOLVER
  !!
  !!   slv   = linear system solver definition
  !!
  type(linSysSlv) :: lss

  !!
  !!
  !! !!!!!!!!!!!!!!!!! VARAIBLE DEFINITION : END

  write(*,*)""
  call paragraph("TUTORIAL EXAMPLE ON: linSysSlv_mod")
  write(*,*)""
  write(*,*)"  RESOLUTION OF THE LINEAR SYSTEM Ax = b"

  !! !!!!!!!!!!!!!!!!!! SYSTEM MATRIX
  !!
  call paragraph("SYSTEM MATRIX, RHS")
  !!
  A = triDiag_SDP()
  b = 1.0_RP


  !! !!!!!!!!!!!!!!!!!! SOLVE, CG, NO PRECONDITIONNING
  !!
  call paragraph("CG (Conjugate Gradient)") 
  !!
  !! solver definition
  call linSysSlv_create(lss, LINSYSSLV_CG, verb=2)
  call set(lss, tol = 1E-8_RP, itMax = 500, verb=2)
  call print(lss)
  !!
  !! linear system inversion
  x = 0.0_RP  !! initial guess
  call linSys_solve(x, A, b, lss, verb=2)


  !! !!!!!!!!!!!!!!!!!! CG,  iC(0) PRECONDITIONNING
  !!
  call paragraph("CG, iC(0) PRECONDITIONNING")
  !!
  !! solver definition
  call linSysSlv_create(lss, LINSYSSLV_CG, verb=2)
  call set(lss, tol = 1E-8_RP, itMax = 500, verb=2)
  call set_precond(lss, type=precond_ICK, verb=2)
  call print(lss)
  !!
  !! linear system inversion
  x = 0.0_RP  !! initial guess
  call linSys_solve(x, A, b, lss, verb=2)

#ifdef WMUMPS
  !! !!!!!!!!!!!!!!!!!! MUMPS LLT
  !!
  call paragraph("MUMPS LLT")
  !!
  !! solver definition
  call linSysSlv_create(lss, MUMPS_LLT, verb=2)
  call print(lss)
  !!
  !! linear system inversion
  call linSys_solve(x, A, b, lss, verb=2)
#endif


  print*
  print*
  print*

contains 

  !> m = finite difference matrix for -Delta u + u = f
  !> with Dirichlet condition.
  !>
  function triDiag_SDP() result(m)
    type(csr) :: m

    real(RP), dimension(N,N) :: A2

    integer :: ii

    !! A2 u = -Delta u
    !!
    A2 = 0.0_RP
    !!
    do ii=2, N-1       
       A2(ii, ii-1) = -1.0_RP 
       A2(ii, ii  ) =  2.0_RP 
       A2(ii, ii+1) = -1.0_RP 
    end do
    !!
    A2(1,1) =  2.0_RP 
    A2(1,2) = -1.0_RP 
    !!
    A2(N, N  ) =  2.0_RP 
    A2(N, N-1) = -1.0_RP 

    call matrix_to_csr(m, A2)

  end function triDiag_SDP


end program linSysSlv_mod_expl