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