From db6500dcfb6c82b1630b7a6d35067cd1ee90cf78 Mon Sep 17 00:00:00 2001
From: William Daniau <william.daniau@femto-st.fr>
Date: Fri, 29 Jun 2012 10:22:41 +0200
Subject: [PATCH] First draft for a multi resolution

---
 fvn_sparse/fvn_sparse.f90 | 83 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 83 insertions(+)

diff --git a/fvn_sparse/fvn_sparse.f90 b/fvn_sparse/fvn_sparse.f90
index fb88415..145fdce 100644
--- a/fvn_sparse/fvn_sparse.f90
+++ b/fvn_sparse/fvn_sparse.f90
@@ -345,6 +345,89 @@ deallocate(A)
 end subroutine
 
 
+subroutine fvn_di_sparse_solve_multi(n,nz,T,Ti,Tj,m,B,x,status,det)
+implicit none
+integer(kind=sp_kind), intent(in) :: n,nz,m
+real(kind=dp_kind),dimension(nz),intent(in) :: T
+integer(kind=sp_kind),dimension(nz),intent(in)  :: Ti,Tj
+real(kind=dp_kind),dimension(n,m),intent(in) :: B
+real(kind=dp_kind),dimension(n,m),intent(out) :: x
+integer(kind=sp_kind), intent(out) :: status
+real(kind=dp_kind), dimension(2), optional, intent(out) :: det
+
+real(kind=dp_kind),dimension(:),allocatable :: A
+integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai
+!integer(kind=dp_kind) :: symbolic,numeric
+integer(kind=sp_kind),dimension(2) :: symbolic,numeric
+! As symbolic and numeric are used to store a C pointer, it is necessary to
+! still use an integer(kind=dp_kind) for 64bits machines
+! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric
+real(kind=dp_kind),dimension(90) :: info
+real(kind=dp_kind),dimension(20) :: control
+integer(kind=sp_kind) :: sys
+integer(kind=sp_kind) :: i
+
+status=0
+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,Ti,Tj,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 : ",trim(umfpack_return_code(status))
+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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
+    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 : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
+    status=info(1)
+endif
+
+! free the C symbolic pointer
+call umfpack_di_free_symbolic (symbolic)
+
+! if parameter det is present, the determinant of the matrix is calculated
+if (present(det) ) then
+  call umfpack_di_get_determinant(det(1),det(2),numeric,info,status)  
+  ! info(1) should be zero
+  if (info(1) /= 0) then
+      if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning
+          write(*,*) "Problem during sparse determinant : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
+      endif
+      status=info(1)
+  endif
+endif
+
+sys=0
+! sys may be used to define type of solving -> see umfpack.h
+do i=1,m
+    ! Solving
+    call umfpack_di_solve (sys, Ap, Ai, A, x(:,i), B(:,i), numeric, control, info)
+    ! info(1) should be zero
+    if (info(1) /= 0) then
+        write(*,*) "Problem during solving : ",trim(umfpack_return_code(int(info(1),kind=sp_kind)))
+        status=info(1)
+    endif
+end do
+! free the C numeric pointer
+call umfpack_di_free_numeric (numeric)
+
+deallocate(A)
+end subroutine
 
 
 
-- 
2.16.4