Commit aee8aaaf183bdbeeb66f800e26aec6d8e9474150
1 parent
9c285563c3
Exists in
master
and in
2 other branches
Added a function to get meaning of umfpack returned code
git-svn-id: https://lxsd.femto-st.fr/svn/fvn@72 b657c933-2333-4658-acf2-d3c7c2708721
Showing 1 changed file with 40 additions and 0 deletions Inline Diff
fvn_sparse/fvn_sparse.f90
| module fvn_sparse | 1 | 1 | module fvn_sparse | |
| use fvn_common | 2 | 2 | use fvn_common | |
| implicit none | 3 | 3 | implicit none | |
| 4 | 4 | |||
| ! Sparse solving | 5 | 5 | ! Sparse solving | |
| interface fvn_sparse_solve | 6 | 6 | interface fvn_sparse_solve | |
| module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve | 7 | 7 | module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve | |
| end interface fvn_sparse_solve | 8 | 8 | end interface fvn_sparse_solve | |
| 9 | 9 | |||
| interface fvn_sparse_det | 10 | 10 | interface fvn_sparse_det | |
| module procedure fvn_zl_sparse_det,fvn_zi_sparse_det,fvn_dl_sparse_det,fvn_di_sparse_det | 11 | 11 | module procedure fvn_zl_sparse_det,fvn_zi_sparse_det,fvn_dl_sparse_det,fvn_di_sparse_det | |
| end interface fvn_sparse_det | 12 | 12 | end interface fvn_sparse_det | |
| contains | 13 | 13 | contains | |
| 14 | ||||
| 15 | function umfpack_return_code(c) | |||
| 16 | implicit none | |||
| 17 | integer(4), intent(in) :: c | |||
| 18 | character(len=80) :: umfpack_return_code | |||
| 19 | select case(c) | |||
| 20 | case(0) | |||
| 21 | umfpack_return_code="UMFPACK_OK" | |||
| 22 | case(1) | |||
| 23 | umfpack_return_code="UMFPACK_WARNING_singular_matrix" | |||
| 24 | case(2) | |||
| 25 | umfpack_return_code="UMFPACK_WARNING_determinant_underflow" | |||
| 26 | case(3) | |||
| 27 | umfpack_return_code="UMFPACK_WARNING_determinant_overflow" | |||
| 28 | case(-1) | |||
| 29 | umfpack_return_code="UMFPACK_ERROR_out_of_memory" | |||
| 30 | case(-3) | |||
| 31 | umfpack_return_code="UMFPACK_ERROR_invalid_Numeric_object" | |||
| 32 | case(-4) | |||
| 33 | umfpack_return_code="UMFPACK_ERROR_invalid_Symbolic_object" | |||
| 34 | case(-5) | |||
| 35 | umfpack_return_code="UMFPACK_ERROR_argument_missing" | |||
| 36 | case(-6) | |||
| 37 | umfpack_return_code="UMFPACK_ERROR_n_nonpositive" | |||
| 38 | case(-8) | |||
| 39 | umfpack_return_code="UMFPACK_ERROR_invalid_matrix" | |||
| 40 | case(-11) | |||
| 41 | umfpack_return_code="UMFPACK_ERROR_different_pattern" | |||
| 42 | case(-13) | |||
| 43 | umfpack_return_code="UMFPACK_ERROR_invalid_system" | |||
| 44 | case(-15) | |||
| 45 | umfpack_return_code="UMFPACK_ERROR_invalid_permutation" | |||
| 46 | case(-911) | |||
| 47 | umfpack_return_code="UMFPACK_ERROR_internal_error" | |||
| 48 | case(-17) | |||
| 49 | umfpack_return_code="UMFPACK_ERROR_file_IO" | |||
| 50 | case default | |||
| 51 | umfpack_return_code="Unknown return code" | |||
| 52 | end select | |||
| 53 | end function | |||
| !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | 14 | 54 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | |
| ! | 15 | 55 | ! | |
| ! SPARSE RESOLUTION | 16 | 56 | ! SPARSE RESOLUTION | |
| ! | 17 | 57 | ! | |
| !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | 18 | 58 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | |
| ! | 19 | 59 | ! | |
| ! Sparse resolution is done by interfaçing Tim Davi's UMFPACK | 20 | 60 | ! Sparse resolution is done by interfaçing Tim Davi's UMFPACK | |
| ! http://www.cise.ufl.edu/research/sparse/SuiteSparse/ | 21 | 61 | ! http://www.cise.ufl.edu/research/sparse/SuiteSparse/ | |
| ! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig | 22 | 62 | ! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig | |
| ! | 23 | 63 | ! | |
| ! Solve Ax=B using UMFPACK | 24 | 64 | ! Solve Ax=B using UMFPACK | |
| ! | 25 | 65 | ! | |
| ! Where A is a sparse matrix given in its triplet form | 26 | 66 | ! Where A is a sparse matrix given in its triplet form | |
| ! T -> non zero elements | 27 | 67 | ! T -> non zero elements | |
| ! Ti,Tj -> row and column index (1-based) of the given elt | 28 | 68 | ! Ti,Tj -> row and column index (1-based) of the given elt | |
| ! n : rank of matrix A | 29 | 69 | ! n : rank of matrix A | |
| ! nz : number of non zero elts | 30 | 70 | ! nz : number of non zero elts | |
| ! | 31 | 71 | ! | |
| ! fvn_*_sparse_solve | 32 | 72 | ! fvn_*_sparse_solve | |
| ! * = zl : double complex + integer(kind=dp_kind) | 33 | 73 | ! * = zl : double complex + integer(kind=dp_kind) | |
| ! * = zi : double complex + integer(kind=sp_kind) | 34 | 74 | ! * = zi : double complex + integer(kind=sp_kind) | |
| ! | 35 | 75 | ! | |
| subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det) | 36 | 76 | subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det) | |
| implicit none | 37 | 77 | implicit none | |
| integer(kind=dp_kind), intent(in) :: n,nz | 38 | 78 | integer(kind=dp_kind), intent(in) :: n,nz | |
| complex(kind=dp_kind),dimension(nz),intent(in) :: T | 39 | 79 | complex(kind=dp_kind),dimension(nz),intent(in) :: T | |
| integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj | 40 | 80 | integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj | |
| complex(kind=dp_kind),dimension(n),intent(in) :: B | 41 | 81 | complex(kind=dp_kind),dimension(n),intent(in) :: B | |
| complex(kind=dp_kind),dimension(n),intent(out) :: x | 42 | 82 | complex(kind=dp_kind),dimension(n),intent(out) :: x | |
| integer(kind=dp_kind), intent(out) :: status | 43 | 83 | integer(kind=dp_kind), intent(out) :: status | |
| real(kind=dp_kind), dimension(3), optional, intent(out) :: det | 44 | 84 | real(kind=dp_kind), dimension(3), optional, intent(out) :: det | |
| 45 | 85 | |||
| integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj | 46 | 86 | integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj | |
| real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz | 47 | 87 | real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz | |
| real(kind=dp_kind),dimension(:),allocatable :: Ax,Az | 48 | 88 | real(kind=dp_kind),dimension(:),allocatable :: Ax,Az | |
| integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai | 49 | 89 | integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai | |
| integer(kind=dp_kind) :: symbolic,numeric | 50 | 90 | integer(kind=dp_kind) :: symbolic,numeric | |
| real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz | 51 | 91 | real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz | |
| real(kind=dp_kind),dimension(90) :: info | 52 | 92 | real(kind=dp_kind),dimension(90) :: info | |
| real(kind=dp_kind),dimension(20) :: control | 53 | 93 | real(kind=dp_kind),dimension(20) :: control | |
| integer(kind=dp_kind) :: sys | 54 | 94 | integer(kind=dp_kind) :: sys | |
| 55 | 95 | |||
| 56 | 96 | |||
| status=0 | 57 | 97 | status=0 | |
| 58 | 98 | |||
| ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation | 59 | 99 | ! 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 | 60 | 100 | ! Tx and Tz are the real and imaginary parts of T | |
| allocate(wTi(nz),wTj(nz)) | 61 | 101 | allocate(wTi(nz),wTj(nz)) | |
| allocate(Tx(nz),Tz(nz)) | 62 | 102 | allocate(Tx(nz),Tz(nz)) | |
| Tx=dble(T) | 63 | 103 | Tx=dble(T) | |
| Tz=aimag(T) | 64 | 104 | Tz=aimag(T) | |
| wTi=Ti-1 | 65 | 105 | wTi=Ti-1 | |
| wTj=Tj-1 | 66 | 106 | wTj=Tj-1 | |
| allocate(Ax(nz),Az(nz)) | 67 | 107 | allocate(Ax(nz),Az(nz)) | |
| allocate(Ap(n+1),Ai(nz)) | 68 | 108 | allocate(Ap(n+1),Ai(nz)) | |
| 69 | 109 | |||
| ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az | 70 | 110 | ! 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) | 71 | 111 | 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 | 72 | 112 | ! if status is not zero a problem has occured | |
| if (status /= 0) then | 73 | 113 | if (status /= 0) then | |
| write(*,*) "Problem during umfpack_zl_triplet_to_col" | 74 | 114 | write(*,*) "Problem during umfpack_zl_triplet_to_col" | |
| endif | 75 | 115 | endif | |
| 76 | 116 | |||
| ! Define defaults control values | 77 | 117 | ! Define defaults control values | |
| call umfpack_zl_defaults(control) | 78 | 118 | call umfpack_zl_defaults(control) | |
| 79 | 119 | |||
| ! Symbolic analysis | 80 | 120 | ! Symbolic analysis | |
| call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) | 81 | 121 | call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) | |
| ! info(1) should be zero | 82 | 122 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 83 | 123 | if (info(1) /= 0) then | |
| write(*,*) "Problem during symbolic analysis" | 84 | 124 | write(*,*) "Problem during symbolic analysis" | |
| status=info(1) | 85 | 125 | status=info(1) | |
| endif | 86 | 126 | endif | |
| 87 | 127 | |||
| ! Numerical factorization | 88 | 128 | ! Numerical factorization | |
| call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) | 89 | 129 | call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) | |
| ! info(1) should be zero | 90 | 130 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 91 | 131 | if (info(1) /= 0) then | |
| write(*,*) "Problem during numerical factorization" | 92 | 132 | write(*,*) "Problem during numerical factorization" | |
| status=info(1) | 93 | 133 | status=info(1) | |
| endif | 94 | 134 | endif | |
| 95 | 135 | |||
| ! free the C symbolic pointer | 96 | 136 | ! free the C symbolic pointer | |
| call umfpack_zl_free_symbolic (symbolic) | 97 | 137 | call umfpack_zl_free_symbolic (symbolic) | |
| 98 | 138 | |||
| ! if parameter det is present, the determinant of the matrix is calculated | 99 | 139 | ! if parameter det is present, the determinant of the matrix is calculated | |
| if (present(det) ) then | 100 | 140 | if (present(det) ) then | |
| call umfpack_zl_get_determinant(det(1),det(2),det(3),numeric,info,status) | 101 | 141 | call umfpack_zl_get_determinant(det(1),det(2),det(3),numeric,info,status) | |
| ! info(1) should be zero | 102 | 142 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 103 | 143 | if (info(1) /= 0) then | |
| if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning | 104 | 144 | if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning | |
| write(*,*) "Problem during sparse determinant, returned code : ",info(1) | 105 | 145 | write(*,*) "Problem during sparse determinant, returned code : ",info(1) | |
| endif | 106 | 146 | endif | |
| status=info(1) | 107 | 147 | status=info(1) | |
| endif | 108 | 148 | endif | |
| endif | 109 | 149 | endif | |
| 110 | 150 | |||
| 111 | 151 | |||
| 112 | 152 | |||
| allocate(bx(n),bz(n),xx(n),xz(n)) | 113 | 153 | allocate(bx(n),bz(n),xx(n),xz(n)) | |
| bx=dble(B) | 114 | 154 | bx=dble(B) | |
| bz=aimag(B) | 115 | 155 | bz=aimag(B) | |
| sys=0 | 116 | 156 | sys=0 | |
| ! sys may be used to define type of solving -> see umfpack.h | 117 | 157 | ! sys may be used to define type of solving -> see umfpack.h | |
| 118 | 158 | |||
| ! Solving | 119 | 159 | ! Solving | |
| call umfpack_zl_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info) | 120 | 160 | call umfpack_zl_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info) | |
| ! info(1) should be zero | 121 | 161 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 122 | 162 | if (info(1) /= 0) then | |
| write(*,*) "Problem during solving" | 123 | 163 | write(*,*) "Problem during solving" | |
| status=info(1) | 124 | 164 | status=info(1) | |
| endif | 125 | 165 | endif | |
| 126 | 166 | |||
| 127 | 167 | |||
| ! free the C numeric pointer | 128 | 168 | ! free the C numeric pointer | |
| call umfpack_zl_free_numeric (numeric) | 129 | 169 | call umfpack_zl_free_numeric (numeric) | |
| 130 | 170 | |||
| x=cmplx(xx,xz,dp_kind) | 131 | 171 | x=cmplx(xx,xz,dp_kind) | |
| 132 | 172 | |||
| deallocate(bx,bz,xx,xz) | 133 | 173 | deallocate(bx,bz,xx,xz) | |
| deallocate(Ax,Az) | 134 | 174 | deallocate(Ax,Az) | |
| deallocate(Tx,Tz) | 135 | 175 | deallocate(Tx,Tz) | |
| deallocate(wTi,wTj) | 136 | 176 | deallocate(wTi,wTj) | |
| end subroutine | 137 | 177 | end subroutine | |
| 138 | 178 | |||
| 139 | 179 | |||
| 140 | 180 | |||
| 141 | 181 | |||
| 142 | 182 | |||
| subroutine fvn_zi_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det) | 143 | 183 | subroutine fvn_zi_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det) | |
| implicit none | 144 | 184 | implicit none | |
| integer(kind=sp_kind), intent(in) :: n,nz | 145 | 185 | integer(kind=sp_kind), intent(in) :: n,nz | |
| complex(kind=dp_kind),dimension(nz),intent(in) :: T | 146 | 186 | complex(kind=dp_kind),dimension(nz),intent(in) :: T | |
| integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj | 147 | 187 | integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj | |
| complex(kind=dp_kind),dimension(n),intent(in) :: B | 148 | 188 | complex(kind=dp_kind),dimension(n),intent(in) :: B | |
| complex(kind=dp_kind),dimension(n),intent(out) :: x | 149 | 189 | complex(kind=dp_kind),dimension(n),intent(out) :: x | |
| integer(kind=sp_kind), intent(out) :: status | 150 | 190 | integer(kind=sp_kind), intent(out) :: status | |
| real(kind=dp_kind), dimension(3), optional, intent(out) :: det | 151 | 191 | real(kind=dp_kind), dimension(3), optional, intent(out) :: det | |
| 152 | 192 | |||
| integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj | 153 | 193 | integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj | |
| real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz | 154 | 194 | real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz | |
| real(kind=dp_kind),dimension(:),allocatable :: Ax,Az | 155 | 195 | real(kind=dp_kind),dimension(:),allocatable :: Ax,Az | |
| integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai | 156 | 196 | integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai | |
| !integer(kind=dp_kind) :: symbolic,numeric | 157 | 197 | !integer(kind=dp_kind) :: symbolic,numeric | |
| integer(kind=sp_kind),dimension(2) :: symbolic,numeric | 158 | 198 | integer(kind=sp_kind),dimension(2) :: symbolic,numeric | |
| ! As symbolic and numeric are used to store a C pointer, it is necessary to | 159 | 199 | ! 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 | 160 | 200 | ! still use an integer(kind=dp_kind) for 64bits machines | |
| ! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric | 161 | 201 | ! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric | |
| real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz | 162 | 202 | real(kind=dp_kind),dimension(:),allocatable :: xx,xz,bx,bz | |
| real(kind=dp_kind),dimension(90) :: info | 163 | 203 | real(kind=dp_kind),dimension(90) :: info | |
| real(kind=dp_kind),dimension(20) :: control | 164 | 204 | real(kind=dp_kind),dimension(20) :: control | |
| integer(kind=sp_kind) :: sys | 165 | 205 | integer(kind=sp_kind) :: sys | |
| 166 | 206 | |||
| status=0 | 167 | 207 | status=0 | |
| ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation | 168 | 208 | ! 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 | 169 | 209 | ! Tx and Tz are the real and imaginary parts of T | |
| allocate(wTi(nz),wTj(nz)) | 170 | 210 | allocate(wTi(nz),wTj(nz)) | |
| allocate(Tx(nz),Tz(nz)) | 171 | 211 | allocate(Tx(nz),Tz(nz)) | |
| Tx=dble(T) | 172 | 212 | Tx=dble(T) | |
| Tz=aimag(T) | 173 | 213 | Tz=aimag(T) | |
| wTi=Ti-1 | 174 | 214 | wTi=Ti-1 | |
| wTj=Tj-1 | 175 | 215 | wTj=Tj-1 | |
| allocate(Ax(nz),Az(nz)) | 176 | 216 | allocate(Ax(nz),Az(nz)) | |
| allocate(Ap(n+1),Ai(nz)) | 177 | 217 | allocate(Ap(n+1),Ai(nz)) | |
| 178 | 218 | |||
| ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az | 179 | 219 | ! 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) | 180 | 220 | 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 | 181 | 221 | ! if status is not zero a problem has occured | |
| if (status /= 0) then | 182 | 222 | if (status /= 0) then | |
| write(*,*) "Problem during umfpack_zl_triplet_to_col" | 183 | 223 | write(*,*) "Problem during umfpack_zl_triplet_to_col" | |
| endif | 184 | 224 | endif | |
| 185 | 225 | |||
| ! Define defaults control values | 186 | 226 | ! Define defaults control values | |
| call umfpack_zi_defaults(control) | 187 | 227 | call umfpack_zi_defaults(control) | |
| 188 | 228 | |||
| ! Symbolic analysis | 189 | 229 | ! Symbolic analysis | |
| call umfpack_zi_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) | 190 | 230 | call umfpack_zi_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) | |
| ! info(1) should be zero | 191 | 231 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 192 | 232 | if (info(1) /= 0) then | |
| write(*,*) "Problem during symbolic analysis" | 193 | 233 | write(*,*) "Problem during symbolic analysis" | |
| status=info(1) | 194 | 234 | status=info(1) | |
| endif | 195 | 235 | endif | |
| 196 | 236 | |||
| ! Numerical factorization | 197 | 237 | ! Numerical factorization | |
| call umfpack_zi_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) | 198 | 238 | call umfpack_zi_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) | |
| ! info(1) should be zero | 199 | 239 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 200 | 240 | if (info(1) /= 0) then | |
| write(*,*) "Problem during numerical factorization" | 201 | 241 | write(*,*) "Problem during numerical factorization" | |
| status=info(1) | 202 | 242 | status=info(1) | |
| endif | 203 | 243 | endif | |
| 204 | 244 | |||
| ! free the C symbolic pointer | 205 | 245 | ! free the C symbolic pointer | |
| call umfpack_zi_free_symbolic (symbolic) | 206 | 246 | call umfpack_zi_free_symbolic (symbolic) | |
| 207 | 247 | |||
| ! if parameter det is present, the determinant of the matrix is calculated | 208 | 248 | ! if parameter det is present, the determinant of the matrix is calculated | |
| if (present(det) ) then | 209 | 249 | if (present(det) ) then | |
| call umfpack_zi_get_determinant(det(1),det(2),det(3),numeric,info,status) | 210 | 250 | call umfpack_zi_get_determinant(det(1),det(2),det(3),numeric,info,status) | |
| ! info(1) should be zero | 211 | 251 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 212 | 252 | if (info(1) /= 0) then | |
| if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning | 213 | 253 | if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning | |
| write(*,*) "Problem during sparse determinant, returned code : ",info(1) | 214 | 254 | write(*,*) "Problem during sparse determinant, returned code : ",info(1) | |
| endif | 215 | 255 | endif | |
| status=info(1) | 216 | 256 | status=info(1) | |
| endif | 217 | 257 | endif | |
| endif | 218 | 258 | endif | |
| 219 | 259 | |||
| 220 | 260 | |||
| 221 | 261 | |||
| 222 | 262 | |||
| allocate(bx(n),bz(n),xx(n),xz(n)) | 223 | 263 | allocate(bx(n),bz(n),xx(n),xz(n)) | |
| bx=dble(B) | 224 | 264 | bx=dble(B) | |
| bz=aimag(B) | 225 | 265 | bz=aimag(B) | |
| sys=0 | 226 | 266 | sys=0 | |
| ! sys may be used to define type of solving -> see umfpack.h | 227 | 267 | ! sys may be used to define type of solving -> see umfpack.h | |
| 228 | 268 | |||
| ! Solving | 229 | 269 | ! Solving | |
| call umfpack_zi_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info) | 230 | 270 | call umfpack_zi_solve (sys, Ap, Ai, Ax,Az, xx,xz, bx,bz, numeric, control, info) | |
| ! info(1) should be zero | 231 | 271 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 232 | 272 | if (info(1) /= 0) then | |
| write(*,*) "Problem during solving" | 233 | 273 | write(*,*) "Problem during solving" | |
| status=info(1) | 234 | 274 | status=info(1) | |
| endif | 235 | 275 | endif | |
| 236 | 276 | |||
| ! free the C numeric pointer | 237 | 277 | ! free the C numeric pointer | |
| call umfpack_zi_free_numeric (numeric) | 238 | 278 | call umfpack_zi_free_numeric (numeric) | |
| 239 | 279 | |||
| x=cmplx(xx,xz,dp_kind) | 240 | 280 | x=cmplx(xx,xz,dp_kind) | |
| 241 | 281 | |||
| deallocate(bx,bz,xx,xz) | 242 | 282 | deallocate(bx,bz,xx,xz) | |
| deallocate(Ax,Az) | 243 | 283 | deallocate(Ax,Az) | |
| deallocate(Tx,Tz) | 244 | 284 | deallocate(Tx,Tz) | |
| deallocate(wTi,wTj) | 245 | 285 | deallocate(wTi,wTj) | |
| end subroutine | 246 | 286 | end subroutine | |
| 247 | 287 | |||
| 248 | 288 | |||
| 249 | 289 | |||
| 250 | 290 | |||
| 251 | 291 | |||
| 252 | 292 | |||
| subroutine fvn_dl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det) | 253 | 293 | subroutine fvn_dl_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det) | |
| implicit none | 254 | 294 | implicit none | |
| integer(kind=dp_kind), intent(in) :: n,nz | 255 | 295 | integer(kind=dp_kind), intent(in) :: n,nz | |
| real(kind=dp_kind),dimension(nz),intent(in) :: T | 256 | 296 | real(kind=dp_kind),dimension(nz),intent(in) :: T | |
| integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj | 257 | 297 | integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj | |
| real(kind=dp_kind),dimension(n),intent(in) :: B | 258 | 298 | real(kind=dp_kind),dimension(n),intent(in) :: B | |
| real(kind=dp_kind),dimension(n),intent(out) :: x | 259 | 299 | real(kind=dp_kind),dimension(n),intent(out) :: x | |
| integer(kind=dp_kind), intent(out) :: status | 260 | 300 | integer(kind=dp_kind), intent(out) :: status | |
| real(kind=dp_kind), dimension(2), optional, intent(out) :: det | 261 | 301 | real(kind=dp_kind), dimension(2), optional, intent(out) :: det | |
| 262 | 302 | |||
| integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj | 263 | 303 | integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj | |
| real(kind=dp_kind),dimension(:),allocatable :: A | 264 | 304 | real(kind=dp_kind),dimension(:),allocatable :: A | |
| integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai | 265 | 305 | integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai | |
| !integer(kind=dp_kind) :: symbolic,numeric | 266 | 306 | !integer(kind=dp_kind) :: symbolic,numeric | |
| integer(kind=dp_kind) :: symbolic,numeric | 267 | 307 | integer(kind=dp_kind) :: symbolic,numeric | |
| real(kind=dp_kind),dimension(90) :: info | 268 | 308 | real(kind=dp_kind),dimension(90) :: info | |
| real(kind=dp_kind),dimension(20) :: control | 269 | 309 | real(kind=dp_kind),dimension(20) :: control | |
| integer(kind=dp_kind) :: sys | 270 | 310 | integer(kind=dp_kind) :: sys | |
| 271 | 311 | |||
| status=0 | 272 | 312 | status=0 | |
| ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation | 273 | 313 | ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation | |
| allocate(wTi(nz),wTj(nz)) | 274 | 314 | allocate(wTi(nz),wTj(nz)) | |
| wTi=Ti-1 | 275 | 315 | wTi=Ti-1 | |
| wTj=Tj-1 | 276 | 316 | wTj=Tj-1 | |
| allocate(A(nz)) | 277 | 317 | allocate(A(nz)) | |
| allocate(Ap(n+1),Ai(nz)) | 278 | 318 | allocate(Ap(n+1),Ai(nz)) | |
| 279 | 319 | |||
| ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az | 280 | 320 | ! 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) | 281 | 321 | 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 | 282 | 322 | ! if status is not zero a problem has occured | |
| if (status /= 0) then | 283 | 323 | if (status /= 0) then | |
| write(*,*) "Problem during umfpack_dl_triplet_to_col" | 284 | 324 | write(*,*) "Problem during umfpack_dl_triplet_to_col" | |
| endif | 285 | 325 | endif | |
| 286 | 326 | |||
| ! Define defaults control values | 287 | 327 | ! Define defaults control values | |
| call umfpack_dl_defaults(control) | 288 | 328 | call umfpack_dl_defaults(control) | |
| 289 | 329 | |||
| ! Symbolic analysis | 290 | 330 | ! Symbolic analysis | |
| call umfpack_dl_symbolic(n,n,Ap,Ai,A,symbolic, control, info) | 291 | 331 | call umfpack_dl_symbolic(n,n,Ap,Ai,A,symbolic, control, info) | |
| ! info(1) should be zero | 292 | 332 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 293 | 333 | if (info(1) /= 0) then | |
| write(*,*) "Problem during symbolic analysis" | 294 | 334 | write(*,*) "Problem during symbolic analysis" | |
| status=info(1) | 295 | 335 | status=info(1) | |
| endif | 296 | 336 | endif | |
| 297 | 337 | |||
| ! Numerical factorization | 298 | 338 | ! Numerical factorization | |
| call umfpack_dl_numeric (Ap, Ai, A, symbolic, numeric, control, info) | 299 | 339 | call umfpack_dl_numeric (Ap, Ai, A, symbolic, numeric, control, info) | |
| ! info(1) should be zero | 300 | 340 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 301 | 341 | if (info(1) /= 0) then | |
| write(*,*) "Problem during numerical factorization" | 302 | 342 | write(*,*) "Problem during numerical factorization" | |
| status=info(1) | 303 | 343 | status=info(1) | |
| endif | 304 | 344 | endif | |
| 305 | 345 | |||
| ! free the C symbolic pointer | 306 | 346 | ! free the C symbolic pointer | |
| call umfpack_dl_free_symbolic (symbolic) | 307 | 347 | call umfpack_dl_free_symbolic (symbolic) | |
| 308 | 348 | |||
| ! if parameter det is present, the determinant of the matrix is calculated | 309 | 349 | ! if parameter det is present, the determinant of the matrix is calculated | |
| if (present(det) ) then | 310 | 350 | if (present(det) ) then | |
| call umfpack_dl_get_determinant(det(1),det(2),numeric,info,status) | 311 | 351 | call umfpack_dl_get_determinant(det(1),det(2),numeric,info,status) | |
| ! info(1) should be zero | 312 | 352 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 313 | 353 | if (info(1) /= 0) then | |
| if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning | 314 | 354 | if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning | |
| write(*,*) "Problem during sparse determinant, returned code : ",info(1) | 315 | 355 | write(*,*) "Problem during sparse determinant, returned code : ",info(1) | |
| endif | 316 | 356 | endif | |
| status=info(1) | 317 | 357 | status=info(1) | |
| endif | 318 | 358 | endif | |
| endif | 319 | 359 | endif | |
| 320 | 360 | |||
| 321 | 361 | |||
| sys=0 | 322 | 362 | sys=0 | |
| ! sys may be used to define type of solving -> see umfpack.h | 323 | 363 | ! sys may be used to define type of solving -> see umfpack.h | |
| 324 | 364 | |||
| ! Solving | 325 | 365 | ! Solving | |
| call umfpack_dl_solve (sys, Ap, Ai, A, x, B, numeric, control, info) | 326 | 366 | call umfpack_dl_solve (sys, Ap, Ai, A, x, B, numeric, control, info) | |
| ! info(1) should be zero | 327 | 367 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 328 | 368 | if (info(1) /= 0) then | |
| write(*,*) "Problem during solving" | 329 | 369 | write(*,*) "Problem during solving" | |
| status=info(1) | 330 | 370 | status=info(1) | |
| endif | 331 | 371 | endif | |
| 332 | 372 | |||
| ! free the C numeric pointer | 333 | 373 | ! free the C numeric pointer | |
| call umfpack_dl_free_numeric (numeric) | 334 | 374 | call umfpack_dl_free_numeric (numeric) | |
| 335 | 375 | |||
| deallocate(A) | 336 | 376 | deallocate(A) | |
| deallocate(wTi,wTj) | 337 | 377 | deallocate(wTi,wTj) | |
| end subroutine | 338 | 378 | end subroutine | |
| 339 | 379 | |||
| 340 | 380 | |||
| 341 | 381 | |||
| 342 | 382 | |||
| 343 | 383 | |||
| 344 | 384 | |||
| subroutine fvn_di_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det) | 345 | 385 | subroutine fvn_di_sparse_solve(n,nz,T,Ti,Tj,B,x,status,det) | |
| implicit none | 346 | 386 | implicit none | |
| integer(kind=sp_kind), intent(in) :: n,nz | 347 | 387 | integer(kind=sp_kind), intent(in) :: n,nz | |
| real(kind=dp_kind),dimension(nz),intent(in) :: T | 348 | 388 | real(kind=dp_kind),dimension(nz),intent(in) :: T | |
| integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj | 349 | 389 | integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj | |
| real(kind=dp_kind),dimension(n),intent(in) :: B | 350 | 390 | real(kind=dp_kind),dimension(n),intent(in) :: B | |
| real(kind=dp_kind),dimension(n),intent(out) :: x | 351 | 391 | real(kind=dp_kind),dimension(n),intent(out) :: x | |
| integer(kind=sp_kind), intent(out) :: status | 352 | 392 | integer(kind=sp_kind), intent(out) :: status | |
| real(kind=dp_kind), dimension(2), optional, intent(out) :: det | 353 | 393 | real(kind=dp_kind), dimension(2), optional, intent(out) :: det | |
| 354 | 394 | |||
| integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj | 355 | 395 | integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj | |
| real(kind=dp_kind),dimension(:),allocatable :: A | 356 | 396 | real(kind=dp_kind),dimension(:),allocatable :: A | |
| integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai | 357 | 397 | integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai | |
| !integer(kind=dp_kind) :: symbolic,numeric | 358 | 398 | !integer(kind=dp_kind) :: symbolic,numeric | |
| integer(kind=sp_kind),dimension(2) :: symbolic,numeric | 359 | 399 | integer(kind=sp_kind),dimension(2) :: symbolic,numeric | |
| ! As symbolic and numeric are used to store a C pointer, it is necessary to | 360 | 400 | ! 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 | 361 | 401 | ! still use an integer(kind=dp_kind) for 64bits machines | |
| ! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric | 362 | 402 | ! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric | |
| real(kind=dp_kind),dimension(90) :: info | 363 | 403 | real(kind=dp_kind),dimension(90) :: info | |
| real(kind=dp_kind),dimension(20) :: control | 364 | 404 | real(kind=dp_kind),dimension(20) :: control | |
| integer(kind=sp_kind) :: sys | 365 | 405 | integer(kind=sp_kind) :: sys | |
| 366 | 406 | |||
| status=0 | 367 | 407 | status=0 | |
| ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation | 368 | 408 | ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation | |
| allocate(wTi(nz),wTj(nz)) | 369 | 409 | allocate(wTi(nz),wTj(nz)) | |
| wTi=Ti-1 | 370 | 410 | wTi=Ti-1 | |
| wTj=Tj-1 | 371 | 411 | wTj=Tj-1 | |
| allocate(A(nz)) | 372 | 412 | allocate(A(nz)) | |
| allocate(Ap(n+1),Ai(nz)) | 373 | 413 | allocate(Ap(n+1),Ai(nz)) | |
| 374 | 414 | |||
| ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az | 375 | 415 | ! 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) | 376 | 416 | 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 | 377 | 417 | ! if status is not zero a problem has occured | |
| if (status /= 0) then | 378 | 418 | if (status /= 0) then | |
| write(*,*) "Problem during umfpack_di_triplet_to_col" | 379 | 419 | write(*,*) "Problem during umfpack_di_triplet_to_col" | |
| endif | 380 | 420 | endif | |
| 381 | 421 | |||
| ! Define defaults control values | 382 | 422 | ! Define defaults control values | |
| call umfpack_di_defaults(control) | 383 | 423 | call umfpack_di_defaults(control) | |
| 384 | 424 | |||
| ! Symbolic analysis | 385 | 425 | ! Symbolic analysis | |
| call umfpack_di_symbolic(n,n,Ap,Ai,A,symbolic, control, info) | 386 | 426 | call umfpack_di_symbolic(n,n,Ap,Ai,A,symbolic, control, info) | |
| ! info(1) should be zero | 387 | 427 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 388 | 428 | if (info(1) /= 0) then | |
| write(*,*) "Problem during symbolic analysis" | 389 | 429 | write(*,*) "Problem during symbolic analysis" | |
| status=info(1) | 390 | 430 | status=info(1) | |
| endif | 391 | 431 | endif | |
| 392 | 432 | |||
| ! Numerical factorization | 393 | 433 | ! Numerical factorization | |
| call umfpack_di_numeric (Ap, Ai, A, symbolic, numeric, control, info) | 394 | 434 | call umfpack_di_numeric (Ap, Ai, A, symbolic, numeric, control, info) | |
| ! info(1) should be zero | 395 | 435 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 396 | 436 | if (info(1) /= 0) then | |
| write(*,*) "Problem during numerical factorization" | 397 | 437 | write(*,*) "Problem during numerical factorization" | |
| status=info(1) | 398 | 438 | status=info(1) | |
| endif | 399 | 439 | endif | |
| 400 | 440 | |||
| ! free the C symbolic pointer | 401 | 441 | ! free the C symbolic pointer | |
| call umfpack_di_free_symbolic (symbolic) | 402 | 442 | call umfpack_di_free_symbolic (symbolic) | |
| 403 | 443 | |||
| ! if parameter det is present, the determinant of the matrix is calculated | 404 | 444 | ! if parameter det is present, the determinant of the matrix is calculated | |
| if (present(det) ) then | 405 | 445 | if (present(det) ) then | |
| call umfpack_di_get_determinant(det(1),det(2),numeric,info,status) | 406 | 446 | call umfpack_di_get_determinant(det(1),det(2),numeric,info,status) | |
| ! info(1) should be zero | 407 | 447 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 408 | 448 | if (info(1) /= 0) then | |
| if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning | 409 | 449 | if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning | |
| write(*,*) "Problem during sparse determinant, returned code : ",info(1) | 410 | 450 | write(*,*) "Problem during sparse determinant, returned code : ",info(1) | |
| endif | 411 | 451 | endif | |
| status=info(1) | 412 | 452 | status=info(1) | |
| endif | 413 | 453 | endif | |
| endif | 414 | 454 | endif | |
| 415 | 455 | |||
| sys=0 | 416 | 456 | sys=0 | |
| ! sys may be used to define type of solving -> see umfpack.h | 417 | 457 | ! sys may be used to define type of solving -> see umfpack.h | |
| ! Solving | 418 | 458 | ! Solving | |
| call umfpack_di_solve (sys, Ap, Ai, A, x, B, numeric, control, info) | 419 | 459 | call umfpack_di_solve (sys, Ap, Ai, A, x, B, numeric, control, info) | |
| ! info(1) should be zero | 420 | 460 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 421 | 461 | if (info(1) /= 0) then | |
| write(*,*) "Problem during solving" | 422 | 462 | write(*,*) "Problem during solving" | |
| status=info(1) | 423 | 463 | status=info(1) | |
| endif | 424 | 464 | endif | |
| 425 | 465 | |||
| ! free the C numeric pointer | 426 | 466 | ! free the C numeric pointer | |
| call umfpack_di_free_numeric (numeric) | 427 | 467 | call umfpack_di_free_numeric (numeric) | |
| 428 | 468 | |||
| deallocate(A) | 429 | 469 | deallocate(A) | |
| deallocate(wTi,wTj) | 430 | 470 | deallocate(wTi,wTj) | |
| end subroutine | 431 | 471 | end subroutine | |
| 432 | 472 | |||
| 433 | 473 | |||
| !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | 434 | 474 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | |
| ! | 435 | 475 | ! | |
| ! SPARSE DETERMINANT | 436 | 476 | ! SPARSE DETERMINANT | |
| ! | 437 | 477 | ! | |
| !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | 438 | 478 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | |
| subroutine fvn_zl_sparse_det(n,nz,T,Ti,Tj,det,status) | 439 | 479 | subroutine fvn_zl_sparse_det(n,nz,T,Ti,Tj,det,status) | |
| implicit none | 440 | 480 | implicit none | |
| integer(kind=dp_kind), intent(in) :: n,nz | 441 | 481 | integer(kind=dp_kind), intent(in) :: n,nz | |
| complex(kind=dp_kind),dimension(nz),intent(in) :: T | 442 | 482 | complex(kind=dp_kind),dimension(nz),intent(in) :: T | |
| integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj | 443 | 483 | integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj | |
| real(kind=dp_kind), dimension(3), intent(out) :: det | 444 | 484 | real(kind=dp_kind), dimension(3), intent(out) :: det | |
| integer(kind=dp_kind), intent(out) :: status | 445 | 485 | integer(kind=dp_kind), intent(out) :: status | |
| 446 | 486 | |||
| integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj | 447 | 487 | integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj | |
| real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz | 448 | 488 | real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz | |
| real(kind=dp_kind),dimension(:),allocatable :: Ax,Az | 449 | 489 | real(kind=dp_kind),dimension(:),allocatable :: Ax,Az | |
| integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai | 450 | 490 | integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai | |
| integer(kind=dp_kind) :: symbolic,numeric | 451 | 491 | integer(kind=dp_kind) :: symbolic,numeric | |
| real(kind=dp_kind),dimension(90) :: info | 452 | 492 | real(kind=dp_kind),dimension(90) :: info | |
| real(kind=dp_kind),dimension(20) :: control | 453 | 493 | real(kind=dp_kind),dimension(20) :: control | |
| real(kind=dp_kind) :: Mx,Mz,Ex | 454 | 494 | real(kind=dp_kind) :: Mx,Mz,Ex | |
| 455 | 495 | |||
| 456 | 496 | |||
| status=0 | 457 | 497 | status=0 | |
| 458 | 498 | |||
| ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation | 459 | 499 | ! 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 | 460 | 500 | ! Tx and Tz are the real and imaginary parts of T | |
| allocate(wTi(nz),wTj(nz)) | 461 | 501 | allocate(wTi(nz),wTj(nz)) | |
| allocate(Tx(nz),Tz(nz)) | 462 | 502 | allocate(Tx(nz),Tz(nz)) | |
| Tx=dble(T) | 463 | 503 | Tx=dble(T) | |
| Tz=aimag(T) | 464 | 504 | Tz=aimag(T) | |
| wTi=Ti-1 | 465 | 505 | wTi=Ti-1 | |
| wTj=Tj-1 | 466 | 506 | wTj=Tj-1 | |
| allocate(Ax(nz),Az(nz)) | 467 | 507 | allocate(Ax(nz),Az(nz)) | |
| allocate(Ap(n+1),Ai(nz)) | 468 | 508 | allocate(Ap(n+1),Ai(nz)) | |
| 469 | 509 | |||
| ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az | 470 | 510 | ! 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) | 471 | 511 | 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 | 472 | 512 | ! if status is not zero a problem has occured | |
| if (status /= 0) then | 473 | 513 | if (status /= 0) then | |
| write(*,*) "Problem during umfpack_zl_triplet_to_col" | 474 | 514 | write(*,*) "Problem during umfpack_zl_triplet_to_col" | |
| endif | 475 | 515 | endif | |
| 476 | 516 | |||
| ! Define defaults control values | 477 | 517 | ! Define defaults control values | |
| call umfpack_zl_defaults(control) | 478 | 518 | call umfpack_zl_defaults(control) | |
| 479 | 519 | |||
| ! Symbolic analysis | 480 | 520 | ! Symbolic analysis | |
| call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) | 481 | 521 | call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) | |
| ! info(1) should be zero | 482 | 522 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 483 | 523 | if (info(1) /= 0) then | |
| write(*,*) "Problem during symbolic analysis" | 484 | 524 | write(*,*) "Problem during symbolic analysis" | |
| status=info(1) | 485 | 525 | status=info(1) | |
| endif | 486 | 526 | endif | |
| 487 | 527 | |||
| ! Numerical factorization | 488 | 528 | ! Numerical factorization | |
| call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) | 489 | 529 | call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) | |
| ! info(1) should be zero | 490 | 530 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 491 | 531 | if (info(1) /= 0) then | |
| write(*,*) "Problem during numerical factorization" | 492 | 532 | write(*,*) "Problem during numerical factorization" | |
| status=info(1) | 493 | 533 | status=info(1) | |
| endif | 494 | 534 | endif | |
| 495 | 535 | |||
| ! free the C symbolic pointer | 496 | 536 | ! free the C symbolic pointer | |
| call umfpack_zl_free_symbolic (symbolic) | 497 | 537 | call umfpack_zl_free_symbolic (symbolic) | |
| 498 | 538 | |||
| call umfpack_zl_get_determinant(det(1),det(2),det(3),numeric,info,status) | 499 | 539 | call umfpack_zl_get_determinant(det(1),det(2),det(3),numeric,info,status) | |
| ! info(1) should be zero | 500 | 540 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 501 | 541 | if (info(1) /= 0) then | |
| if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning | 502 | 542 | if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning | |
| write(*,*) "Problem during sparse determinant, returned code : ",info(1) | 503 | 543 | write(*,*) "Problem during sparse determinant, returned code : ",info(1) | |
| endif | 504 | 544 | endif | |
| status=info(1) | 505 | 545 | status=info(1) | |
| endif | 506 | 546 | endif | |
| 507 | 547 | |||
| ! free the C numeric pointer | 508 | 548 | ! free the C numeric pointer | |
| call umfpack_zl_free_numeric (numeric) | 509 | 549 | call umfpack_zl_free_numeric (numeric) | |
| 510 | 550 | |||
| deallocate(Ax,Az) | 511 | 551 | deallocate(Ax,Az) | |
| deallocate(Tx,Tz) | 512 | 552 | deallocate(Tx,Tz) | |
| deallocate(wTi,wTj) | 513 | 553 | deallocate(wTi,wTj) | |
| end subroutine | 514 | 554 | end subroutine | |
| 515 | 555 | |||
| 516 | 556 | |||
| subroutine fvn_zi_sparse_det(n,nz,T,Ti,Tj,det,status) | 517 | 557 | subroutine fvn_zi_sparse_det(n,nz,T,Ti,Tj,det,status) | |
| implicit none | 518 | 558 | implicit none | |
| integer(kind=sp_kind), intent(in) :: n,nz | 519 | 559 | integer(kind=sp_kind), intent(in) :: n,nz | |
| complex(kind=dp_kind),dimension(nz),intent(in) :: T | 520 | 560 | complex(kind=dp_kind),dimension(nz),intent(in) :: T | |
| integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj | 521 | 561 | integer(kind=sp_kind),dimension(nz),intent(in) :: Ti,Tj | |
| integer(kind=sp_kind), intent(out) :: status | 522 | 562 | integer(kind=sp_kind), intent(out) :: status | |
| real(kind=dp_kind), dimension(3), intent(out) :: det | 523 | 563 | real(kind=dp_kind), dimension(3), intent(out) :: det | |
| 524 | 564 | |||
| integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj | 525 | 565 | integer(kind=sp_kind),dimension(:),allocatable :: wTi,wTj | |
| real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz | 526 | 566 | real(kind=dp_kind),dimension(:),allocatable :: Tx,Tz | |
| real(kind=dp_kind),dimension(:),allocatable :: Ax,Az | 527 | 567 | real(kind=dp_kind),dimension(:),allocatable :: Ax,Az | |
| integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai | 528 | 568 | integer(kind=sp_kind),dimension(:),allocatable :: Ap,Ai | |
| !integer(kind=dp_kind) :: symbolic,numeric | 529 | 569 | !integer(kind=dp_kind) :: symbolic,numeric | |
| integer(kind=sp_kind),dimension(2) :: symbolic,numeric | 530 | 570 | integer(kind=sp_kind),dimension(2) :: symbolic,numeric | |
| ! As symbolic and numeric are used to store a C pointer, it is necessary to | 531 | 571 | ! 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 | 532 | 572 | ! still use an integer(kind=dp_kind) for 64bits machines | |
| ! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric | 533 | 573 | ! An other possibility : integer(kind=sp_kind),dimension(2) :: symbolic,numeric | |
| real(kind=dp_kind),dimension(90) :: info | 534 | 574 | real(kind=dp_kind),dimension(90) :: info | |
| real(kind=dp_kind),dimension(20) :: control | 535 | 575 | real(kind=dp_kind),dimension(20) :: control | |
| 536 | 576 | |||
| status=0 | 537 | 577 | status=0 | |
| ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation | 538 | 578 | ! 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 | 539 | 579 | ! Tx and Tz are the real and imaginary parts of T | |
| allocate(wTi(nz),wTj(nz)) | 540 | 580 | allocate(wTi(nz),wTj(nz)) | |
| allocate(Tx(nz),Tz(nz)) | 541 | 581 | allocate(Tx(nz),Tz(nz)) | |
| Tx=dble(T) | 542 | 582 | Tx=dble(T) | |
| Tz=aimag(T) | 543 | 583 | Tz=aimag(T) | |
| wTi=Ti-1 | 544 | 584 | wTi=Ti-1 | |
| wTj=Tj-1 | 545 | 585 | wTj=Tj-1 | |
| allocate(Ax(nz),Az(nz)) | 546 | 586 | allocate(Ax(nz),Az(nz)) | |
| allocate(Ap(n+1),Ai(nz)) | 547 | 587 | allocate(Ap(n+1),Ai(nz)) | |
| 548 | 588 | |||
| ! perform the triplet to compressed column form -> Ap,Ai,Ax,Az | 549 | 589 | ! 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) | 550 | 590 | 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 | 551 | 591 | ! if status is not zero a problem has occured | |
| if (status /= 0) then | 552 | 592 | if (status /= 0) then | |
| write(*,*) "Problem during umfpack_zl_triplet_to_col" | 553 | 593 | write(*,*) "Problem during umfpack_zl_triplet_to_col" | |
| endif | 554 | 594 | endif | |
| 555 | 595 | |||
| ! Define defaults control values | 556 | 596 | ! Define defaults control values | |
| call umfpack_zi_defaults(control) | 557 | 597 | call umfpack_zi_defaults(control) | |
| 558 | 598 | |||
| ! Symbolic analysis | 559 | 599 | ! Symbolic analysis | |
| call umfpack_zi_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) | 560 | 600 | call umfpack_zi_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info) | |
| ! info(1) should be zero | 561 | 601 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 562 | 602 | if (info(1) /= 0) then | |
| write(*,*) "Problem during symbolic analysis" | 563 | 603 | write(*,*) "Problem during symbolic analysis" | |
| status=info(1) | 564 | 604 | status=info(1) | |
| endif | 565 | 605 | endif | |
| 566 | 606 | |||
| ! Numerical factorization | 567 | 607 | ! Numerical factorization | |
| call umfpack_zi_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) | 568 | 608 | call umfpack_zi_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info) | |
| ! info(1) should be zero | 569 | 609 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 570 | 610 | if (info(1) /= 0) then | |
| write(*,*) "Problem during numerical factorization" | 571 | 611 | write(*,*) "Problem during numerical factorization" | |
| status=info(1) | 572 | 612 | status=info(1) | |
| endif | 573 | 613 | endif | |
| 574 | 614 | |||
| ! free the C symbolic pointer | 575 | 615 | ! free the C symbolic pointer | |
| call umfpack_zi_free_symbolic (symbolic) | 576 | 616 | call umfpack_zi_free_symbolic (symbolic) | |
| 577 | 617 | |||
| call umfpack_zi_get_determinant(det(1),det(2),det(3),numeric,info,status) | 578 | 618 | call umfpack_zi_get_determinant(det(1),det(2),det(3),numeric,info,status) | |
| ! info(1) should be zero | 579 | 619 | ! info(1) should be zero | |
| if (info(1) /= 0) then | 580 | 620 | if (info(1) /= 0) then | |
| if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning | 581 | 621 | if ( (info(1) < 1) .or. (info(1) >3) ) then ! not a warning | |
| write(*,*) "Problem during sparse determinant, returned code : ",info(1) | 582 | 622 | write(*,*) "Problem during sparse determinant, returned code : ",info(1) | |
| endif | 583 | 623 | endif | |
| status=info(1) | 584 | 624 | status=info(1) | |
| endif | 585 | 625 | endif | |
| 586 | 626 | |||
| ! free the C numeric pointer | 587 | 627 | ! free the C numeric pointer | |
| call umfpack_zi_free_numeric (numeric) | 588 | 628 | call umfpack_zi_free_numeric (numeric) | |
| 589 | 629 | |||
| deallocate(Ax,Az) | 590 | 630 | deallocate(Ax,Az) | |
| deallocate(Tx,Tz) | 591 | 631 | deallocate(Tx,Tz) | |
| deallocate(wTi,wTj) | 592 | 632 | deallocate(wTi,wTj) | |
| end subroutine | 593 | 633 | end subroutine | |
| 594 | 634 | |||
| subroutine fvn_dl_sparse_det(n,nz,T,Ti,Tj,det,status) | 595 | 635 | subroutine fvn_dl_sparse_det(n,nz,T,Ti,Tj,det,status) | |
| implicit none | 596 | 636 | implicit none | |
| integer(kind=dp_kind), intent(in) :: n,nz | 597 | 637 | integer(kind=dp_kind), intent(in) :: n,nz | |
| real(kind=dp_kind),dimension(nz),intent(in) :: T | 598 | 638 | real(kind=dp_kind),dimension(nz),intent(in) :: T | |
| integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj | 599 | 639 | integer(kind=dp_kind),dimension(nz),intent(in) :: Ti,Tj | |
| integer(kind=dp_kind), intent(out) :: status | 600 | 640 | integer(kind=dp_kind), intent(out) :: status | |
| real(kind=dp_kind), dimension(2), intent(out) :: det | 601 | 641 | real(kind=dp_kind), dimension(2), intent(out) :: det | |
| 602 | 642 | |||
| integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj | 603 | 643 | integer(kind=dp_kind),dimension(:),allocatable :: wTi,wTj | |
| real(kind=dp_kind),dimension(:),allocatable :: A | 604 | 644 | real(kind=dp_kind),dimension(:),allocatable :: A | |
| integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai | 605 | 645 | integer(kind=dp_kind),dimension(:),allocatable :: Ap,Ai | |
| !integer(kind=dp_kind) :: symbolic,numeric | 606 | 646 | !integer(kind=dp_kind) :: symbolic,numeric | |
| integer(kind=dp_kind) :: symbolic,numeric | 607 | 647 | integer(kind=dp_kind) :: symbolic,numeric | |
| real(kind=dp_kind),dimension(90) :: info | 608 | 648 | real(kind=dp_kind),dimension(90) :: info | |
| real(kind=dp_kind),dimension(20) :: control | 609 | 649 | real(kind=dp_kind),dimension(20) :: control | |
| 610 | 650 | |||
| status=0 | 611 | 651 | status=0 |