[SPARSE] use of umfpack

systems_in_parallel/CMakeLists.txt

 
 cmake_minimum_required(VERSION 2.8.0)
 
-project(Systems Fortran)
+project(Systems Fortran C)
 
 find_package(LAPACK REQUIRED)
 
+set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fopenmp")
+
 # looking for mumps
 if( MP )
   message(STATUS "looking for mumps in ${MP}")
@@ -31,6 +33,15 @@ if( PP )
   message(STATUS "PETSc INCLUDE  : ${PETSC_C_INCLUDE_DIR} ${PETSC_F_INCLUDE_DIR}")
 endif( PP )
 
+# looking for umfpack
+if( UP )
+  message(STATUS "looking for umfpack")
+  find_library(UMFPACK_FOUND umfpack REQUIRED)
+  find_path(UMFPACK_INCLUDE_DIR umfpack.h)
+  message(STATUS "UMFPACK LIBRARIES: ${UMFPACK_FOUND}")
+  message(STATUS "UMFPACK INCLUDE  : ${UMFPACK_INCLUDE_DIR}")
+endif( UP )
+
 # mumps example
 if( MUMPS_FOUND )
   include_directories(${MUMPS_INCLUDE_DIR})
@@ -56,3 +67,25 @@ if( PETSC_FOUND )
   add_executable( petscTest reference_test.f90 )
   target_link_libraries( petscTest petscAlgebra )
 endif( PETSC_FOUND )
+
+# umfpack example
+if( UMFPACK_FOUND )
+
+  include_directories(${UMFPACK_INCLUDE_DIR})
+  include_directories(${CMAKE_CURRENT_SOURCE_DIR})
+
+  add_library(cbind umfpack_cbinding.c)
+  target_link_libraries( cbind ${UMFPACK_FOUND})
+
+  add_library( umfpackAlgebra dataStructure.f90 umfpack_binding.f90)
+  target_link_libraries( umfpackAlgebra cbind )
+
+  add_executable( umfpackTest reference_test.f90 )
+  target_link_libraries( umfpackTest umfpackAlgebra )
+
+  # simple C example of use
+  add_executable( umfpackCTest umfpack_test.c )
+  target_link_libraries( umfpackCTest ${UMFPACK_FOUND})
+
+endif( UMFPACK_FOUND )
+

systems_in_parallel/dataStructure.f90

@@ -6,7 +6,7 @@ module dataStructure
   implicit none
 
   type, public :: dstruc
-    logical :: sym
+    logical :: sym = .false.
     integer(kind=4) :: nb_dofs, nb_nz
     integer(kind=4), dimension(:), pointer :: i_ind => null()
     integer(kind=4), dimension(:), pointer :: j_ind => null()
@@ -36,9 +36,9 @@ contains
       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)%i_ind = (/  1,   1,   1,   2,    2,   2,    3,   3,   3,   4,   5,   5/)
+      systems(i_sys)%j_ind = (/  1,   2,   3,   1,    3,   5,    2,   3,   4,   3,   2,   5/)
+      systems(i_sys)%val   = (/2.0, 3.0, 4.0, 3.0, -3.0, 6.0, -1.0, 1.0, 2.0, 2.0, 4.0, 1.0/)
       systems(i_sys)%rhs   = (/20.0, 24.0, 9.0, 6.0,13.0 /)
 
     end do

systems_in_parallel/umfpack_binding.f90

+
+module LinearAlgebra
+
+  use iso_c_binding
+
+  implicit none 
+
+  private
+
+  interface
+    function umfpack_symbolic(nrow,ncol,ap,ai,ax,symbolic,control,info) bind(c, name="c_symbolic")
+      import c_int, c_double, c_ptr
+
+      integer(c_int), intent(in), value :: nrow
+      integer(c_int), intent(in), value :: ncol
+
+      type(c_ptr), value :: ap, ai, ax
+      type(c_ptr) :: symbolic, control
+
+      integer(c_int) :: umfpack_symbolic
+
+      real(c_double), dimension(90) :: info
+
+    end function
+  end interface
+
+  interface
+    function umfpack_numeric(ap,ai,ax,symbolic,numeric,control,info) bind(c, name="c_numeric")
+      import c_int, c_double, c_ptr
+
+      type(c_ptr), value :: ap, ai, ax
+      type(c_ptr) :: symbolic, numeric, control
+
+      integer(c_int) :: umfpack_numeric
+
+      real(c_double), dimension(90) :: info
+
+    end function
+  end interface
+
+  interface
+    function umfpack_solve(ap,ai,ax,x,b,numeric,control,info) bind(c, name="c_solve")
+      import c_int, c_double, c_ptr
+
+      type(c_ptr), value :: ap, ai, ax
+      type(c_ptr) :: x, b, numeric, control
+
+      integer(c_int) :: umfpack_solve
+
+      real(c_double), dimension(90) :: info
+
+    end function
+  end interface
+
+  interface
+    subroutine umfpack_free_symbolic(symbolic) bind(c, name="c_free_symbolic")
+      import c_ptr
+      type(c_ptr) :: symbolic
+    end subroutine
+  end interface
+
+  interface
+    subroutine umfpack_free_numeric(numeric) bind(c, name="c_free_numeric")
+      import c_ptr
+      type(c_ptr) :: numeric
+    end subroutine
+  end interface
+
+
+  type, public :: sparse_mat
+    private
+
+    integer(c_int) :: n
+
+    type(c_ptr) :: symbolic = c_null_ptr
+    type(c_ptr) :: numeric  = c_null_ptr
+
+    type(c_ptr) :: Ap = c_null_ptr
+    type(c_ptr) :: Ai = c_null_ptr
+    type(c_ptr) :: Ax = c_null_ptr
+
+    type(c_ptr) :: control = c_null_ptr
+    real(kind=c_double), dimension(90) :: info
+
+    real(kind=c_double), dimension(:), pointer :: sol => null()
+
+  end type sparse_mat
+
+
+  public declare, build, &
+         solve  , erase, &
+         init   , finalize
+
+
+  contains
+
+
+  subroutine init()
+    implicit none
+
+    print *, 'empty initialization'
+
+  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
+
+    ! Ap ~ system%cc_adj_dof
+    ! Ai ~ system%adj_dof
+
+    integer(kind=c_int), dimension(:), pointer :: Ap, Ai
+    real(kind=c_double), dimension(:), pointer :: control
+
+    integer(kind=4) :: idx, i
+
+    if( is_sym ) then
+      print *, 'ERROR : umfpack available only for non-symmetrical systems'
+    end if
+
+    allocate( Ap(nb_dofs+1  ) )
+    allocate( Ai(nb_non_zero) )
+
+    Ap(:) = 0
+    Ai(:) = j_indices(:)-1
+
+    ! assume i_indices is ordered in an increasing order
+
+    ! counting non zero on a line
+    idx = 0
+
+    do i = 2, nb_non_zero
+      if( idx /= i_indices(i) ) idx = i_indices(i)
+      Ap(idx+1) = Ap(idx+1) + 1
+    end do
+
+    ! cumulating count
+    Ap(1) = 1
+    do i = 1, nb_dofs
+      Ap(i+1) = Ap(i) + Ap(i+1)
+    end do
+    Ap(1) = 0
+
+    matrix%Ap = c_loc( Ap(1) )
+    matrix%Ai = c_loc( Ai(1) )
+
+    matrix%n  = nb_dofs
+
+    allocate( control(20) )
+    control(:) = 0.d0
+    matrix%control  = c_loc( control(1) )
+
+
+    info = umfpack_symbolic(nb_dofs,nb_dofs,matrix%Ap,matrix%Ai,matrix%Ax,matrix%symbolic,matrix%control,matrix%info)
+
+    allocate( matrix%sol(nb_dofs) )
+    matrix%sol = 0.d0
+
+  end subroutine
+
+
+  subroutine build(matrix, val, info)
+    implicit none
+    type(sparse_mat) :: matrix
+    real(kind=8)   , dimension(:), pointer :: val
+    integer(kind=4) :: info
+
+    matrix%Ax = c_loc( val(1) )
+
+    ! Factorize
+    info = umfpack_numeric(matrix%Ap,matrix%Ai,matrix%Ax,matrix%symbolic,matrix%numeric,matrix%control,matrix%info)
+
+  end subroutine
+
+  subroutine solve(matrix, rhs, info)
+    implicit none
+
+    type(sparse_mat) :: matrix
+
+    real(kind=8)   , dimension(:), pointer :: rhs
+    integer(kind=4) :: info
+
+    ! Solve
+    info = umfpack_solve( matrix%Ap, matrix%Ai, matrix%Ax, c_loc(matrix%sol(1)), c_loc(rhs(1)), &
+                          matrix%numeric, matrix%control, matrix%info )
+
+    rhs = matrix%sol
+
+  end subroutine
+
+  subroutine erase(matrix)
+    implicit none
+    type(sparse_mat) :: matrix
+
+    integer(kind=c_int), dimension(:), pointer :: Ap, Ai
+    real(kind=c_double), dimension(:), pointer :: Ax, sol, control
+
+    call umfpack_free_symbolic( matrix%symbolic )
+    call umfpack_free_numeric(  matrix%numeric  )
+
+    matrix%symbolic = c_null_ptr
+    matrix%numeric  = c_null_ptr
+
+    call c_f_pointer(matrix%Ap, ap, (/matrix%n/))
+
+    call c_f_pointer(matrix%Ai, ai, (/ap(matrix%n+1)/))
+
+    if(associated(ap)) deallocate(ap)
+    if(associated(ai)) deallocate(ai)
+
+    matrix%Ap = c_null_ptr
+    matrix%Ai = c_null_ptr
+    matrix%Ax = c_null_ptr
+
+    call c_f_pointer(matrix%control, control, (/20/))
+    deallocate(control)
+    matrix%control = c_null_ptr
+
+    if( associated(matrix%sol) ) deallocate(matrix%sol)
+    matrix%sol => null()
+
+  end subroutine
+
+  subroutine finalize()
+    implicit none
+
+    print *, 'empty finalization'
+
+  end subroutine
+
+end module

systems_in_parallel/umfpack_cbinding.c

+
+#include "stdio.h"
+
+#include "umfpack_cbinding.h"
+
+
+int c_symbolic(int nrow, int ncol, int * ap, int * ai, double * ax, void ** symbolic, double ** control, double * info)
+{
+  int i;
+  umfpack_di_defaults(*control);
+  return umfpack_di_symbolic (nrow, ncol, ap, ai, ax, symbolic, *control, info) ;
+}
+
+int c_numeric(int * ap, int * ai, double * ax, void ** symbolic, void ** numeric, double ** control, double * info)
+{
+  return umfpack_di_numeric(ap, ai, ax, *symbolic, numeric, *control, info) ;
+}
+
+int c_solve(int * ap, int * ai, double * ax, double ** x, double ** b, void ** numeric, double ** control, double * info)
+{
+  return umfpack_di_solve(UMFPACK_Aat, ap, ai, ax, *x, *b, *numeric, *control, info) ;
+}
+
+void c_free_symbolic(void ** symbolic)
+{
+  return umfpack_di_free_symbolic(symbolic);
+}
+
+void c_free_numeric(void ** numeric)
+{
+  return umfpack_di_free_numeric(numeric);
+}
+

systems_in_parallel/umfpack_cbinding.h

+
+#ifndef umfpack_cbind
+#define umfpack_cbind
+
+#include "umfpack.h"
+
+int c_symbolic(int nrow, int ncol, int * ap, int * ai, double * ax, void ** symbolic, double ** control, double * info);
+
+int c_numeric(int * ap, int * ai, double * ax, void ** symbolic, void ** numeric, double ** control, double * info);
+
+int c_solve(int * ap, int * ai, double * ax, double ** x, double ** b, void ** numeric, double ** control, double * info);
+
+void c_free_symbolic(void ** symbolic);
+
+void c_free_numeric(void ** numeric);
+
+#endif // umfpack_cbind
+

systems_in_parallel/umfpack_test.c

+
+#include "umfpack.h"
+#include "stdio.h"
+
+int n=5;
+int Ap[]={0,3,6,9,10,12} ;
+int Ai[]={0, 1, 2, 0, 2, 4, 1, 2, 3, 2, 1, 4};
+double Ax [ ] = {2., 3., 4., 3.,-3., 6.,-1., 1., 2., 2., 4., 1.};
+double b [ ] = {20., 24., 9., 6., 13.} ;
+double x [5] ;
+
+int main (void)
+{
+    double *null = (double *) NULL ;
+    int i ;
+    void *Symbolic, *Numeric ;
+    i = umfpack_di_symbolic (n, n, Ap, Ai, null, &Symbolic, null, null) ;
+    printf("result of symbolic %d\n",i);
+    (void) umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, null, null) ;
+    umfpack_di_free_symbolic (&Symbolic) ;
+    (void) umfpack_di_solve (UMFPACK_Aat, Ap, Ai, Ax, x, b, Numeric, null, null) ;
+    umfpack_di_free_numeric (&Numeric) ;
+    for (i = 0 ; i < n ; i++) printf ("x [%d] = %g\n", i, x [i]) ;
+    return (0) ;
+}
+