!>
!> @file zpardiso_ex1.f
!>
!> @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
      USE wsmp_bsplines
      USE pardiso_bsplines
      IMPLICIT NONE
c
      INTEGER :: n=9
      INTEGER ia(10)
      INTEGER ja(29)
      COMPLEX*16 avals(29)
      COMPLEX*16 b(9), sol(9), arow(9)
c
c$$$      type(zwsmp_mat) :: mat
      type(zpardiso_mat) :: mat
      integer :: i, k
c
      DATA ia /1,5,9,13,17,21,25,27,29,30/
      data ja 
     1        /1,          3,                      7,    8,
     2	             2,    3,                            8,    9,
     3                     3,                      7,    8,    9,
     4                           4,          6,    7,    8,
     5                                 5,    6,          8,    9,
     6                                       6,    7,    8,    9,
     7                                             7,    8,
     8                                                   8,    9,
     9                                                         9/
      data avals 
     1 /(14.d0,0.d0), (-1.d0,-5.d-2), (-1.d0,-5.0d-2), (-3.d0,-1.5d-1),
     2 (14.d0,0.d0), (-1.d0,-5.d-2), (-3.d0,-1.5d-1), (-1.d0,-5.0d-2),
     3 (16.d0,0.d0), (-2.d0,-1.d-1), (-4.d0,-2.0d-1), (-2.d0,-1.0d-1),
     4 (14.d0,0.d0), (-1.d0,-5.d-2), (-1.d0,-5.0d-2), (-3.d0,-1.5d-1),
     5 (14.d0,0.d0), (-1.d0,-5.d-2), (-3.d0,-1.5d-1), (-1.d0,-5.0d-2),
     6 (16.d0,0.d0), (-2.d0,-1.d-1), (-4.d0,-2.0d-1), (-2.d0,-1.0d-1),
     7 (16.d0,0.d0), (-4.d0,-2.d-1),
     8 (71.d0,0.d0), (-4.d0,-2.d-1),
     9 (16.d0,0.d0)/ 
c
c$$$      call init(n, 1, mat, nlherm=.false., nlsym=.true., nlpos=.true.)
      call init(n, 1, mat, nlherm=.true., nlpos=.true.)
      do i=1,n
         do k=ia(i),ia(i+1)-1
            call putele(mat, i, ja(k), avals(k))
         end do
      end do
c
      call factor(mat)
c
      print*, 'diff of val', cnorm2(avals-mat%val)
      print*, 'diff of ia', ia-mat%irow
      print*,' diff ja', ja-mat%cols
c
      print*, 'The RHS:'
      do i = 1, n
         call getrow(mat,i, arow)
         b(i) = sum(arow)
         print *,i,' : ',b(i)
      end do
      call bsolve(mat,b,sol)
      print *,'The solution of the system is as follows:'
      do i = 1, n
         print *,i,' : ',sol(i)
      end do
      print*, 'Residue =', cnorm2(vmx(mat,sol)-b)
      contains
      FUNCTION cnorm2(x)
      DOUBLE COMPLEX, INTENT(in) :: x(:)
      DOUBLE PRECISION :: cnorm2
      cnorm2 = SQRT(DOT_PRODUCT(x,x))    
      END FUNCTION cnorm2
c
      END PROGRAM main