module fvn_sparse use fvn_common implicit none ! Sparse solving interface fvn_sparse_solve module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve end interface fvn_sparse_solve contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! SPARSE RESOLUTION ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Sparse resolution is done by interfaçing Tim Davi's UMFPACK ! http://www.cise.ufl.edu/research/sparse/SuiteSparse/ ! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig ! ! Solve Ax=B using UMFPACK ! ! Where A is a sparse matrix given in its triplet form ! T -> non zero elements ! Ti,Tj -> row and column index (1-based) of the given elt ! n : rank of matrix A ! nz : number of non zero elts ! ! fvn_*_sparse_solve ! * = zl : double complex + integer(8) ! * = zi : double complex + integer(4) ! subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status) implicit none integer(8), intent(in) :: n,nz complex(8),dimension(nz),intent(in) :: T integer(8),dimension(nz),intent(in) :: Ti,Tj complex(8),dimension(n),intent(in) :: B complex(8),dimension(n),intent(out) :: x integer(8), intent(out) :: status integer(8),dimension(:),allocatable :: wTi,wTj real(8),dimension(:),allocatable :: Tx,Tz real(8),dimension(:),allocatable :: Ax,Az integer(8),dimension(:),allocatable :: Ap,Ai integer(8) :: symbolic,numeric real(8),dimension(:),allocatable :: xx,xz,bx,bz real(8),dimension(90) :: info real(8),dimension(20) :: control integer(8) :: sys status=0 ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation ! Tx and Tz are the real and imaginary parts of T allocate(wTi(nz),wTj(nz)) allocate(Tx(nz),Tz(nz)) Tx=dble(T) Tz=aimag(T) wTi=Ti-1 wTj=Tj-1 allocate(Ax(nz),Az(nz)) allocate(Ap(n+1),Ai(nz)) ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az call umfpack_zl_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status) ! if status is not zero a problem has occured if (status /= 0) then write(*,*) "Problem during umfpack_zl_triplet_to_col" endif ! Define defaults control values call umfpack_zl_defaults(control) ! Symbolic analysis call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during symbolic analysis" status=info(1) endif ! Numerical factorization call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during numerical factorization" status=info(1) endif ! free the C symbolic pointer call umfpack_zl_free_symbolic (symbolic) allocate(bx(n),bz(n),xx(n),xz(n)) bx=dble(B) bz=aimag(B) sys=0 ! sys may be used to define type of solving -> see umfpack.h ! Solving call umfpack_zl_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during solving" status=info(1) endif ! free the C numeric pointer call umfpack_zl_free_numeric (numeric) x=dcmplx(xx,xz) deallocate(bx,bz,xx,xz) deallocate(Ax,Az) deallocate(Tx,Tz) deallocate(wTi,wTj) end subroutine subroutine fvn_zi_sparse_solve(n,nz,T,Ti,Tj,B,x,status) implicit none integer(4), intent(in) :: n,nz complex(8),dimension(nz),intent(in) :: T integer(4),dimension(nz),intent(in) :: Ti,Tj complex(8),dimension(n),intent(in) :: B complex(8),dimension(n),intent(out) :: x integer(4), intent(out) :: status integer(4),dimension(:),allocatable :: wTi,wTj real(8),dimension(:),allocatable :: Tx,Tz real(8),dimension(:),allocatable :: Ax,Az integer(4),dimension(:),allocatable :: Ap,Ai !integer(8) :: symbolic,numeric integer(4),dimension(2) :: symbolic,numeric ! As symbolic and numeric are used to store a C pointer, it is necessary to ! still use an integer(8) for 64bits machines ! An other possibility : integer(4),dimension(2) :: symbolic,numeric real(8),dimension(:),allocatable :: xx,xz,bx,bz real(8),dimension(90) :: info real(8),dimension(20) :: control integer(4) :: sys status=0 ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation ! Tx and Tz are the real and imaginary parts of T allocate(wTi(nz),wTj(nz)) allocate(Tx(nz),Tz(nz)) Tx=dble(T) Tz=aimag(T) wTi=Ti-1 wTj=Tj-1 allocate(Ax(nz),Az(nz)) allocate(Ap(n+1),Ai(nz)) ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az call umfpack_zi_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status) ! if status is not zero a problem has occured if (status /= 0) then write(*,*) "Problem during umfpack_zl_triplet_to_col" endif ! Define defaults control values call umfpack_zi_defaults(control) ! Symbolic analysis call umfpack_zi_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during symbolic analysis" status=info(1) endif ! Numerical factorization call umfpack_zi_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during numerical factorization" status=info(1) endif ! free the C symbolic pointer call umfpack_zi_free_symbolic (symbolic) allocate(bx(n),bz(n),xx(n),xz(n)) bx=dble(B) bz=aimag(B) sys=0 ! sys may be used to define type of solving -> see umfpack.h ! Solving call umfpack_zi_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during solving" status=info(1) endif ! free the C numeric pointer call umfpack_zi_free_numeric (numeric) x=dcmplx(xx,xz) deallocate(bx,bz,xx,xz) deallocate(Ax,Az) deallocate(Tx,Tz) deallocate(wTi,wTj) end subroutine subroutine fvn_dl_sparse_solve(n,nz,T,Ti,Tj,B,x,status) implicit none integer(8), intent(in) :: n,nz real(8),dimension(nz),intent(in) :: T integer(8),dimension(nz),intent(in) :: Ti,Tj real(8),dimension(n),intent(in) :: B real(8),dimension(n),intent(out) :: x integer(8), intent(out) :: status integer(8),dimension(:),allocatable :: wTi,wTj real(8),dimension(:),allocatable :: A integer(8),dimension(:),allocatable :: Ap,Ai !integer(8) :: symbolic,numeric integer(8) :: symbolic,numeric real(8),dimension(90) :: info real(8),dimension(20) :: control integer(8) :: sys status=0 ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation allocate(wTi(nz),wTj(nz)) wTi=Ti-1 wTj=Tj-1 allocate(A(nz)) allocate(Ap(n+1),Ai(nz)) ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az call umfpack_dl_triplet_to_col(n,n,nz,wTi,wTj,T,Ap,Ai,A,status) ! if status is not zero a problem has occured if (status /= 0) then write(*,*) "Problem during umfpack_dl_triplet_to_col" endif ! Define defaults control values call umfpack_dl_defaults(control) ! Symbolic analysis call umfpack_dl_symbolic(n,n,Ap,Ai,A,symbolic, control, info) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during symbolic analysis" status=info(1) endif ! Numerical factorization call umfpack_dl_numeric (Ap, Ai, A, symbolic, numeric, control, info) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during numerical factorization" status=info(1) endif ! free the C symbolic pointer call umfpack_dl_free_symbolic (symbolic) sys=0 ! sys may be used to define type of solving -> see umfpack.h ! Solving call umfpack_dl_solve (sys, Ap, Ai, A, x, B, numeric, control, info) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during solving" status=info(1) endif ! free the C numeric pointer call umfpack_dl_free_numeric (numeric) deallocate(A) deallocate(wTi,wTj) end subroutine subroutine fvn_di_sparse_solve(n,nz,T,Ti,Tj,B,x,status) implicit none integer(4), intent(in) :: n,nz real(8),dimension(nz),intent(in) :: T integer(4),dimension(nz),intent(in) :: Ti,Tj real(8),dimension(n),intent(in) :: B real(8),dimension(n),intent(out) :: x integer(4), intent(out) :: status integer(4),dimension(:),allocatable :: wTi,wTj real(8),dimension(:),allocatable :: A integer(4),dimension(:),allocatable :: Ap,Ai !integer(8) :: symbolic,numeric integer(4),dimension(2) :: symbolic,numeric ! As symbolic and numeric are used to store a C pointer, it is necessary to ! still use an integer(8) for 64bits machines ! An other possibility : integer(4),dimension(2) :: symbolic,numeric real(8),dimension(90) :: info real(8),dimension(20) :: control integer(4) :: sys status=0 ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation allocate(wTi(nz),wTj(nz)) wTi=Ti-1 wTj=Tj-1 allocate(A(nz)) allocate(Ap(n+1),Ai(nz)) ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az call umfpack_di_triplet_to_col(n,n,nz,wTi,wTj,T,Ap,Ai,A,status) ! if status is not zero a problem has occured if (status /= 0) then write(*,*) "Problem during umfpack_di_triplet_to_col" endif ! Define defaults control values call umfpack_di_defaults(control) ! Symbolic analysis call umfpack_di_symbolic(n,n,Ap,Ai,A,symbolic, control, info) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during symbolic analysis" status=info(1) endif ! Numerical factorization call umfpack_di_numeric (Ap, Ai, A, symbolic, numeric, control, info) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during numerical factorization" status=info(1) endif ! free the C symbolic pointer call umfpack_di_free_symbolic (symbolic) sys=0 ! sys may be used to define type of solving -> see umfpack.h ! Solving call umfpack_di_solve (sys, Ap, Ai, A, x, B, numeric, control, info) ! info(1) should be zero if (info(1) /= 0) then write(*,*) "Problem during solving" status=info(1) endif ! free the C numeric pointer call umfpack_di_free_numeric (numeric) deallocate(A) deallocate(wTi,wTj) end subroutine end module fvn_sparse