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