Commit c1c82a0a authored by mozul's avatar mozul

[SPARSE] add tests of sparse linear algebra libraries

parent 8b9fb763
cmake_minimum_required(VERSION 2.8.0)
project(Systems Fortran)
find_package(LAPACK REQUIRED)
# looking for mumps
message(STATUS "looking for mumps in ${MP}")
find_library(MUMPS_FOUND dmumps HINTS ${MP} PATH_SUFFIXES lib NO_DEFAULT_PATH REQUIRED)
find_library(MUMPS_COMMON_FOUND mumps_common HINTS ${MP} PATH_SUFFIXES lib NO_DEFAULT_PATH REQUIRED)
find_library(PORD_FOUND pord HINTS ${MP} PATH_SUFFIXES lib NO_DEFAULT_PATH REQUIRED)
find_library(MPISEQ_FOUND mpiseq HINTS ${MP} PATH_SUFFIXES lib libseq REQUIRED)
find_path(MUMPS_INCLUDE_DIRS dmumps_struc.h NO_DEFAULT_PATH HINTS ${MP} ${MD} PATH_SUFFIXES include)
set(MUMPS_LIBRARIES ${MUMPS_FOUND} ${MUMPS_COMMON_FOUND} ${PORD_FOUND} ${MPISEQ_FOUND})
# looking for petsc
# mumps example
if( ${MP} )
include_directories(${MUMPS_INCLUDE_DIRS})
add_library( mumpsAlgebra dataStructure.f90 mumps_binding.f90)
target_link_libraries( mumpsAlgebra ${MUMPS_LIBRARIES} ${LAPACK_LIBRARIES})
add_executable( mumpsTest reference_test.f90 )
target_link_libraries( mumpsTest mumpsAlgebra )
endif( ${MP} )
# petsc example
if( ${PP} )
include_directories(${PETSC_INCLUDE_DIRS})
add_library( petscAlgebra dataStructure.f90 petsc_binding.f90)
target_link_libraries( petscAlgebra ${PETSC_LIBRARIES} ${LAPACK_LIBRARIES})
add_executable( petscTest reference_test.f90 )
target_link_libraries( petscTest petscAlgebra )
endif( ${PP} )
Le but de ce petit projet est de tester des bibliothèques
de résolutions de systèmes linéaires creux dans
l'utilisation qu'on en veut.
C'est à dire, avoir une collection de système où la
résolution des systèmes se fait dans une boucle openMP.
Descriptions des fichiers:
- reference_test.f90 : programmes principales avec les boucles
openMP pour la résolution
- dataStructure.f90 : défintion d'une structure de données de
système linéaire et initialisation des systmèmes
- xxxx_binding.f90 : binding sur une bibliothèque de résolution
de système creux, quasiment identique aux modules de binding de LMGC90
module dataStructure
use LinearAlgebra, only : sparse_mat
implicit none
type, public :: dstruc
logical :: sym
integer(kind=4) :: nb_dofs, nb_nz
integer(kind=4), dimension(:), pointer :: i_ind => null()
integer(kind=4), dimension(:), pointer :: j_ind => null()
real(kind=8) , dimension(:), pointer :: val => null()
real(kind=8) , dimension(:), pointer :: rhs => null()
type(sparse_mat) :: sparse_sys
end type
contains
subroutine systems_generation(systems)
implicit none
type(dstruc), dimension(:) :: systems
!
integer(kind=4) :: i_sys
do i_sys = 1, size(systems)
systems(i_sys)%nb_dofs = 5
systems(i_sys)%nb_nz = 12
allocate(systems(i_sys)%i_ind(systems(i_sys)%nb_nz))
allocate(systems(i_sys)%j_ind(systems(i_sys)%nb_nz))
allocate(systems(i_sys)%val(systems(i_sys)%nb_nz))
allocate(systems(i_sys)%rhs(systems(i_sys)%nb_dofs))
systems(i_sys)%i_ind = (/1, 2, 4, 5, 2, 1, 5, 3, 2, 3, 1, 3/)
systems(i_sys)%j_ind = (/2, 3, 3, 5, 1, 1, 2, 4, 5, 2, 3, 3/)
systems(i_sys)%val = (/3.0, -3.0, 2.0, 1.0, 3.0, 2.0, 4.0, 2.0, 6.0, -1.0, 4.0, 1.0/)
systems(i_sys)%rhs = (/20.0, 24.0, 9.0, 6.0,13.0 /)
end do
end subroutine
end module
module LinearAlgebra
implicit none
private
include 'dmumps_struc.h'
type, public :: sparse_mat
private
type(dmumps_struc) :: mumps_par
end type sparse_mat
public declare, build, &
solve, erase
contains
subroutine declare(matrix, nb_dofs, nb_non_zero, i_indices, j_indices, is_sym, info)
implicit none
type(sparse_mat) :: matrix
integer(kind=4) :: nb_dofs, nb_non_zero, info
integer(kind=4), dimension(:), pointer :: i_indices
integer(kind=4), dimension(:), pointer :: j_indices
logical :: is_sym
! initializing a mumps instance
matrix%mumps_par%job = -1
matrix%mumps_par%comm = 0
! is input symetric ?
if (is_sym) then
matrix%mumps_par%sym = 1
else
matrix%mumps_par%sym = 0
endif
! is parallelism used ?
matrix%mumps_par%par = 1
! let's do it
call DMUMPS( matrix%mumps_par )
! Analysis
matrix%mumps_par%job = 1
! setting verbosity to minimum
matrix%mumps_par%ICNTL(1)=0
matrix%mumps_par%ICNTL(2)=0
matrix%mumps_par%ICNTL(3)=0
matrix%mumps_par%ICNTL(4)=0
! assembled format input matrix
matrix%mumps_par%ICNTL(5 ) = 0
! triplet coordinates
matrix%mumps_par%ICNTL(18) = 0
matrix%mumps_par%N = nb_dofs
matrix%mumps_par%NZ = nb_non_zero
matrix%mumps_par%IRN => i_indices
matrix%mumps_par%JCN => j_indices
! scaling method
!matrix%mumps_par%ICNTL(8) = 8
! let's do it
call DMUMPS( matrix%mumps_par )
info = matrix%mumps_par%INFOG(1)
end subroutine
subroutine build(matrix, val, info)
implicit none
type(sparse_mat) :: matrix
real(kind=8) , dimension(:), pointer :: val
integer(kind=4) :: info
! Factorize
matrix%mumps_par%job = 2
matrix%mumps_par%A => val
! let's do it
call DMUMPS( matrix%mumps_par )
info = matrix%mumps_par%INFOG(1)
end subroutine
subroutine solve(matrix, rhs, info)
implicit none
type(sparse_mat) :: matrix
real(kind=8) , dimension(:), pointer :: rhs
integer(kind=4) :: info
! solve
matrix%mumps_par%job = 3
! rhs
matrix%mumps_par%RHS => rhs
! let's do it
call DMUMPS( matrix%mumps_par )
info = -matrix%mumps_par%INFOG(1)
end subroutine
subroutine erase(matrix)
implicit none
type(sparse_mat) :: matrix
matrix%mumps_par%job = -2
! let's do it
call DMUMPS( matrix%mumps_par )
end subroutine
end module
module LinearAlgebra
implicit none
private
#include <finclude/petsc.h90>
type, public :: sparse_mat
private
!> solution
Vec :: x
!> righ hand side
Vec :: b
!> linear system matrix
Mat :: A
! linear solver context
KSP :: ksp
! preconditioner context
PC :: pc
logical :: sym = .false.
integer(kind=4), dimension(:), pointer :: i_indices => null()
integer(kind=4), dimension(:), pointer :: j_indices => null()
end type sparse_mat
public sparse_init , &
sparse_declare, &
sparse_build , &
sparse_solve , &
sparse_erase , &
sparse_finalize
contains
subroutine sparse_init()
implicit none
!
integer(kind=4) :: ierr
character(len=32) :: IAM
! 12345678901234567890123456789012
IAM = 'SparseLinearAlgebra::sparse_init'
call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
if( ierr > 0 ) then
print *, '[', IAM, '] PetscInitalize failed with value : ', ierr
stop
end if
end subroutine
subroutine declare(matrix,nb_dofs, nb_non_zero, i_indices, j_indices, is_sym, info)
implicit none
type(sparse_mat) :: matrix
integer(kind=4) :: nb_dofs, nb_non_zero, info
integer(kind=4), dimension(:), pointer :: i_indices
integer(kind=4), dimension(:), pointer :: j_indices
logical :: is_sym
!
PetscBool :: is_spd
real(kind=8) :: tol
character(len=35) :: IAM
! 12345678901234567890123456789012345
IAM = 'SparseLinearAlgebra::sparse_declare'
is_spd = PETSC_FALSE
tol = 1.e-5
! create right hand side
call VecCreate(PETSC_COMM_WORLD,matrix%b,info)
call VecSetSizes(matrix%b,PETSC_DECIDE,nb_dofs,info)
call VecSetFromOptions(matrix%b,info)
! create solution
call VecDuplicate(matrix%b,matrix%x,info)
! create left hand side
call MatCreate(PETSC_COMM_WORLD,matrix%A,info)
call MatSetSizes(matrix%A,PETSC_DECIDE,PETSC_DECIDE,nb_dofs,nb_dofs,info)
call MatSetFromOptions(matrix%A,info)
!if( is_sym ) then
! is_spd = PETSC_TRUE
! matrix%sym = .true.
!end if
!call MatSetOption(matrix%A,MAT_SPD,is_spd)
call MatSetUp(matrix%A,info)
! create linear solver
call KSPCreate(PETSC_COMM_WORLD,matrix%ksp,info)
call KSPGetPC(matrix%ksp,matrix%pc,info)
call PCSetType(matrix%pc,PCJACOBI,info)
!call KSPSetTolerances(matrix%ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,info)
matrix%i_indices => i_indices
matrix%j_indices => j_indices
end subroutine
subroutine build(matrix, val, info)
implicit none
type(sparse_mat) :: matrix
real(kind=8) , dimension(:), pointer :: val
integer(kind=4) :: info
!
integer(kind=4) :: i
character(len=33) :: IAM
! 123456789012345678901234567890123
IAM = 'SparseLinearAlgebra::sparse_build'
do i = 1, size(val)
call MatSetValues(matrix%A,1,matrix%i_indices(i)-1,1,matrix%j_indices(i)-1, &
val(i),INSERT_VALUES, info)
end do
if( .not. matrix%sym ) then
do i = 1, size(val)
call MatSetValues(matrix%A,1,matrix%j_indices(i)-1,1,matrix%i_indices(i)-1, &
val(i),INSERT_VALUES, info)
end do
end if
call MatAssemblyBegin(matrix%A,MAT_FINAL_ASSEMBLY,info)
call MatAssemblyEnd(matrix%A,MAT_FINAL_ASSEMBLY,info)
call KSPSetOperators(matrix%ksp,matrix%A,matrix%A,info)
call KSPSetFromOptions(matrix%ksp,info)
end subroutine
subroutine solve(matrix, rhs, info)
implicit none
type(sparse_mat) :: matrix
real(kind=8) , dimension(:), pointer :: rhs
integer(kind=4) :: info
!
integer(kind=4) :: i
character(len=33) :: IAM
! 123456789012345678901234567890123
IAM = 'SparseLinearAlgebra::sparse_solve'
do i = 1, size(rhs)
call VecSetValues(matrix%b,1,i-1,rhs(i),INSERT_VALUES,info)
end do
call VecAssemblyBegin(matrix%b,info)
call VecAssemblyEnd(matrix%b,info)
call KSPSolve(matrix%ksp,matrix%b,matrix%x,info)
!activate only if verbose...
!call KSPView(matrix%ksp,PETSC_VIEWER_STDOUT_WORLD,info)
do i = 1, size(rhs)
call VecGetValues(matrix%x,1,i-1,rhs(i),info)
end do
end subroutine
subroutine erase(matrix)
implicit none
type(sparse_mat) :: matrix
!
integer(kind=4) :: info
character(len=33) :: IAM
! 123456789012345678901234567890123
IAM = 'SparseLinearAlgebra::sparse_erase'
call VecDestroy(matrix%x,info)
call VecDestroy(matrix%b,info)
call MatDestroy(matrix%A,info)
call KSPDestroy(matrix%ksp,info)
end subroutine
subroutine sparse_finalize()
implicit none
!
integer(kind=4) :: ierr
character(len=36) :: IAM
! 123456789012345678901234567890123456
IAM = 'SparseLinearAlgebra::sparse_finalize'
call PetsCFinalize(ierr)
if( ierr > 0 ) then
print *, '[', IAM, '] PetscFinalize failed with value : ', ierr
stop
end if
end subroutine
end module
program test_systems
use dataStructure, only : dstruc, systems_generation
use LinearAlgebra, only : declare, build, solve, erase
implicit none
! number of systems
integer(kind=4), parameter :: nb_systems = 16
! list of sparse systems
type(dstruc), dimension(nb_systems) :: systems
!
integer(kind=4) :: i_sys, info
call systems_generation(systems)
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i_sys)
!$OMP DO SCHEDULE(RUNTIME)
do i_sys = 1, nb_systems
call declare( systems(i_sys)%sparse_sys, systems(i_sys)%nb_dofs, &
systems(i_sys)%nb_nz, systems(i_sys)%i_ind , &
systems(i_sys)%j_ind, systems(i_sys)%sym, info )
if( info /= 0 ) then
print *, 'Error when declaring system : ', i_sys
stop 1
end if
end do
!$OMP END DO
!$OMP END PARALLEL
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i_sys)
!$OMP DO SCHEDULE(RUNTIME)
do i_sys = 1, nb_systems
call build( systems(i_sys)%sparse_sys, systems(i_sys)%val, info)
if( info /= 0 ) then
print *, 'Error when building system : ', i_sys
stop 2
end if
end do
!$OMP END DO
!$OMP END PARALLEL
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i_sys)
!$OMP DO SCHEDULE(RUNTIME)
do i_sys = 1, nb_systems
call solve( systems(i_sys)%sparse_sys, systems(i_sys)%rhs, info)
print *, i_sys, systems(i_sys)%rhs
if( info /= 0 ) then
print *, 'Error when solving system : ', i_sys
stop 3
end if
end do
!$OMP END DO
!$OMP END PARALLEL
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i_sys)
!$OMP DO SCHEDULE(RUNTIME)
do i_sys = 1, nb_systems
deallocate(systems(i_sys)%i_ind)
deallocate(systems(i_sys)%j_ind)
deallocate(systems(i_sys)%val)
deallocate(systems(i_sys)%rhs)
call erase( systems(i_sys)%sparse_sys)
end do
!$OMP END DO
!$OMP END PARALLEL
end program
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment