program build_nullbasis_petsc_mumps
implicit none
#include <petsc/finclude/petscsys.h>
#include "petsc/finclude/petscvec.h"
#include "petsc/finclude/petscmat.h"
#include "petsc/finclude/petscpc.h"
#include "petsc/finclude/petscksp.h"
!---------------------------------------------!
!Variables declarations
!---------------------------------------------!
   Mat :: mat_c, mat_c_temp, F, x, b, mat_temp, mat_temp_2, mat_temp_3
   PC :: pc
   KSP :: ksp
   PetscInt :: n, m, dim_mat, nbcol, ival, icntl, infog28, II, Istart, Iend, nb_procs, i
   PetscInt, dimension(:), allocatable :: colonnes
   PetscScalar, dimension(:), allocatable :: valeurs
   PetscReal, dimension(:), allocatable :: norms_c, norms_ct
   PetscReal :: tol, eps
   PetscBool :: base_ok
   PetscErrorCode :: ierr
   PetscMPIInt :: rank
   PetscViewer :: viewer
!---------------------------------------------!
!MPI initialization
!---------------------------------------------!
   call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
   call MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr)
   call MPI_Comm_size(PETSC_COMM_WORLD, nb_procs, ierr)
!---------------------------------------------!
!Matrices creation
!---------------------------------------------!
   do i = 1,4
   write(*,*) 'BEGIN PROC', rank
   if (rank == 0) then
      write(*,*) 'ITERATION', i
   end if
!----------Read file (mat_c)----------!
   call MatCreate(PETSC_COMM_WORLD, mat_c, ierr)
   call MatSetType(mat_c, MATAIJ, ierr)
   call PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mat_c_bin.txt",0,viewer, ierr) 
   call MatLoad(mat_c, viewer,ierr)
   call PetscViewerDestroy(viewer, ierr)
!   call MatView(mat_c,PETSC_VIEWER_STDOUT_WORLD,ierr) 

   call MatGetSize(mat_c, m, n, ierr)
   write(*,*) 'm,n', m,n

   dim_mat = max(m,n)
   allocate(colonnes(dim_mat))
   allocate(valeurs(dim_mat))
!----------Create mat_c_temp (square matrix) for MUMPS----------!
   call MatCreate(PETSC_COMM_WORLD, mat_c_temp, ierr)
   call MatSetSizes(mat_c_temp, PETSC_DECIDE, PETSC_DECIDE, dim_mat, dim_mat, ierr)
   call MatSetType(mat_c_temp, MATAIJ, ierr)
!   call MatMPIAIJSetPreallocation(mat_c_temp, dim_mat, PETSC_NULL_INTEGER, dim_mat, PETSC_NULL_INTEGER, ierr)
   call MatSetUp(mat_c_temp, ierr)
!   call MatZeroEntries(mat_c_temp, ierr)
   call MatGetOwnershipRange(mat_c, Istart, Iend, ierr)
   do II = Istart, Iend - 1
         colonnes(:) = 0
         valeurs(:) = 0
         nbcol = 0
      call MatGetRow(mat_c, II, nbcol, colonnes, valeurs, ierr)
      call MatSetValues(mat_c_temp, 1, II, nbcol, colonnes(1:nbcol:1), valeurs(1:nbcol:1), INSERT_VALUES, ierr)
      call MatRestoreRow(mat_c, II, nbcol, colonnes, valeurs, ierr)
   end do
   call MatAssemblyBegin(mat_c_temp, MAT_FINAL_ASSEMBLY, ierr)
   call MatAssemblyEnd(mat_c_temp, MAT_FINAL_ASSEMBLY, ierr)
!---------------------------------------------!
!KSP
!---------------------------------------------!
   call KSPCreate(PETSC_COMM_WORLD, ksp, ierr)
   call KSPSetOperators(ksp, mat_c_temp, mat_c_temp, ierr)
   tol = 1.e-7
   call KSPSetTolerances(ksp, tol, PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL, PETSC_DEFAULT_INTEGER, ierr)
   call KSPSetType(ksp, KSPPREONLY, ierr)
   call KSPGetPC(ksp, pc, ierr)
   call PCSetType(pc, PCLU, ierr)
   call PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS, ierr)
   call PCFactorSetUpMatSolverPackage(pc, ierr)
   call PCFactorGetMatrix(pc, F, ierr)
!---------------------------------------------!
!MUMPS
!---------------------------------------------!
   icntl = 6
   ival = 0
   call MatMumpsSetIcntl(F, icntl, ival, ierr)
!----------------------------------------------
   icntl = 8
   ival = 77
   call MatMumpsSetIcntl(F, icntl, ival, ierr)
!----------------------------------------------
   icntl = 21
   if (nb_procs == 1) then
      ival = 0
   else
      ival = 1
   end if
   call MatMumpsSetIcntl(F, icntl, ival, ierr)
!----------------------------------------------
   icntl = 24
   ival = 1
   call MatMumpsSetIcntl(F, icntl, ival, ierr)
!----------------------------------------------
   icntl = 25
   ival = -1
   call MatMumpsSetIcntl(F, icntl, ival, ierr)
!----------------------------------------------
   if (rank == 0) then
      write(*,*) 'ECHO 1'
   end if
   call KSPSetUp(ksp, ierr)
   if (rank == 0) then
      write(*,*) 'ECHO 2'
   end if
!----------------------------------------------
   call MatMumpsGetInfog(F, 28, infog28, ierr)
   if (rank == 0) then
      write(*,*) 'INFOG(28):', infog28
   end if
!----------------------------------------------
   call MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, dim_mat, infog28, PETSC_NULL_SCALAR, x, ierr) !pour pouvoir utiliser MatMatSolve, il faut que X soit dense
   call MatSetType(x, MATDENSE, ierr)
   call MatSetUp(x, ierr)
   call MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY,ierr)
   call MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY,ierr)
   call MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, dim_mat, infog28, PETSC_NULL_SCALAR, b, ierr)
   call MatSetType(b, MATDENSE, ierr)
   call MatSetUp(b, ierr)
   call MatZeroEntries(b, ierr)
   call MatAssemblyBegin(b,MAT_FINAL_ASSEMBLY,ierr)
   call MatAssemblyEnd(b,MAT_FINAL_ASSEMBLY,ierr)
!---------------------------------------------!
!Solve
!---------------------------------------------!
   call MatMatSolve(F, b, x, ierr)
!---------------------------------------------!
!Write solution in .txt file
!---------------------------------------------!
   call PetscViewerASCIIOpen(PETSC_COMM_WORLD,"mat_x.txt",viewer, ierr)
   call MatView(x, viewer, ierr)
   call PetscViewerDestroy(viewer, ierr)
!---------------------------------------------!
! Verification
!---------------------------------------------!
   base_ok = PETSC_TRUE
   eps = 0.0001
   allocate(norms_c(n)) 
   allocate(norms_ct(n))
!----------------------------------------------
   call MatTranspose(mat_c_temp, MAT_INITIAL_MATRIX, mat_temp, ierr)
   call MatGetColumnNorms(mat_temp, NORM_2, norms_c, ierr)
!----------------------------------------------
   call MatConvert(x,MATAIJ, MAT_INPLACE_MATRIX, x, ierr)
   call MatMatMult(mat_c_temp, x, MAT_INITIAL_MATRIX, PETSC_DEFAULT_REAL, mat_temp_2, ierr)
   call MatTranspose(mat_temp_2, MAT_INITIAL_MATRIX, mat_temp_3, ierr)
   call MatGetColumnNorms(mat_temp_3, NORM_2, norms_ct, ierr)
!----------------------------------------------
   if (rank == 0) then
      do II = 1, m
         if (norms_c(II) < eps ) then 
            write(*, '(A11,I3,A3,E11.4)') 'LINE ', II, ' : ', norms_ct(II)/norms_c(II)
            base_ok = PETSC_FALSE
         end if
      end do
!----------------------------------------------
      if (base_ok) then
         write(*,*) 'BASIS OK', rank
      else
         write(*,*) 'BASIS NO OK', rank
      end if
   end if

!---------------------------------------------!
!End
!--------------------------------------------!
   call MatDestroy(mat_temp, ierr)
   call MatDestroy(mat_temp_2, ierr)
   call MatDestroy(mat_temp_3, ierr)
   call MatDestroy(x, ierr)
   call MatDestroy(b, ierr)
   call MatDestroy(mat_c, ierr)
   call MatDestroy(mat_c_temp, ierr)
   call KSPDestroy(ksp, ierr)
   deallocate(colonnes)
   deallocate(valeurs)
   deallocate(norms_c)
   deallocate(norms_ct)
   write(*,*) 'END PROC  ',rank
   end do


   call PetscFinalize(ierr)
end program build_nullbasis_petsc_mumps