Hi,I'm a new user of PETSc and I try to use it with MUMPS functionalities to compute a nullbasis. I wrote a code where I compute 4 times the same nullbasis. It does work well when I run it with several procs but with only one processor I get an error on the 2nd iteration when KSPSetUp is called. Furthermore when it is run with a debugger ( --with-debugging=yes), it works fine with one or several processors. Have you got any idea about why it doesn't work with one processor and no debugger?
Thanks. Constantin. PS: You can find the code and the files required to run it enclosed.
mat_c_bin.txt
Description: Binary data
-matload_block_size 1
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)
call PetscViewerDestroy(viewer, ierr)
call MatGetSize(mat_c, m, n, ierr)
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, MATMPIAIJ, 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, MATMPIDENSE, 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, MATMPIDENSE, 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, MATMPIAIJ, MAT_REUSE_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
