!>
!> @file tsparse2.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs 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 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs 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 Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
!   Simple 2D Poisson using 5 points FD
!
  USE pardiso_bsplines
  IMPLICIT NONE
!
  TYPE(pardiso_mat) :: amat, bmat
  DOUBLE PRECISION, ALLOCATABLE :: arow(:), rhs(:), sol(:)
  DOUBLE PRECISION, ALLOCATABLE :: arown(:,:), rhsn(:,:), soln(:,:)
  INTEGER :: nx=5, ny=4, n
  INTEGER :: i, j, irow, jcol
!
  DOUBLE PRECISION   :: mem, seconds
!
  n = nx*ny
  CALL init(n, 1, amat)                ! Non-symmetric matrix
  CALL init(n, 1, bmat, nlsym=.TRUE.)  ! Symmetric matrix
  ALLOCATE(arow(n), arown(n,2))
  ALLOCATE(rhs(n), rhsn(n,2))
  ALLOCATE(sol(n), soln(n,2))
!
!  Construct the FD matrix amat, using sparse rows (linked lists)
!
  DO j=1,ny
     DO i=1,nx
        arow = 0.0d0
        irow = numb(i,j)
        arow(irow) = 4.0d0
        IF(i.GT.1)  arow(numb(i-1,j)) = -1.0d0
        IF(i.LT.nx) arow(numb(i+1,j)) = -1.0d0
        IF(j.GT.1)  arow(numb(i,j-1)) = -1.0d0
        IF(j.LT.ny) arow(numb(i,j+1)) = -1.0d0
        rhs(irow) = SUM(arow)
        CALL putrow(amat, irow, arow) ! General matrix
        CALL putrow(bmat, irow, arow) ! Symmetric matrix
     END DO
  END DO
!
!   Print the matrices
!
  WRITE(*,'(/a)') 'Matrix A'
  DO i=1,n
     CALL getrow(amat, i, arow)
     WRITE(*,'(30f4.0)') arow
  END DO
  PRINT*, 'nnz from get_count', get_count(amat)
!
  WRITE(*,'(/a)') 'Matrix B'
  DO i=1,n
     CALL getrow(bmat, i, arow)
     WRITE(*,'(30f4.0)') arow
  END DO
  PRINT*, 'nnz from get_count', get_count(bmat)
!
!   Factor the matrix using Pardiso
!
  CALL factor(amat, nlmetis=.TRUE.)
  WRITE(*,'(/a,i5)') 'Number of nonzeros in factors of A  = ',amat%p%iparm(18)
  WRITE(*,'(a,i5)')  'Number of factorization MFLOPS      = ',amat%p%iparm(19)
!
  CALL factor(bmat, nlmetis=.TRUE.)
  WRITE(*,'(/a,i5)') 'Number of nonzeros in factors of B  = ',bmat%p%iparm(18)
  WRITE(*,'(a,i5)')  'Number of factorization MFLOPS      = ',bmat%p%iparm(19)
!
  WRITE(*,'(/a/(10f8.4))') 'rhs', rhs
!
!   Backsolve Ax = b, using Pardiso
!
  sol = rhs
  CALL bsolve(amat, sol, debug=.FALSE.)
  WRITE(*,'(/a/(10f8.4))') 'sol (non-sym)', sol
  WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(sol-1.0d0))
!
!   Backsolve Bx = b, using Pardiso
!
  sol = rhs
  CALL bsolve(bmat, sol, debug=.FALSE.)
  WRITE(*,'(/a/(10f8.4))') 'sol (sym)', sol
  WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(sol-1.0d0))
!
  arow = vmx(amat, sol)
  WRITE(*,'(/a/(10f8.4))') 'A*x (non-sym)', arow
  WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(arow-rhs))
!
  arow = vmx(bmat, sol)
  WRITE(*,'(/a/(10f8.4))') 'B*x (sym)', arow
  WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(arow-rhs))
!
!   Multiple RHS
!
  rhsn(:,1) = -rhs(:)
  rhsn(:,2) = 2.d0*rhs(:)
  CALL bsolve(amat, rhsn, soln)
  WRITE(*,'(/a/(10f8.4))') 'soln (non-sym)', soln
  WRITE(*,'(a,2(1pe12.3))') 'Error', MAXVAL(ABS(soln(:,1)+1.0d0)), &
       &                  MAXVAL(ABS(soln(:,2)-2.0d0))
!
  CALL bsolve(bmat, rhsn, soln)
  WRITE(*,'(/a/(10f8.4))') 'soln (sym)', soln
  WRITE(*,'(a,2(1pe12.3))') 'Error', MAXVAL(ABS(soln(:,1)+1.0d0)), &
       &                  MAXVAL(ABS(soln(:,2)-2.0d0))
!
  arown = vmx(amat, soln)
  WRITE(*,'(/a/(10f8.4))') 'A*x (non-sym)', arown
  WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(arown-rhsn))
!
  arown = vmx(bmat, soln)
  WRITE(*,'(/a/(10f8.4))') 'A*x (sym)', arown
  WRITE(*,'(a,1pe12.3)') 'Error', MAXVAL(ABS(arown-rhsn))
!
!   Clean up
!
  DEALLOCATE(arow,arown)
  DEALLOCATE(rhs,rhsn)
  DEALLOCATE(sol,soln)
  CALL destroy(amat)
  CALL destroy(bmat)
CONTAINS
  INTEGER FUNCTION numb(i,j)
    INTEGER, INTENT(in) :: i, j
    numb = (j-1)*nx + i
  END FUNCTION numb
!
END PROGRAM main