program test

use f90_kind
use petsc_stuff

! Include PETSc files

#include "finclude/petscvecdef.h"
#include "finclude/petscmatdef.h"
#include "finclude/petsckspdef.h"
use petscvec
use petscmat
use petscksp
implicit none

Mat :: D_petsc  ! Note D not square
Mat :: M_petsc
Mat :: P_petsc  ! Our matrix of interest (P = D * M * D^T)
Vec :: x_petsc  ! Solution
Vec :: b_petsc  ! Rhs
KSP :: ksp_M
PetscErrorCode :: ierr

integer :: no_rows,iter_done,row,col,max_iter,final_iter,option
real(dp) :: tol,mat_entry

call init_petsc

max_iter = 10
no_rows = 4
tol = 1.d-5
option = 1

! Initialize the Matrices and vectors

call init_petsc_matrix(M_petsc,6*no_rows,6*no_rows,6*no_rows)
call init_petsc_matrix(D_petsc,no_rows,6*no_rows,6*no_rows)
call init_petsc_matrix(P_petsc,no_rows,no_rows,no_rows)
call init_petsc_vector(x_petsc,no_rows)
call init_petsc_vector(b_petsc,no_rows)

! Fill and view matrix D (Removing the zeros leads to solvable problem)

call addcoef_petsc_matrix(D_petsc,1,          1,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,1,          2, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,          3,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,1,          4,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,1,          5,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,1,          6, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,          7, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,          8, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,          9, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,         10,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,1,         11, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,1,         12, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,         13,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,1,         14, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,         15, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,1,         16, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,         17, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,         18, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,         19, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,         20, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,         21, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,         22, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,         23, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,1,         24, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,          1, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,          2, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,          3, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,          4, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,2,          5, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,2,          6, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,          7,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,2,          8, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,          9,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,2,         10, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,2,         11,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,2,         12, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,         13, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,         14, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,         15, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,         16, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,         17, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,         18, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,         19,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,2,         20, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,         21, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,2,         22, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,         23, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,2,         24, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,          1, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,3,          2, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,          3, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,3,          4, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,          5, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,          6, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,          7, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,          8, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,          9, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,         10, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,         11, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,         12, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,         13, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,3,         14, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,         15,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,3,         16,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,3,         17,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,3,         18, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,         19, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,         20, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,         21, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,3,         22,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,3,         23, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,3,         24, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,          1, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,          2, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,          3, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,          4, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,          5, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,          6, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,          7, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,4,          8, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,          9, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,4,         10, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,         11, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,         12, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,         13, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,         14, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,         15, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,         16, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,4,         17, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,4,         18, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,         19, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,4,         20, 0.000000000000000d+000)
call addcoef_petsc_matrix(D_petsc,4,         21,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,4,         22, 0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,4,         23,-0.250000000000000d0)
call addcoef_petsc_matrix(D_petsc,4,         24, 0.000000000000000d+000)
call finalize_petsc_matrix(D_petsc)

! Fill and view matrix M

call addcoef_petsc_matrix(M_petsc,1,1,4.d0)
call addcoef_petsc_matrix(M_petsc,2,2,12.d0)
call addcoef_petsc_matrix(M_petsc,3,3,12.d0)
call addcoef_petsc_matrix(M_petsc,4,4,4.d0)
call addcoef_petsc_matrix(M_petsc,5,5,12.d0)
call addcoef_petsc_matrix(M_petsc,6,6,12.d0)
call addcoef_petsc_matrix(M_petsc,7,7,4.d0)
call addcoef_petsc_matrix(M_petsc,8,8,12.d0)
call addcoef_petsc_matrix(M_petsc,9,9,12.d0)
call addcoef_petsc_matrix(M_petsc,10,10,4.d0)
call addcoef_petsc_matrix(M_petsc,11,11,12.d0)
call addcoef_petsc_matrix(M_petsc,12,12,12.d0)
call addcoef_petsc_matrix(M_petsc,13,13,4.d0)
call addcoef_petsc_matrix(M_petsc,14,14,12.d0)
call addcoef_petsc_matrix(M_petsc,15,15,12.d0)
call addcoef_petsc_matrix(M_petsc,16,16,4.d0)
call addcoef_petsc_matrix(M_petsc,17,17,12.d0)
call addcoef_petsc_matrix(M_petsc,18,18,12.d0)
call addcoef_petsc_matrix(M_petsc,19,19,4.d0)
call addcoef_petsc_matrix(M_petsc,20,20,12.d0)
call addcoef_petsc_matrix(M_petsc,21,21,12.d0)
call addcoef_petsc_matrix(M_petsc,22,22,4.d0)
call addcoef_petsc_matrix(M_petsc,23,23,12.d0)
call addcoef_petsc_matrix(M_petsc,24,24,12.d0)
call finalize_petsc_matrix(M_petsc)

if (option == 1) then

! Construct final matrix = 1.d-2 * D * M^-1 * D^T

  call triple_mat_prod(D_petsc,M_petsc,P_petsc,1.d-2)

else

! Fill by hand

  do row=1,4
    call addcoef_petsc_matrix(P_petsc,row,row,4.d-2)
  enddo
  call addcoef_petsc_matrix(P_petsc,1,2,-2.d-2)
  call addcoef_petsc_matrix(P_petsc,1,3,-2.d-2)
  call addcoef_petsc_matrix(P_petsc,2,1,-2.d-2)
  call addcoef_petsc_matrix(P_petsc,2,4,-2.d-2)
  call addcoef_petsc_matrix(P_petsc,3,1,-2.d-2)
  call addcoef_petsc_matrix(P_petsc,3,4,-2.d-2)
  call addcoef_petsc_matrix(P_petsc,4,2,-2.d-2)
  call addcoef_petsc_matrix(P_petsc,4,3,-2.d-2)
  call finalize_petsc_matrix(P_petsc)
endif

call MatView(P_petsc,PETSC_VIEWER_STDOUT_SELF,ierr)

! Print in more detail by using MatGetValues

do row=1,4
do col=1,4
  call MatGetValues(P_petsc,1,row-1,1,col-1,mat_entry,ierr)
  print *,row,col,mat_entry
enddo
enddo

! Fill and view rhs vector

call addcoef_petsc_vector(b_petsc,1,+1.d0)
call addcoef_petsc_vector(b_petsc,2,-1.d0)
call addcoef_petsc_vector(b_petsc,3, 0.d0)
call addcoef_petsc_vector(b_petsc,4, 0.d0)
call finalize_petsc_vector(b_petsc)
call VecView(b_petsc,PETSC_VIEWER_STDOUT_SELF,ierr)

! Finalize the solver and set options

call finalize_petsc_solver(P_petsc,ksp_M)

! Solve

call solve_petsc_system(ksp_M,b_petsc,x_petsc,tol,max_iter,final_iter)

if (final_iter >= 0) then
  print *,'convergence in',final_iter
else
  print *,'Maximum no iterations reached'
endif

! Final result

call VecView(x_petsc,PETSC_VIEWER_STDOUT_SELF,ierr)

call stop_petsc

end program test