lxsd / fvn

Download as
Browse Code ยป

Commit 42627b85e3786e526a4c2c75e81ec4bd80890a42

Authored by daniau
1 parent fdbb9bb0a5
Exists in master and in 3 other branches geevx, geevx_2020, multi

git-svn-id: https://lxsd.femto-st.fr/svn/fvn@48 b657c933-2333-4658-acf2-d3c7c2708721

Showing 1 changed file with 16 additions and 8 deletions Inline Diff

fvn_linear/fvn_linear.f90
Diff comments   View file @ 42627b8
module fvn_linear 1 1 module fvn_linear
use fvn_common 2 2 use fvn_common
implicit none 3 3 implicit none
4 4
5 5
! 6 6 !
! Interfaces for matrix operators 7 7 ! Interfaces for matrix operators
8 8
! .x. matrix multiplication 9 9 ! .x. matrix multiplication
interface operator(.x.) 10 10 interface operator(.x.)
module procedure fvn_op_s_matmul,fvn_op_d_matmul,fvn_op_c_matmul,fvn_op_z_matmul 11 11 module procedure fvn_op_s_matmul,fvn_op_d_matmul,fvn_op_c_matmul,fvn_op_z_matmul
end interface 12 12 end interface
13 13
! .t. matrix transposition 14 14 ! .t. matrix transposition
interface operator(.t.) 15 15 interface operator(.t.)
module procedure fvn_op_s_transpose,fvn_op_d_transpose,fvn_op_c_transpose,fvn_op_z_transpose 16 16 module procedure fvn_op_s_transpose,fvn_op_d_transpose,fvn_op_c_transpose,fvn_op_z_transpose
end interface 17 17 end interface
18 18
! .tx. transpose first operand and multiply 19 19 ! .tx. transpose first operand and multiply
interface operator(.tx.) 20 20 interface operator(.tx.)
module procedure fvn_op_s_tx,fvn_op_d_tx,fvn_op_c_tx,fvn_op_z_tx 21 21 module procedure fvn_op_s_tx,fvn_op_d_tx,fvn_op_c_tx,fvn_op_z_tx
end interface 22 22 end interface
23 23
! .xt. transpose second operand and multiply 24 24 ! .xt. transpose second operand and multiply
interface operator(.xt.) 25 25 interface operator(.xt.)
module procedure fvn_op_s_xt,fvn_op_d_xt,fvn_op_c_xt,fvn_op_z_xt 26 26 module procedure fvn_op_s_xt,fvn_op_d_xt,fvn_op_c_xt,fvn_op_z_xt
end interface 27 27 end interface
28 28
! .i. inverse matrix 29 29 ! .i. inverse matrix
interface operator(.i.) 30 30 interface operator(.i.)
module procedure fvn_op_s_matinv,fvn_op_d_matinv,fvn_op_c_matinv,fvn_op_z_matinv 31 31 module procedure fvn_op_s_matinv,fvn_op_d_matinv,fvn_op_c_matinv,fvn_op_z_matinv
end interface 32 32 end interface
33 33
! .ix. inverse first operand and multiply 34 34 ! .ix. inverse first operand and multiply
interface operator(.ix.) 35 35 interface operator(.ix.)
module procedure fvn_op_s_ix,fvn_op_d_ix,fvn_op_c_ix,fvn_op_z_ix 36 36 module procedure fvn_op_s_ix,fvn_op_d_ix,fvn_op_c_ix,fvn_op_z_ix
end interface 37 37 end interface
38 38
! .xi. inverse second operand and multiply 39 39 ! .xi. inverse second operand and multiply
interface operator(.xi.) 40 40 interface operator(.xi.)
module procedure fvn_op_s_xi,fvn_op_d_xi,fvn_op_c_xi,fvn_op_z_xi 41 41 module procedure fvn_op_s_xi,fvn_op_d_xi,fvn_op_c_xi,fvn_op_z_xi
end interface 42 42 end interface
43 43
! .h. transpose conjugate (adjoint) 44 44 ! .h. transpose conjugate (adjoint)
interface operator(.h.) 45 45 interface operator(.h.)
module procedure fvn_op_s_transpose,fvn_op_d_transpose,fvn_op_c_conj_transpose,fvn_op_z_conj_transpose 46 46 module procedure fvn_op_s_transpose,fvn_op_d_transpose,fvn_op_c_conj_transpose,fvn_op_z_conj_transpose
end interface 47 47 end interface
48 48
! .hx. transpose conjugate first operand and multiply 49 49 ! .hx. transpose conjugate first operand and multiply
interface operator(.hx.) 50 50 interface operator(.hx.)
module procedure fvn_op_s_tx,fvn_op_d_tx,fvn_op_c_hx,fvn_op_z_hx 51 51 module procedure fvn_op_s_tx,fvn_op_d_tx,fvn_op_c_hx,fvn_op_z_hx
end interface 52 52 end interface
53 53
! .xh. transpose conjugate second operand and multiply 54 54 ! .xh. transpose conjugate second operand and multiply
interface operator(.xh.) 55 55 interface operator(.xh.)
module procedure fvn_op_s_xt,fvn_op_d_xt,fvn_op_c_xh,fvn_op_z_xh 56 56 module procedure fvn_op_s_xt,fvn_op_d_xt,fvn_op_c_xh,fvn_op_z_xh
end interface 57 57 end interface
58 58
59 59
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 60 60 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 61 61 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Generic interface Definition 62 62 ! Generic interface Definition
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 63 63 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 64 64 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
65 65
! Matrix inversion 66 66 ! Matrix inversion
interface fvn_matinv 67 67 interface fvn_matinv
module procedure fvn_s_matinv,fvn_d_matinv,fvn_c_matinv,fvn_z_matinv 68 68 module procedure fvn_s_matinv,fvn_d_matinv,fvn_c_matinv,fvn_z_matinv
end interface fvn_matinv 69 69 end interface fvn_matinv
70 70
! Determinant 71 71 ! Determinant
interface fvn_det 72 72 interface fvn_det
module procedure fvn_s_det,fvn_d_det,fvn_c_det,fvn_z_det 73 73 module procedure fvn_s_det,fvn_d_det,fvn_c_det,fvn_z_det
end interface fvn_det 74 74 end interface fvn_det
75 75
! Condition 76 76 ! Condition
interface fvn_matcon 77 77 interface fvn_matcon
module procedure fvn_s_matcon,fvn_d_matcon,fvn_c_matcon,fvn_z_matcon 78 78 module procedure fvn_s_matcon,fvn_d_matcon,fvn_c_matcon,fvn_z_matcon
end interface fvn_matcon 79 79 end interface fvn_matcon
80 80
! Eigen 81 81 ! Eigen
interface fvn_matev 82 82 interface fvn_matev
module procedure fvn_s_matev,fvn_d_matev,fvn_c_matev,fvn_z_matev 83 83 module procedure fvn_s_matev,fvn_d_matev,fvn_c_matev,fvn_z_matev
end interface fvn_matev 84 84 end interface fvn_matev
85 85
! Least square polynomial 86 86 ! Least square polynomial
interface fvn_lspoly 87 87 interface fvn_lspoly
module procedure fvn_s_lspoly,fvn_d_lspoly 88 88 module procedure fvn_s_lspoly,fvn_d_lspoly
end interface fvn_lspoly 89 89 end interface fvn_lspoly
90 90
91 91
contains 92 92 contains
93 93
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 94 94 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 95 95 !
! Linear operators 96 96 ! Linear operators
! 97 97 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 98 98 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
99 99
! 100 100 !
! .x. 101 101 ! .x.
! 102 102 !
function fvn_op_s_matmul(a,b) 103 103 function fvn_op_s_matmul(a,b)
implicit none 104 104 implicit none
real(4), dimension(:,:),intent(in) :: a,b 105 105 real(4), dimension(:,:),intent(in) :: a,b
real(4), dimension(size(a,1),size(b,2)) :: fvn_op_s_matmul 106 106 real(4), dimension(size(a,1),size(b,2)) :: fvn_op_s_matmul
fvn_op_s_matmul=matmul(a,b) 107 107 fvn_op_s_matmul=matmul(a,b)
end function 108 108 end function
function fvn_op_d_matmul(a,b) 109 109 function fvn_op_d_matmul(a,b)
implicit none 110 110 implicit none
real(8), dimension(:,:),intent(in) :: a,b 111 111 real(8), dimension(:,:),intent(in) :: a,b
real(8), dimension(size(a,1),size(b,2)) :: fvn_op_d_matmul 112 112 real(8), dimension(size(a,1),size(b,2)) :: fvn_op_d_matmul
fvn_op_d_matmul=matmul(a,b) 113 113 fvn_op_d_matmul=matmul(a,b)
end function 114 114 end function
function fvn_op_c_matmul(a,b) 115 115 function fvn_op_c_matmul(a,b)
implicit none 116 116 implicit none
complex(4), dimension(:,:),intent(in) :: a,b 117 117 complex(4), dimension(:,:),intent(in) :: a,b
complex(4), dimension(size(a,1),size(b,2)) :: fvn_op_c_matmul 118 118 complex(4), dimension(size(a,1),size(b,2)) :: fvn_op_c_matmul
fvn_op_c_matmul=matmul(a,b) 119 119 fvn_op_c_matmul=matmul(a,b)
end function 120 120 end function
function fvn_op_z_matmul(a,b) 121 121 function fvn_op_z_matmul(a,b)
implicit none 122 122 implicit none
complex(8), dimension(:,:), intent(in) :: a,b 123 123 complex(8), dimension(:,:), intent(in) :: a,b
complex(8), dimension(size(a,1),size(b,2)) :: fvn_op_z_matmul 124 124 complex(8), dimension(size(a,1),size(b,2)) :: fvn_op_z_matmul
fvn_op_z_matmul=matmul(a,b) 125 125 fvn_op_z_matmul=matmul(a,b)
end function 126 126 end function
127 127
! 128 128 !
! .tx. 129 129 ! .tx.
! 130 130 !
function fvn_op_s_tx(a,b) 131 131 function fvn_op_s_tx(a,b)
implicit none 132 132 implicit none
real(4), dimension(:,:), intent(in) :: a,b 133 133 real(4), dimension(:,:), intent(in) :: a,b
real(4), dimension(size(a,2),size(b,2)) :: fvn_op_s_tx 134 134 real(4), dimension(size(a,2),size(b,2)) :: fvn_op_s_tx
fvn_op_s_tx=matmul(transpose(a),b) 135 135 fvn_op_s_tx=matmul(transpose(a),b)
end function 136 136 end function
function fvn_op_d_tx(a,b) 137 137 function fvn_op_d_tx(a,b)
implicit none 138 138 implicit none
real(8), dimension(:,:), intent(in) :: a,b 139 139 real(8), dimension(:,:), intent(in) :: a,b
real(8), dimension(size(a,2),size(b,2)) :: fvn_op_d_tx 140 140 real(8), dimension(size(a,2),size(b,2)) :: fvn_op_d_tx
fvn_op_d_tx=matmul(transpose(a),b) 141 141 fvn_op_d_tx=matmul(transpose(a),b)
end function 142 142 end function
function fvn_op_c_tx(a,b) 143 143 function fvn_op_c_tx(a,b)
implicit none 144 144 implicit none
complex(4), dimension(:,:), intent(in) :: a,b 145 145 complex(4), dimension(:,:), intent(in) :: a,b
complex(4), dimension(size(a,2),size(b,2)) :: fvn_op_c_tx 146 146 complex(4), dimension(size(a,2),size(b,2)) :: fvn_op_c_tx
fvn_op_c_tx=matmul(transpose(a),b) 147 147 fvn_op_c_tx=matmul(transpose(a),b)
end function 148 148 end function
function fvn_op_z_tx(a,b) 149 149 function fvn_op_z_tx(a,b)
implicit none 150 150 implicit none
complex(8), dimension(:,:), intent(in) :: a,b 151 151 complex(8), dimension(:,:), intent(in) :: a,b
complex(8), dimension(size(a,2),size(b,2)) :: fvn_op_z_tx 152 152 complex(8), dimension(size(a,2),size(b,2)) :: fvn_op_z_tx
fvn_op_z_tx=matmul(transpose(a),b) 153 153 fvn_op_z_tx=matmul(transpose(a),b)
end function 154 154 end function
155 155
156 156
! 157 157 !
! .xt. 158 158 ! .xt.
! 159 159 !
function fvn_op_s_xt(a,b) 160 160 function fvn_op_s_xt(a,b)
implicit none 161 161 implicit none
real(4), dimension(:,:), intent(in) :: a,b 162 162 real(4), dimension(:,:), intent(in) :: a,b
real(4), dimension(size(a,1),size(b,1)) :: fvn_op_s_xt 163 163 real(4), dimension(size(a,1),size(b,1)) :: fvn_op_s_xt
fvn_op_s_xt=matmul(a,transpose(b)) 164 164 fvn_op_s_xt=matmul(a,transpose(b))
end function 165 165 end function
function fvn_op_d_xt(a,b) 166 166 function fvn_op_d_xt(a,b)
implicit none 167 167 implicit none
real(8), dimension(:,:), intent(in) :: a,b 168 168 real(8), dimension(:,:), intent(in) :: a,b
real(8), dimension(size(a,1),size(b,1)) :: fvn_op_d_xt 169 169 real(8), dimension(size(a,1),size(b,1)) :: fvn_op_d_xt
fvn_op_d_xt=matmul(a,transpose(b)) 170 170 fvn_op_d_xt=matmul(a,transpose(b))
end function 171 171 end function
function fvn_op_c_xt(a,b) 172 172 function fvn_op_c_xt(a,b)
implicit none 173 173 implicit none
complex(4), dimension(:,:), intent(in) :: a,b 174 174 complex(4), dimension(:,:), intent(in) :: a,b
complex(4), dimension(size(a,1),size(b,1)) :: fvn_op_c_xt 175 175 complex(4), dimension(size(a,1),size(b,1)) :: fvn_op_c_xt
fvn_op_c_xt=matmul(a,transpose(b)) 176 176 fvn_op_c_xt=matmul(a,transpose(b))
end function 177 177 end function
function fvn_op_z_xt(a,b) 178 178 function fvn_op_z_xt(a,b)
implicit none 179 179 implicit none
complex(8), dimension(:,:), intent(in) :: a,b 180 180 complex(8), dimension(:,:), intent(in) :: a,b
complex(8), dimension(size(a,1),size(b,1)) :: fvn_op_z_xt 181 181 complex(8), dimension(size(a,1),size(b,1)) :: fvn_op_z_xt
fvn_op_z_xt=matmul(a,transpose(b)) 182 182 fvn_op_z_xt=matmul(a,transpose(b))
end function 183 183 end function
184 184
! 185 185 !
! .t. 186 186 ! .t.
! 187 187 !
function fvn_op_s_transpose(a) 188 188 function fvn_op_s_transpose(a)
implicit none 189 189 implicit none
real(4),dimension(:,:),intent(in) :: a 190 190 real(4),dimension(:,:),intent(in) :: a
real(4),dimension(size(a,2),size(a,1)) :: fvn_op_s_transpose 191 191 real(4),dimension(size(a,2),size(a,1)) :: fvn_op_s_transpose
fvn_op_s_transpose=transpose(a) 192 192 fvn_op_s_transpose=transpose(a)
end function 193 193 end function
function fvn_op_d_transpose(a) 194 194 function fvn_op_d_transpose(a)
implicit none 195 195 implicit none
real(8),dimension(:,:),intent(in) :: a 196 196 real(8),dimension(:,:),intent(in) :: a
real(8),dimension(size(a,2),size(a,1)) :: fvn_op_d_transpose 197 197 real(8),dimension(size(a,2),size(a,1)) :: fvn_op_d_transpose
fvn_op_d_transpose=transpose(a) 198 198 fvn_op_d_transpose=transpose(a)
end function 199 199 end function
function fvn_op_c_transpose(a) 200 200 function fvn_op_c_transpose(a)
implicit none 201 201 implicit none
complex(4),dimension(:,:),intent(in) :: a 202 202 complex(4),dimension(:,:),intent(in) :: a
complex(4),dimension(size(a,2),size(a,1)) :: fvn_op_c_transpose 203 203 complex(4),dimension(size(a,2),size(a,1)) :: fvn_op_c_transpose
fvn_op_c_transpose=transpose(a) 204 204 fvn_op_c_transpose=transpose(a)
end function 205 205 end function
function fvn_op_z_transpose(a) 206 206 function fvn_op_z_transpose(a)
implicit none 207 207 implicit none
complex(8),dimension(:,:),intent(in) :: a 208 208 complex(8),dimension(:,:),intent(in) :: a
complex(8),dimension(size(a,2),size(a,1)) :: fvn_op_z_transpose 209 209 complex(8),dimension(size(a,2),size(a,1)) :: fvn_op_z_transpose
fvn_op_z_transpose=transpose(a) 210 210 fvn_op_z_transpose=transpose(a)
end function 211 211 end function
212 212
! 213 213 !
! .i. 214 214 ! .i.
! 215 215 !
216 ! It seems that there's a problem with automatic arrays with gfortran
217 ! in some circumstances. To allow compilation with gfortran we use here a temporary array
218 ! for the call. Without that there's a warning at compile time and a segmentation fault
219 ! during execution. This is odd as we double memory use.
function fvn_op_s_matinv(a) 216 220 function fvn_op_s_matinv(a)
implicit none 217 221 implicit none
real(4),dimension(:,:),intent(in) :: a 218 222 real(4),dimension(:,:),intent(in) :: a
real(4),dimension(size(a,1),size(a,1)) :: fvn_op_s_matinv 219 223 real(4),dimension(size(a,1),size(a,1)) :: fvn_op_s_matinv,tmp_array
call fvn_s_matinv(size(a,1),a,fvn_op_s_matinv,fvn_status) 220 224 call fvn_s_matinv(size(a,1),a,tmp_array,fvn_status)
225 fvn_op_s_matinv=tmp_array
end function 221 226 end function
function fvn_op_d_matinv(a) 222 227 function fvn_op_d_matinv(a)
implicit none 223 228 implicit none
real(8),dimension(:,:),intent(in) :: a 224 229 real(8),dimension(:,:),intent(in) :: a
real(8),dimension(size(a,1),size(a,1)) :: fvn_op_d_matinv 225 230 real(8),dimension(size(a,1),size(a,1)) :: fvn_op_d_matinv,tmp_array
call fvn_d_matinv(size(a,1),a,fvn_op_d_matinv,fvn_status) 226 231 call fvn_d_matinv(size(a,1),a,tmp_array,fvn_status)
232 fvn_op_d_matinv=tmp_array
end function 227 233 end function
function fvn_op_c_matinv(a) 228 234 function fvn_op_c_matinv(a)
implicit none 229 235 implicit none
complex(4),dimension(:,:),intent(in) :: a 230 236 complex(4),dimension(:,:),intent(in) :: a
complex(4),dimension(size(a,1),size(a,1)) :: fvn_op_c_matinv 231 237 complex(4),dimension(size(a,1),size(a,1)) :: fvn_op_c_matinv,tmp_array
call fvn_c_matinv(size(a,1),a,fvn_op_c_matinv,fvn_status) 232 238 call fvn_c_matinv(size(a,1),a,tmp_array,fvn_status)
239 fvn_op_c_matinv=tmp_array
end function 233 240 end function
function fvn_op_z_matinv(a) 234 241 function fvn_op_z_matinv(a)
implicit none 235 242 implicit none
complex(8),dimension(:,:),intent(in) :: a 236 243 complex(8),dimension(:,:),intent(in) :: a
complex(8),dimension(size(a,1),size(a,1)) :: fvn_op_z_matinv 237 244 complex(8),dimension(size(a,1),size(a,1)) :: fvn_op_z_matinv,tmp_array
call fvn_z_matinv(size(a,1),a,fvn_op_z_matinv,fvn_status) 238 245 call fvn_z_matinv(size(a,1),a,tmp_array,fvn_status)
246 fvn_op_z_matinv=tmp_array
end function 239 247 end function
240 248
! 241 249 !
! .ix. 242 250 ! .ix.
! 243 251 !
function fvn_op_s_ix(a,b) 244 252 function fvn_op_s_ix(a,b)
implicit none 245 253 implicit none
real(4), dimension(:,:), intent(in) :: a,b 246 254 real(4), dimension(:,:), intent(in) :: a,b
real(4), dimension(size(a,1),size(b,2)) :: fvn_op_s_ix 247 255 real(4), dimension(size(a,1),size(b,2)) :: fvn_op_s_ix
fvn_op_s_ix=matmul(fvn_op_s_matinv(a),b) 248 256 fvn_op_s_ix=matmul(fvn_op_s_matinv(a),b)
end function 249 257 end function
function fvn_op_d_ix(a,b) 250 258 function fvn_op_d_ix(a,b)
implicit none 251 259 implicit none
real(8), dimension(:,:), intent(in) :: a,b 252 260 real(8), dimension(:,:), intent(in) :: a,b
real(8), dimension(size(a,1),size(b,2)) :: fvn_op_d_ix 253 261 real(8), dimension(size(a,1),size(b,2)) :: fvn_op_d_ix
fvn_op_d_ix=matmul(fvn_op_d_matinv(a),b) 254 262 fvn_op_d_ix=matmul(fvn_op_d_matinv(a),b)
end function 255 263 end function
function fvn_op_c_ix(a,b) 256 264 function fvn_op_c_ix(a,b)
implicit none 257 265 implicit none
complex(4), dimension(:,:), intent(in) :: a,b 258 266 complex(4), dimension(:,:), intent(in) :: a,b
complex(4), dimension(size(a,1),size(b,2)) :: fvn_op_c_ix 259 267 complex(4), dimension(size(a,1),size(b,2)) :: fvn_op_c_ix
fvn_op_c_ix=matmul(fvn_op_c_matinv(a),b) 260 268 fvn_op_c_ix=matmul(fvn_op_c_matinv(a),b)
end function 261 269 end function
function fvn_op_z_ix(a,b) 262 270 function fvn_op_z_ix(a,b)
implicit none 263 271 implicit none
complex(8), dimension(:,:), intent(in) :: a,b 264 272 complex(8), dimension(:,:), intent(in) :: a,b
complex(8), dimension(size(a,1),size(b,2)) :: fvn_op_z_ix 265 273 complex(8), dimension(size(a,1),size(b,2)) :: fvn_op_z_ix
fvn_op_z_ix=matmul(fvn_op_z_matinv(a),b) 266 274 fvn_op_z_ix=matmul(fvn_op_z_matinv(a),b)
end function 267 275 end function
268 276
! 269 277 !
! .xi. 270 278 ! .xi.
! 271 279 !
function fvn_op_s_xi(a,b) 272 280 function fvn_op_s_xi(a,b)
implicit none 273 281 implicit none
real(4), dimension(:,:), intent(in) :: a,b 274 282 real(4), dimension(:,:), intent(in) :: a,b
real(4), dimension(size(a,1),size(b,2)) :: fvn_op_s_xi 275 283 real(4), dimension(size(a,1),size(b,2)) :: fvn_op_s_xi
fvn_op_s_xi=matmul(a,fvn_op_s_matinv(b)) 276 284 fvn_op_s_xi=matmul(a,fvn_op_s_matinv(b))
end function 277 285 end function
function fvn_op_d_xi(a,b) 278 286 function fvn_op_d_xi(a,b)
implicit none 279 287 implicit none
real(8), dimension(:,:), intent(in) :: a,b 280 288 real(8), dimension(:,:), intent(in) :: a,b
real(8), dimension(size(a,1),size(b,2)) :: fvn_op_d_xi 281 289 real(8), dimension(size(a,1),size(b,2)) :: fvn_op_d_xi
fvn_op_d_xi=matmul(a,fvn_op_d_matinv(b)) 282 290 fvn_op_d_xi=matmul(a,fvn_op_d_matinv(b))
end function 283 291 end function
function fvn_op_c_xi(a,b) 284 292 function fvn_op_c_xi(a,b)
implicit none 285 293 implicit none
complex(4), dimension(:,:), intent(in) :: a,b 286 294 complex(4), dimension(:,:), intent(in) :: a,b
complex(4), dimension(size(a,1),size(b,2)) :: fvn_op_c_xi 287 295 complex(4), dimension(size(a,1),size(b,2)) :: fvn_op_c_xi
fvn_op_c_xi=matmul(a,fvn_op_c_matinv(b)) 288 296 fvn_op_c_xi=matmul(a,fvn_op_c_matinv(b))
end function 289 297 end function
function fvn_op_z_xi(a,b) 290 298 function fvn_op_z_xi(a,b)
implicit none 291 299 implicit none
complex(8), dimension(:,:), intent(in) :: a,b 292 300 complex(8), dimension(:,:), intent(in) :: a,b
complex(8), dimension(size(a,1),size(b,2)) :: fvn_op_z_xi 293 301 complex(8), dimension(size(a,1),size(b,2)) :: fvn_op_z_xi
fvn_op_z_xi=matmul(a,fvn_op_z_matinv(b)) 294 302 fvn_op_z_xi=matmul(a,fvn_op_z_matinv(b))
end function 295 303 end function
296 304
! 297 305 !
! .h. 298 306 ! .h.
! 299 307 !
function fvn_op_c_conj_transpose(a) 300 308 function fvn_op_c_conj_transpose(a)
implicit none 301 309 implicit none
complex(4),dimension(:,:),intent(in) :: a 302 310 complex(4),dimension(:,:),intent(in) :: a
complex(4),dimension(size(a,2),size(a,1)) :: fvn_op_c_conj_transpose 303 311 complex(4),dimension(size(a,2),size(a,1)) :: fvn_op_c_conj_transpose
fvn_op_c_conj_transpose=transpose(conjg(a)) 304 312 fvn_op_c_conj_transpose=transpose(conjg(a))
end function 305 313 end function
function fvn_op_z_conj_transpose(a) 306 314 function fvn_op_z_conj_transpose(a)
implicit none 307 315 implicit none
complex(8),dimension(:,:),intent(in) :: a 308 316 complex(8),dimension(:,:),intent(in) :: a
complex(8),dimension(size(a,2),size(a,1)) :: fvn_op_z_conj_transpose 309 317 complex(8),dimension(size(a,2),size(a,1)) :: fvn_op_z_conj_transpose
fvn_op_z_conj_transpose=transpose(conjg(a)) 310 318 fvn_op_z_conj_transpose=transpose(conjg(a))
end function 311 319 end function
312 320
313 321
! 314 322 !
! .hx. 315 323 ! .hx.
! 316 324 !
function fvn_op_c_hx(a,b) 317 325 function fvn_op_c_hx(a,b)
implicit none 318 326 implicit none
complex(4), dimension(:,:), intent(in) :: a,b 319 327 complex(4), dimension(:,:), intent(in) :: a,b
complex(4), dimension(size(a,2),size(b,2)) :: fvn_op_c_hx 320 328 complex(4), dimension(size(a,2),size(b,2)) :: fvn_op_c_hx
fvn_op_c_hx=matmul(transpose(conjg(a)),b) 321 329 fvn_op_c_hx=matmul(transpose(conjg(a)),b)
end function 322 330 end function
function fvn_op_z_hx(a,b) 323 331 function fvn_op_z_hx(a,b)
implicit none 324 332 implicit none
complex(8), dimension(:,:), intent(in) :: a,b 325 333 complex(8), dimension(:,:), intent(in) :: a,b
complex(8), dimension(size(a,2),size(b,2)) :: fvn_op_z_hx 326 334 complex(8), dimension(size(a,2),size(b,2)) :: fvn_op_z_hx
fvn_op_z_hx=matmul(transpose(conjg(a)),b) 327 335 fvn_op_z_hx=matmul(transpose(conjg(a)),b)
end function 328 336 end function
329 337
330 338
! 331 339 !
! .xh. 332 340 ! .xh.
! 333 341 !
function fvn_op_c_xh(a,b) 334 342 function fvn_op_c_xh(a,b)
implicit none 335 343 implicit none
complex(4), dimension(:,:), intent(in) :: a,b 336 344 complex(4), dimension(:,:), intent(in) :: a,b
complex(4), dimension(size(a,1),size(b,1)) :: fvn_op_c_xh 337 345 complex(4), dimension(size(a,1),size(b,1)) :: fvn_op_c_xh
fvn_op_c_xh=matmul(a,transpose(conjg(b))) 338 346 fvn_op_c_xh=matmul(a,transpose(conjg(b)))
end function 339 347 end function
function fvn_op_z_xh(a,b) 340 348 function fvn_op_z_xh(a,b)
implicit none 341 349 implicit none
complex(8), dimension(:,:), intent(in) :: a,b 342 350 complex(8), dimension(:,:), intent(in) :: a,b
complex(8), dimension(size(a,1),size(b,1)) :: fvn_op_z_xh 343 351 complex(8), dimension(size(a,1),size(b,1)) :: fvn_op_z_xh
fvn_op_z_xh=matmul(a,transpose(conjg(b))) 344 352 fvn_op_z_xh=matmul(a,transpose(conjg(b)))
end function 345 353 end function
346 354
347 355
348 356
349 357
350 358
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 351 359 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 352 360 !
! Identity Matrix 353 361 ! Identity Matrix
! 354 362 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 355 363 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function fvn_d_ident(n) 356 364 function fvn_d_ident(n)
implicit none 357 365 implicit none
integer(kind=4) :: n 358 366 integer(kind=4) :: n
real(kind=8), dimension(n,n) :: fvn_d_ident 359 367 real(kind=8), dimension(n,n) :: fvn_d_ident
360 368
real(kind=8),dimension(n*n) :: vect 361 369 real(kind=8),dimension(n*n) :: vect
integer(kind=4) :: i 362 370 integer(kind=4) :: i
363 371
vect=0._8 364 372 vect=0._8
vect(1:n*n:n+1) = 1._8 365 373 vect(1:n*n:n+1) = 1._8
fvn_d_ident=reshape(vect, shape = (/ n,n /)) 366 374 fvn_d_ident=reshape(vect, shape = (/ n,n /))
end function 367 375 end function
368 376
function fvn_s_ident(n) 369 377 function fvn_s_ident(n)
implicit none 370 378 implicit none
integer(kind=4) :: n 371 379 integer(kind=4) :: n
real(kind=4), dimension(n,n) :: fvn_s_ident 372 380 real(kind=4), dimension(n,n) :: fvn_s_ident
373 381
real(kind=4),dimension(n*n) :: vect 374 382 real(kind=4),dimension(n*n) :: vect
integer(kind=4) :: i 375 383 integer(kind=4) :: i
376 384
vect=0._4 377 385 vect=0._4
vect(1:n*n:n+1) = 1._4 378 386 vect(1:n*n:n+1) = 1._4
fvn_s_ident=reshape(vect, shape = (/ n,n /)) 379 387 fvn_s_ident=reshape(vect, shape = (/ n,n /))
end function 380 388 end function
381 389
function fvn_c_ident(n) 382 390 function fvn_c_ident(n)
implicit none 383 391 implicit none
integer(kind=4) :: n 384 392 integer(kind=4) :: n
complex(kind=4), dimension(n,n) :: fvn_c_ident 385 393 complex(kind=4), dimension(n,n) :: fvn_c_ident
386 394
complex(kind=4),dimension(n*n) :: vect 387 395 complex(kind=4),dimension(n*n) :: vect
integer(kind=4) :: i 388 396 integer(kind=4) :: i
389 397
vect=(0._4,0._4) 390 398 vect=(0._4,0._4)
vect(1:n*n:n+1) = (1._4,0._4) 391 399 vect(1:n*n:n+1) = (1._4,0._4)
fvn_c_ident=reshape(vect, shape = (/ n,n /)) 392 400 fvn_c_ident=reshape(vect, shape = (/ n,n /))
end function 393 401 end function
394 402
function fvn_z_ident(n) 395 403 function fvn_z_ident(n)
implicit none 396 404 implicit none
integer(kind=4) :: n 397 405 integer(kind=4) :: n
complex(kind=8), dimension(n,n) :: fvn_z_ident 398 406 complex(kind=8), dimension(n,n) :: fvn_z_ident
399 407
complex(kind=8),dimension(n*n) :: vect 400 408 complex(kind=8),dimension(n*n) :: vect
integer(kind=4) :: i 401 409 integer(kind=4) :: i
402 410
vect=(0._8,0._8) 403 411 vect=(0._8,0._8)
vect(1:n*n:n+1) = (1._8,0._8) 404 412 vect(1:n*n:n+1) = (1._8,0._8)
fvn_z_ident=reshape(vect, shape = (/ n,n /)) 405 413 fvn_z_ident=reshape(vect, shape = (/ n,n /))
end function 406 414 end function
407 415
408 416
409 417
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 410 418 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 411 419 !
! Matrix inversion subroutines 412 420 ! Matrix inversion subroutines
! 413 421 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 414 422 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
415 423
subroutine fvn_s_matinv(d,a,inva,status) 416 424 subroutine fvn_s_matinv(d,a,inva,status)
! 417 425 !
! Matrix inversion of a real matrix using BLAS and LAPACK 418 426 ! Matrix inversion of a real matrix using BLAS and LAPACK
! 419 427 !
! d (in) : matrix rank 420 428 ! d (in) : matrix rank
! a (in) : input matrix 421 429 ! a (in) : input matrix
! inva (out) : inversed matrix 422 430 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 423 431 ! status (ou) : =0 if something failed
! 424 432 !
implicit none 425 433 implicit none
integer, intent(in) :: d 426 434 integer, intent(in) :: d
real, intent(in) :: a(d,d) 427 435 real, intent(in) :: a(d,d)
real, intent(out) :: inva(d,d) 428 436 real, intent(out) :: inva(d,d)
integer, intent(out),optional :: status 429 437 integer, intent(out),optional :: status
430 438
integer, allocatable :: ipiv(:) 431 439 integer, allocatable :: ipiv(:)
real, allocatable :: work(:) 432 440 real, allocatable :: work(:)
real twork(1) 433 441 real twork(1)
integer :: info 434 442 integer :: info
integer :: lwork 435 443 integer :: lwork
436 444
if (present(status)) status=1 437 445 if (present(status)) status=1
438 446
allocate(ipiv(d)) 439 447 allocate(ipiv(d))
! copy a into inva using BLAS 440 448 ! copy a into inva using BLAS
!call scopy(d*d,a,1,inva,1) 441 449 !call scopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 442 450 inva(:,:)=a(:,:)
! LU factorization using LAPACK 443 451 ! LU factorization using LAPACK
call sgetrf(d,d,inva,d,ipiv,info) 444 452 call sgetrf(d,d,inva,d,ipiv,info)
! if info is not equal to 0, something went wrong we exit setting status to 0 445 453 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 446 454 if (info /= 0) then
if (present(status)) status=0 447 455 if (present(status)) status=0
deallocate(ipiv) 448 456 deallocate(ipiv)
return 449 457 return
end if 450 458 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 451 459 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call sgetri(d,inva,d,ipiv,twork,-1,info) 452 460 call sgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 453 461 lwork=int(twork(1))
allocate(work(lwork)) 454 462 allocate(work(lwork))
! Matrix inversion using LAPACK 455 463 ! Matrix inversion using LAPACK
call sgetri(d,inva,d,ipiv,work,lwork,info) 456 464 call sgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 457 465 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 458 466 if (info /= 0) then
if (present(status)) status=0 459 467 if (present(status)) status=0
end if 460 468 end if
deallocate(work) 461 469 deallocate(work)
deallocate(ipiv) 462 470 deallocate(ipiv)
end subroutine 463 471 end subroutine
464 472
subroutine fvn_d_matinv(d,a,inva,status) 465 473 subroutine fvn_d_matinv(d,a,inva,status)
! 466 474 !
! Matrix inversion of a double precision matrix using BLAS and LAPACK 467 475 ! Matrix inversion of a double precision matrix using BLAS and LAPACK
! 468 476 !
! d (in) : matrix rank 469 477 ! d (in) : matrix rank
! a (in) : input matrix 470 478 ! a (in) : input matrix
! inva (out) : inversed matrix 471 479 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 472 480 ! status (ou) : =0 if something failed
! 473 481 !
implicit none 474 482 implicit none
integer, intent(in) :: d 475 483 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 476 484 double precision, intent(in) :: a(d,d)
double precision, intent(out) :: inva(d,d) 477 485 double precision, intent(out) :: inva(d,d)
integer, intent(out),optional :: status 478 486 integer, intent(out),optional :: status
479 487
integer, allocatable :: ipiv(:) 480 488 integer, allocatable :: ipiv(:)
double precision, allocatable :: work(:) 481 489 double precision, allocatable :: work(:)
double precision :: twork(1) 482 490 double precision :: twork(1)
integer :: info 483 491 integer :: info
integer :: lwork 484 492 integer :: lwork
485 493
if (present(status)) status=1 486 494 if (present(status)) status=1
487 495
allocate(ipiv(d)) 488 496 allocate(ipiv(d))
! copy a into inva using BLAS 489 497 ! copy a into inva using BLAS
!call dcopy(d*d,a,1,inva,1) 490 498 !call dcopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 491 499 inva(:,:)=a(:,:)
! LU factorization using LAPACK 492 500 ! LU factorization using LAPACK
call dgetrf(d,d,inva,d,ipiv,info) 493 501 call dgetrf(d,d,inva,d,ipiv,info)
! if info is not equal to 0, something went wrong we exit setting status to 0 494 502 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 495 503 if (info /= 0) then
if (present(status)) status=0 496 504 if (present(status)) status=0
deallocate(ipiv) 497 505 deallocate(ipiv)
return 498 506 return
end if 499 507 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 500 508 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call dgetri(d,inva,d,ipiv,twork,-1,info) 501 509 call dgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 502 510 lwork=int(twork(1))
allocate(work(lwork)) 503 511 allocate(work(lwork))
! Matrix inversion using LAPACK 504 512 ! Matrix inversion using LAPACK
call dgetri(d,inva,d,ipiv,work,lwork,info) 505 513 call dgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 506 514 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 507 515 if (info /= 0) then
if (present(status)) status=0 508 516 if (present(status)) status=0
end if 509 517 end if
deallocate(work) 510 518 deallocate(work)
deallocate(ipiv) 511 519 deallocate(ipiv)
end subroutine 512 520 end subroutine
513 521
subroutine fvn_c_matinv(d,a,inva,status) 514 522 subroutine fvn_c_matinv(d,a,inva,status)
! 515 523 !
! Matrix inversion of a complex matrix using BLAS and LAPACK 516 524 ! Matrix inversion of a complex matrix using BLAS and LAPACK
! 517 525 !
! d (in) : matrix rank 518 526 ! d (in) : matrix rank
! a (in) : input matrix 519 527 ! a (in) : input matrix
! inva (out) : inversed matrix 520 528 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 521 529 ! status (ou) : =0 if something failed
! 522 530 !
implicit none 523 531 implicit none
integer, intent(in) :: d 524 532 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 525 533 complex, intent(in) :: a(d,d)
complex, intent(out) :: inva(d,d) 526 534 complex, intent(out) :: inva(d,d)
integer, intent(out),optional :: status 527 535 integer, intent(out),optional :: status
528 536
integer, allocatable :: ipiv(:) 529 537 integer, allocatable :: ipiv(:)
complex, allocatable :: work(:) 530 538 complex, allocatable :: work(:)
complex :: twork(1) 531 539 complex :: twork(1)
integer :: info 532 540 integer :: info
integer :: lwork 533 541 integer :: lwork
534 542
if (present(status)) status=1 535 543 if (present(status)) status=1
536 544
allocate(ipiv(d)) 537 545 allocate(ipiv(d))
! copy a into inva using BLAS 538 546 ! copy a into inva using BLAS
!call ccopy(d*d,a,1,inva,1) 539 547 !call ccopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 540 548 inva(:,:)=a(:,:)
541 549
! LU factorization using LAPACK 542 550 ! LU factorization using LAPACK
call cgetrf(d,d,inva,d,ipiv,info) 543 551 call cgetrf(d,d,inva,d,ipiv,info)
! if info is not equal to 0, something went wrong we exit setting status to 0 544 552 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 545 553 if (info /= 0) then
if (present(status)) status=0 546 554 if (present(status)) status=0
deallocate(ipiv) 547 555 deallocate(ipiv)
return 548 556 return
end if 549 557 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 550 558 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call cgetri(d,inva,d,ipiv,twork,-1,info) 551 559 call cgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 552 560 lwork=int(twork(1))
allocate(work(lwork)) 553 561 allocate(work(lwork))
! Matrix inversion using LAPACK 554 562 ! Matrix inversion using LAPACK
call cgetri(d,inva,d,ipiv,work,lwork,info) 555 563 call cgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 556 564 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 557 565 if (info /= 0) then
if (present(status)) status=0 558 566 if (present(status)) status=0
end if 559 567 end if
deallocate(work) 560 568 deallocate(work)
deallocate(ipiv) 561 569 deallocate(ipiv)
end subroutine 562 570 end subroutine
563 571
subroutine fvn_z_matinv(d,a,inva,status) 564 572 subroutine fvn_z_matinv(d,a,inva,status)
! 565 573 !
! Matrix inversion of a double complex matrix using BLAS and LAPACK 566 574 ! Matrix inversion of a double complex matrix using BLAS and LAPACK
! 567 575 !
! d (in) : matrix rank 568 576 ! d (in) : matrix rank
! a (in) : input matrix 569 577 ! a (in) : input matrix
! inva (out) : inversed matrix 570 578 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 571 579 ! status (ou) : =0 if something failed
! 572 580 !
implicit none 573 581 implicit none
integer, intent(in) :: d 574 582 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 575 583 double complex, intent(in) :: a(d,d)
double complex, intent(out) :: inva(d,d) 576 584 double complex, intent(out) :: inva(d,d)
integer, intent(out),optional :: status 577 585 integer, intent(out),optional :: status
578 586
integer, allocatable :: ipiv(:) 579 587 integer, allocatable :: ipiv(:)
double complex, allocatable :: work(:) 580 588 double complex, allocatable :: work(:)
double complex :: twork(1) 581 589 double complex :: twork(1)
integer :: info 582 590 integer :: info
integer :: lwork 583 591 integer :: lwork
584 592
if (present(status)) status=1 585 593 if (present(status)) status=1
586 594
allocate(ipiv(d)) 587 595 allocate(ipiv(d))
! copy a into inva using BLAS 588 596 ! copy a into inva using BLAS
!call zcopy(d*d,a,1,inva,1) 589 597 !call zcopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 590 598 inva(:,:)=a(:,:)
591 599
! LU factorization using LAPACK 592 600 ! LU factorization using LAPACK
call zgetrf(d,d,inva,d,ipiv,info) 593 601 call zgetrf(d,d,inva,d,ipiv,info)
! if info is not equal to 0, something went wrong we exit setting status to 0 594 602 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 595 603 if (info /= 0) then
if (present(status)) status=0 596 604 if (present(status)) status=0
deallocate(ipiv) 597 605 deallocate(ipiv)
return 598 606 return
end if 599 607 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 600 608 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call zgetri(d,inva,d,ipiv,twork,-1,info) 601 609 call zgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 602 610 lwork=int(twork(1))
allocate(work(lwork)) 603 611 allocate(work(lwork))
! Matrix inversion using LAPACK 604 612 ! Matrix inversion using LAPACK
call zgetri(d,inva,d,ipiv,work,lwork,info) 605 613 call zgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 606 614 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 607 615 if (info /= 0) then
if (present(status)) status=0 608 616 if (present(status)) status=0
end if 609 617 end if
deallocate(work) 610 618 deallocate(work)
deallocate(ipiv) 611 619 deallocate(ipiv)
end subroutine 612 620 end subroutine
613 621
614 622
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 615 623 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 616 624 !
! Determinants 617 625 ! Determinants
! 618 626 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 619 627 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function fvn_s_det(d,a,status) 620 628 function fvn_s_det(d,a,status)
! 621 629 !
! Evaluate the determinant of a square matrix using lapack LU factorization 622 630 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 623 631 !
! d (in) : matrix rank 624 632 ! d (in) : matrix rank
! a (in) : The Matrix 625 633 ! a (in) : The Matrix
! status (out) : =0 if LU factorization failed 626 634 ! status (out) : =0 if LU factorization failed
! 627 635 !
implicit none 628 636 implicit none
integer, intent(in) :: d 629 637 integer, intent(in) :: d
real, intent(in) :: a(d,d) 630 638 real, intent(in) :: a(d,d)
integer, intent(out), optional :: status 631 639 integer, intent(out), optional :: status
real :: fvn_s_det 632 640 real :: fvn_s_det
633 641
real, allocatable :: wc_a(:,:) 634 642 real, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 635 643 integer, allocatable :: ipiv(:)
integer :: info,i 636 644 integer :: info,i
637 645
if (present(status)) status=1 638 646 if (present(status)) status=1
allocate(wc_a(d,d)) 639 647 allocate(wc_a(d,d))
allocate(ipiv(d)) 640 648 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 641 649 wc_a(:,:)=a(:,:)
call sgetrf(d,d,wc_a,d,ipiv,info) 642 650 call sgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 643 651 if (info/= 0) then
if (present(status)) status=0 644 652 if (present(status)) status=0
fvn_s_det=0.e0 645 653 fvn_s_det=0.e0
deallocate(ipiv) 646 654 deallocate(ipiv)
deallocate(wc_a) 647 655 deallocate(wc_a)
return 648 656 return
end if 649 657 end if
fvn_s_det=1.e0 650 658 fvn_s_det=1.e0
do i=1,d 651 659 do i=1,d
if (ipiv(i)==i) then 652 660 if (ipiv(i)==i) then
fvn_s_det=fvn_s_det*wc_a(i,i) 653 661 fvn_s_det=fvn_s_det*wc_a(i,i)
else 654 662 else
fvn_s_det=-fvn_s_det*wc_a(i,i) 655 663 fvn_s_det=-fvn_s_det*wc_a(i,i)
end if 656 664 end if
end do 657 665 end do
deallocate(ipiv) 658 666 deallocate(ipiv)
deallocate(wc_a) 659 667 deallocate(wc_a)
660 668
end function 661 669 end function
662 670
function fvn_d_det(d,a,status) 663 671 function fvn_d_det(d,a,status)
! 664 672 !
! Evaluate the determinant of a square matrix using lapack LU factorization 665 673 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 666 674 !
! d (in) : matrix rank 667 675 ! d (in) : matrix rank
! a (in) : The Matrix 668 676 ! a (in) : The Matrix
! status (out) : =0 if LU factorization failed 669 677 ! status (out) : =0 if LU factorization failed
! 670 678 !
implicit none 671 679 implicit none
integer, intent(in) :: d 672 680 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 673 681 double precision, intent(in) :: a(d,d)
integer, intent(out), optional :: status 674 682 integer, intent(out), optional :: status
double precision :: fvn_d_det 675 683 double precision :: fvn_d_det
676 684
double precision, allocatable :: wc_a(:,:) 677 685 double precision, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 678 686 integer, allocatable :: ipiv(:)
integer :: info,i 679 687 integer :: info,i
680 688
if (present(status)) status=1 681 689 if (present(status)) status=1
allocate(wc_a(d,d)) 682 690 allocate(wc_a(d,d))
allocate(ipiv(d)) 683 691 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 684 692 wc_a(:,:)=a(:,:)
call dgetrf(d,d,wc_a,d,ipiv,info) 685 693 call dgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 686 694 if (info/= 0) then
if (present(status)) status=0 687 695 if (present(status)) status=0
fvn_d_det=0.d0 688 696 fvn_d_det=0.d0
deallocate(ipiv) 689 697 deallocate(ipiv)
deallocate(wc_a) 690 698 deallocate(wc_a)
return 691 699 return
end if 692 700 end if
fvn_d_det=1.d0 693 701 fvn_d_det=1.d0
do i=1,d 694 702 do i=1,d
if (ipiv(i)==i) then 695 703 if (ipiv(i)==i) then
fvn_d_det=fvn_d_det*wc_a(i,i) 696 704 fvn_d_det=fvn_d_det*wc_a(i,i)
else 697 705 else
fvn_d_det=-fvn_d_det*wc_a(i,i) 698 706 fvn_d_det=-fvn_d_det*wc_a(i,i)
end if 699 707 end if
end do 700 708 end do
deallocate(ipiv) 701 709 deallocate(ipiv)
deallocate(wc_a) 702 710 deallocate(wc_a)
703 711
end function 704 712 end function
705 713
function fvn_c_det(d,a,status) ! 706 714 function fvn_c_det(d,a,status) !
! Evaluate the determinant of a square matrix using lapack LU factorization 707 715 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 708 716 !
! d (in) : matrix rank 709 717 ! d (in) : matrix rank
! a (in) : The Matrix 710 718 ! a (in) : The Matrix
! status (out) : =0 if LU factorization failed 711 719 ! status (out) : =0 if LU factorization failed
! 712 720 !
implicit none 713 721 implicit none
integer, intent(in) :: d 714 722 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 715 723 complex, intent(in) :: a(d,d)
integer, intent(out), optional :: status 716 724 integer, intent(out), optional :: status
complex :: fvn_c_det 717 725 complex :: fvn_c_det
718 726
complex, allocatable :: wc_a(:,:) 719 727 complex, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 720 728 integer, allocatable :: ipiv(:)
integer :: info,i 721 729 integer :: info,i
722 730
if (present(status)) status=1 723 731 if (present(status)) status=1
allocate(wc_a(d,d)) 724 732 allocate(wc_a(d,d))
allocate(ipiv(d)) 725 733 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 726 734 wc_a(:,:)=a(:,:)
call cgetrf(d,d,wc_a,d,ipiv,info) 727 735 call cgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 728 736 if (info/= 0) then
if (present(status)) status=0 729 737 if (present(status)) status=0
fvn_c_det=(0.e0,0.e0) 730 738 fvn_c_det=(0.e0,0.e0)
deallocate(ipiv) 731 739 deallocate(ipiv)
deallocate(wc_a) 732 740 deallocate(wc_a)
return 733 741 return
end if 734 742 end if
fvn_c_det=(1.e0,0.e0) 735 743 fvn_c_det=(1.e0,0.e0)
do i=1,d 736 744 do i=1,d
if (ipiv(i)==i) then 737 745 if (ipiv(i)==i) then
fvn_c_det=fvn_c_det*wc_a(i,i) 738 746 fvn_c_det=fvn_c_det*wc_a(i,i)
else 739 747 else
fvn_c_det=-fvn_c_det*wc_a(i,i) 740 748 fvn_c_det=-fvn_c_det*wc_a(i,i)
end if 741 749 end if
end do 742 750 end do
deallocate(ipiv) 743 751 deallocate(ipiv)
deallocate(wc_a) 744 752 deallocate(wc_a)
745 753
end function 746 754 end function
747 755
function fvn_z_det(d,a,status) 748 756 function fvn_z_det(d,a,status)
! 749 757 !
! Evaluate the determinant of a square matrix using lapack LU factorization 750 758 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 751 759 !
! d (in) : matrix rank 752 760 ! d (in) : matrix rank
! a (in) : The Matrix 753 761 ! a (in) : The Matrix
! det (out) : determinant 754 762 ! det (out) : determinant
! status (out) : =0 if LU factorization failed 755 763 ! status (out) : =0 if LU factorization failed
! 756 764 !
implicit none 757 765 implicit none
integer, intent(in) :: d 758 766 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 759 767 double complex, intent(in) :: a(d,d)
integer, intent(out), optional :: status 760 768 integer, intent(out), optional :: status
double complex :: fvn_z_det 761 769 double complex :: fvn_z_det
762 770
double complex, allocatable :: wc_a(:,:) 763 771 double complex, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 764 772 integer, allocatable :: ipiv(:)
integer :: info,i 765 773 integer :: info,i
766 774
if (present(status)) status=1 767 775 if (present(status)) status=1
allocate(wc_a(d,d)) 768 776 allocate(wc_a(d,d))
allocate(ipiv(d)) 769 777 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 770 778 wc_a(:,:)=a(:,:)
call zgetrf(d,d,wc_a,d,ipiv,info) 771 779 call zgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 772 780 if (info/= 0) then
if (present(status)) status=0 773 781 if (present(status)) status=0
fvn_z_det=(0.d0,0.d0) 774 782 fvn_z_det=(0.d0,0.d0)
deallocate(ipiv) 775 783 deallocate(ipiv)
deallocate(wc_a) 776 784 deallocate(wc_a)
return 777 785 return
end if 778 786 end if
fvn_z_det=(1.d0,0.d0) 779 787 fvn_z_det=(1.d0,0.d0)
do i=1,d 780 788 do i=1,d
if (ipiv(i)==i) then 781 789 if (ipiv(i)==i) then
fvn_z_det=fvn_z_det*wc_a(i,i) 782 790 fvn_z_det=fvn_z_det*wc_a(i,i)
else 783 791 else
fvn_z_det=-fvn_z_det*wc_a(i,i) 784 792 fvn_z_det=-fvn_z_det*wc_a(i,i)
end if 785 793 end if
end do 786 794 end do
deallocate(ipiv) 787 795 deallocate(ipiv)
deallocate(wc_a) 788 796 deallocate(wc_a)
789 797
end function 790 798 end function
791 799
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 792 800 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 793 801 !
! Condition test 794 802 ! Condition test
! 795 803 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 796 804 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1-norm 797 805 ! 1-norm
! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm 798 806 ! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm
! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond 799 807 ! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond
! 800 808 !
subroutine fvn_s_matcon(d,a,rcond,status) 801 809 subroutine fvn_s_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 802 810 ! Matrix condition (reciprocal of condition number)
! 803 811 !
! d (in) : matrix rank 804 812 ! d (in) : matrix rank
! a (in) : The Matrix 805 813 ! a (in) : The Matrix
! rcond (out) : guess what 806 814 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 807 815 ! status (out) : =0 if something went wrong
! 808 816 !
implicit none 809 817 implicit none
integer, intent(in) :: d 810 818 integer, intent(in) :: d
real, intent(in) :: a(d,d) 811 819 real, intent(in) :: a(d,d)
real, intent(out) :: rcond 812 820 real, intent(out) :: rcond
integer, intent(out), optional :: status 813 821 integer, intent(out), optional :: status
814 822
real, allocatable :: work(:) 815 823 real, allocatable :: work(:)
integer, allocatable :: iwork(:) 816 824 integer, allocatable :: iwork(:)
real :: anorm 817 825 real :: anorm
real, allocatable :: wc_a(:,:) ! working copy of a 818 826 real, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 819 827 integer :: info
integer, allocatable :: ipiv(:) 820 828 integer, allocatable :: ipiv(:)
821 829
real, external :: slange 822 830 real, external :: slange
823 831
824 832
if (present(status)) status=1 825 833 if (present(status)) status=1
826 834
anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm 827 835 anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
828 836
allocate(wc_a(d,d)) 829 837 allocate(wc_a(d,d))
!call scopy(d*d,a,1,wc_a,1) 830 838 !call scopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 831 839 wc_a(:,:)=a(:,:)
832 840
allocate(ipiv(d)) 833 841 allocate(ipiv(d))
call sgetrf(d,d,wc_a,d,ipiv,info) 834 842 call sgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 835 843 if (info /= 0) then
if (present(status)) status=0 836 844 if (present(status)) status=0
deallocate(ipiv) 837 845 deallocate(ipiv)
deallocate(wc_a) 838 846 deallocate(wc_a)
return 839 847 return
end if 840 848 end if
allocate(work(4*d)) 841 849 allocate(work(4*d))
allocate(iwork(d)) 842 850 allocate(iwork(d))
call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info) 843 851 call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
if (info /= 0) then 844 852 if (info /= 0) then
if (present(status)) status=0 845 853 if (present(status)) status=0
end if 846 854 end if
deallocate(iwork) 847 855 deallocate(iwork)
deallocate(work) 848 856 deallocate(work)
deallocate(ipiv) 849 857 deallocate(ipiv)
deallocate(wc_a) 850 858 deallocate(wc_a)
851 859
end subroutine 852 860 end subroutine
853 861
subroutine fvn_d_matcon(d,a,rcond,status) 854 862 subroutine fvn_d_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 855 863 ! Matrix condition (reciprocal of condition number)
! 856 864 !
! d (in) : matrix rank 857 865 ! d (in) : matrix rank
! a (in) : The Matrix 858 866 ! a (in) : The Matrix
! rcond (out) : guess what 859 867 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 860 868 ! status (out) : =0 if something went wrong
! 861 869 !
implicit none 862 870 implicit none
integer, intent(in) :: d 863 871 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 864 872 double precision, intent(in) :: a(d,d)
double precision, intent(out) :: rcond 865 873 double precision, intent(out) :: rcond
integer, intent(out), optional :: status 866 874 integer, intent(out), optional :: status
867 875
double precision, allocatable :: work(:) 868 876 double precision, allocatable :: work(:)
integer, allocatable :: iwork(:) 869 877 integer, allocatable :: iwork(:)
double precision :: anorm 870 878 double precision :: anorm
double precision, allocatable :: wc_a(:,:) ! working copy of a 871 879 double precision, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 872 880 integer :: info
integer, allocatable :: ipiv(:) 873 881 integer, allocatable :: ipiv(:)
874 882
double precision, external :: dlange 875 883 double precision, external :: dlange
876 884
877 885
if (present(status)) status=1 878 886 if (present(status)) status=1
879 887
anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm 880 888 anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
881 889
allocate(wc_a(d,d)) 882 890 allocate(wc_a(d,d))
!call dcopy(d*d,a,1,wc_a,1) 883 891 !call dcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 884 892 wc_a(:,:)=a(:,:)
885 893
allocate(ipiv(d)) 886 894 allocate(ipiv(d))
call dgetrf(d,d,wc_a,d,ipiv,info) 887 895 call dgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 888 896 if (info /= 0) then
if (present(status)) status=0 889 897 if (present(status)) status=0
deallocate(ipiv) 890 898 deallocate(ipiv)
deallocate(wc_a) 891 899 deallocate(wc_a)
return 892 900 return
end if 893 901 end if
894 902
allocate(work(4*d)) 895 903 allocate(work(4*d))
allocate(iwork(d)) 896 904 allocate(iwork(d))
call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info) 897 905 call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
if (info /= 0) then 898 906 if (info /= 0) then
if (present(status)) status=0 899 907 if (present(status)) status=0
end if 900 908 end if
deallocate(iwork) 901 909 deallocate(iwork)
deallocate(work) 902 910 deallocate(work)
deallocate(ipiv) 903 911 deallocate(ipiv)
deallocate(wc_a) 904 912 deallocate(wc_a)
905 913
end subroutine 906 914 end subroutine
907 915
subroutine fvn_c_matcon(d,a,rcond,status) 908 916 subroutine fvn_c_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 909 917 ! Matrix condition (reciprocal of condition number)
! 910 918 !
! d (in) : matrix rank 911 919 ! d (in) : matrix rank
! a (in) : The Matrix 912 920 ! a (in) : The Matrix
! rcond (out) : guess what 913 921 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 914 922 ! status (out) : =0 if something went wrong
! 915 923 !
implicit none 916 924 implicit none
integer, intent(in) :: d 917 925 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 918 926 complex, intent(in) :: a(d,d)
real, intent(out) :: rcond 919 927 real, intent(out) :: rcond
integer, intent(out), optional :: status 920 928 integer, intent(out), optional :: status
921 929
real, allocatable :: rwork(:) 922 930 real, allocatable :: rwork(:)
complex, allocatable :: work(:) 923 931 complex, allocatable :: work(:)
integer, allocatable :: iwork(:) 924 932 integer, allocatable :: iwork(:)
real :: anorm 925 933 real :: anorm
complex, allocatable :: wc_a(:,:) ! working copy of a 926 934 complex, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 927 935 integer :: info
integer, allocatable :: ipiv(:) 928 936 integer, allocatable :: ipiv(:)
929 937
real, external :: clange 930 938 real, external :: clange
931 939
932 940
if (present(status)) status=1 933 941 if (present(status)) status=1
934 942
anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm 935 943 anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
936 944
allocate(wc_a(d,d)) 937 945 allocate(wc_a(d,d))
!call ccopy(d*d,a,1,wc_a,1) 938 946 !call ccopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 939 947 wc_a(:,:)=a(:,:)
940 948
allocate(ipiv(d)) 941 949 allocate(ipiv(d))
call cgetrf(d,d,wc_a,d,ipiv,info) 942 950 call cgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 943 951 if (info /= 0) then
if (present(status)) status=0 944 952 if (present(status)) status=0
deallocate(ipiv) 945 953 deallocate(ipiv)
deallocate(wc_a) 946 954 deallocate(wc_a)
return 947 955 return
end if 948 956 end if
allocate(work(2*d)) 949 957 allocate(work(2*d))
allocate(rwork(2*d)) 950 958 allocate(rwork(2*d))
call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info) 951 959 call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
if (info /= 0) then 952 960 if (info /= 0) then
if (present(status)) status=0 953 961 if (present(status)) status=0
end if 954 962 end if
deallocate(rwork) 955 963 deallocate(rwork)
deallocate(work) 956 964 deallocate(work)
deallocate(ipiv) 957 965 deallocate(ipiv)
deallocate(wc_a) 958 966 deallocate(wc_a)
end subroutine 959 967 end subroutine
960 968
subroutine fvn_z_matcon(d,a,rcond,status) 961 969 subroutine fvn_z_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 962 970 ! Matrix condition (reciprocal of condition number)
! 963 971 !
! d (in) : matrix rank 964 972 ! d (in) : matrix rank
! a (in) : The Matrix 965 973 ! a (in) : The Matrix
! rcond (out) : guess what 966 974 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 967 975 ! status (out) : =0 if something went wrong
! 968 976 !
implicit none 969 977 implicit none
integer, intent(in) :: d 970 978 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 971 979 double complex, intent(in) :: a(d,d)
double precision, intent(out) :: rcond 972 980 double precision, intent(out) :: rcond
integer, intent(out), optional :: status 973 981 integer, intent(out), optional :: status
974 982
double complex, allocatable :: work(:) 975 983 double complex, allocatable :: work(:)
double precision, allocatable :: rwork(:) 976 984 double precision, allocatable :: rwork(:)
double precision :: anorm 977 985 double precision :: anorm
double complex, allocatable :: wc_a(:,:) ! working copy of a 978 986 double complex, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 979 987 integer :: info
integer, allocatable :: ipiv(:) 980 988 integer, allocatable :: ipiv(:)
981 989
double precision, external :: zlange 982 990 double precision, external :: zlange
983 991
984 992
if (present(status)) status=1 985 993 if (present(status)) status=1
986 994
anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm 987 995 anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
988 996
allocate(wc_a(d,d)) 989 997 allocate(wc_a(d,d))
!call zcopy(d*d,a,1,wc_a,1) 990 998 !call zcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 991 999 wc_a(:,:)=a(:,:)
992 1000
allocate(ipiv(d)) 993 1001 allocate(ipiv(d))
call zgetrf(d,d,wc_a,d,ipiv,info) 994 1002 call zgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 995 1003 if (info /= 0) then
if (present(status)) status=0 996 1004 if (present(status)) status=0
deallocate(ipiv) 997 1005 deallocate(ipiv)
deallocate(wc_a) 998 1006 deallocate(wc_a)
return 999 1007 return
end if 1000 1008 end if
1001 1009
allocate(work(2*d)) 1002 1010 allocate(work(2*d))
allocate(rwork(2*d)) 1003 1011 allocate(rwork(2*d))
call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info) 1004 1012 call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
if (info /= 0) then 1005 1013 if (info /= 0) then
if (present(status)) status=0 1006 1014 if (present(status)) status=0
end if 1007 1015 end if
deallocate(rwork) 1008 1016 deallocate(rwork)
deallocate(work) 1009 1017 deallocate(work)
deallocate(ipiv) 1010 1018 deallocate(ipiv)
deallocate(wc_a) 1011 1019 deallocate(wc_a)
end subroutine 1012 1020 end subroutine
1013 1021
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1014 1022 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1015 1023 !
! Valeurs propres/ Vecteurs propre 1016 1024 ! Valeurs propres/ Vecteurs propre
! 1017 1025 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1018 1026 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1019 1027
subroutine fvn_s_matev(d,a,evala,eveca,status) 1020 1028 subroutine fvn_s_matev(d,a,evala,eveca,status)
! 1021 1029 !
! integer d (in) : matrice rank 1022 1030 ! integer d (in) : matrice rank
! real a(d,d) (in) : The Matrix 1023 1031 ! real a(d,d) (in) : The Matrix
! complex evala(d) (out) : eigenvalues 1024 1032 ! complex evala(d) (out) : eigenvalues
! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1025 1033 ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 1026 1034 ! integer (out) : status =0 if something went wrong
! 1027 1035 !
! interfacing Lapack routine SGEEV 1028 1036 ! interfacing Lapack routine SGEEV
implicit none 1029 1037 implicit none
integer, intent(in) :: d 1030 1038 integer, intent(in) :: d
real, intent(in) :: a(d,d) 1031 1039 real, intent(in) :: a(d,d)
complex, intent(out) :: evala(d) 1032 1040 complex, intent(out) :: evala(d)
complex, intent(out) :: eveca(d,d) 1033 1041 complex, intent(out) :: eveca(d,d)
integer, intent(out), optional :: status 1034 1042 integer, intent(out), optional :: status
1035 1043
real, allocatable :: wc_a(:,:) ! a working copy of a 1036 1044 real, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 1037 1045 integer :: info
integer :: lwork 1038 1046 integer :: lwork
real, allocatable :: wr(:),wi(:) 1039 1047 real, allocatable :: wr(:),wi(:)
real :: vl ! unused but necessary for the call 1040 1048 real :: vl ! unused but necessary for the call
real, allocatable :: vr(:,:) 1041 1049 real, allocatable :: vr(:,:)
real, allocatable :: work(:) 1042 1050 real, allocatable :: work(:)
real :: twork(1) 1043 1051 real :: twork(1)
integer i 1044 1052 integer i
integer j 1045 1053 integer j
1046 1054
if (present(status)) status=1 1047 1055 if (present(status)) status=1
1048 1056
! making a working copy of a 1049 1057 ! making a working copy of a
allocate(wc_a(d,d)) 1050 1058 allocate(wc_a(d,d))
!call scopy(d*d,a,1,wc_a,1) 1051 1059 !call scopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1052 1060 wc_a(:,:)=a(:,:)
1053 1061
allocate(wr(d)) 1054 1062 allocate(wr(d))
allocate(wi(d)) 1055 1063 allocate(wi(d))
allocate(vr(d,d)) 1056 1064 allocate(vr(d,d))
! query optimal work size 1057 1065 ! query optimal work size
call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) 1058 1066 call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
lwork=int(twork(1)) 1059 1067 lwork=int(twork(1))
allocate(work(lwork)) 1060 1068 allocate(work(lwork))
call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) 1061 1069 call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
1062 1070
if (info /= 0) then 1063 1071 if (info /= 0) then
if (present(status)) status=0 1064 1072 if (present(status)) status=0
deallocate(work) 1065 1073 deallocate(work)
deallocate(vr) 1066 1074 deallocate(vr)
deallocate(wi) 1067 1075 deallocate(wi)
deallocate(wr) 1068 1076 deallocate(wr)
deallocate(wc_a) 1069 1077 deallocate(wc_a)
return 1070 1078 return
end if 1071 1079 end if
1072 1080
! now fill in the results 1073 1081 ! now fill in the results
i=1 1074 1082 i=1
do while(i<=d) 1075 1083 do while(i<=d)
evala(i)=cmplx(wr(i),wi(i)) 1076 1084 evala(i)=cmplx(wr(i),wi(i))
if (wi(i) == 0.) then ! eigenvalue is real 1077 1085 if (wi(i) == 0.) then ! eigenvalue is real
eveca(:,i)=cmplx(vr(:,i),0.) 1078 1086 eveca(:,i)=cmplx(vr(:,i),0.)
else ! eigenvalue is complex 1079 1087 else ! eigenvalue is complex
evala(i+1)=cmplx(wr(i+1),wi(i+1)) 1080 1088 evala(i+1)=cmplx(wr(i+1),wi(i+1))
eveca(:,i)=cmplx(vr(:,i),vr(:,i+1)) 1081 1089 eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))
eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1)) 1082 1090 eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))
i=i+1 1083 1091 i=i+1
end if 1084 1092 end if
i=i+1 1085 1093 i=i+1
enddo 1086 1094 enddo
deallocate(work) 1087 1095 deallocate(work)
deallocate(vr) 1088 1096 deallocate(vr)
deallocate(wi) 1089 1097 deallocate(wi)
deallocate(wr) 1090 1098 deallocate(wr)
deallocate(wc_a) 1091 1099 deallocate(wc_a)
1092 1100
end subroutine 1093 1101 end subroutine
1094 1102
subroutine fvn_d_matev(d,a,evala,eveca,status) 1095 1103 subroutine fvn_d_matev(d,a,evala,eveca,status)
! 1096 1104 !
! integer d (in) : matrice rank 1097 1105 ! integer d (in) : matrice rank
! double precision a(d,d) (in) : The Matrix 1098 1106 ! double precision a(d,d) (in) : The Matrix
! double complex evala(d) (out) : eigenvalues 1099 1107 ! double complex evala(d) (out) : eigenvalues
! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1100 1108 ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 1101 1109 ! integer (out) : status =0 if something went wrong
! 1102 1110 !
! interfacing Lapack routine DGEEV 1103 1111 ! interfacing Lapack routine DGEEV
implicit none 1104 1112 implicit none
integer, intent(in) :: d 1105 1113 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 1106 1114 double precision, intent(in) :: a(d,d)
double complex, intent(out) :: evala(d) 1107 1115 double complex, intent(out) :: evala(d)
double complex, intent(out) :: eveca(d,d) 1108 1116 double complex, intent(out) :: eveca(d,d)
integer, intent(out), optional :: status 1109 1117 integer, intent(out), optional :: status
1110 1118
double precision, allocatable :: wc_a(:,:) ! a working copy of a 1111 1119 double precision, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 1112 1120 integer :: info
integer :: lwork 1113 1121 integer :: lwork
double precision, allocatable :: wr(:),wi(:) 1114 1122 double precision, allocatable :: wr(:),wi(:)
double precision :: vl ! unused but necessary for the call 1115 1123 double precision :: vl ! unused but necessary for the call
double precision, allocatable :: vr(:,:) 1116 1124 double precision, allocatable :: vr(:,:)
double precision, allocatable :: work(:) 1117 1125 double precision, allocatable :: work(:)
double precision :: twork(1) 1118 1126 double precision :: twork(1)
integer i 1119 1127 integer i
integer j 1120 1128 integer j
1121 1129
if (present(status)) status=1 1122 1130 if (present(status)) status=1
1123 1131
! making a working copy of a 1124 1132 ! making a working copy of a
allocate(wc_a(d,d)) 1125 1133 allocate(wc_a(d,d))
!call dcopy(d*d,a,1,wc_a,1) 1126 1134 !call dcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1127 1135 wc_a(:,:)=a(:,:)
1128 1136
allocate(wr(d)) 1129 1137 allocate(wr(d))
allocate(wi(d)) 1130 1138 allocate(wi(d))
allocate(vr(d,d)) 1131 1139 allocate(vr(d,d))
! query optimal work size 1132 1140 ! query optimal work size
call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) 1133 1141 call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
lwork=int(twork(1)) 1134 1142 lwork=int(twork(1))
allocate(work(lwork)) 1135 1143 allocate(work(lwork))
call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) 1136 1144 call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
1137 1145
if (info /= 0) then 1138 1146 if (info /= 0) then
if (present(status)) status=0 1139 1147 if (present(status)) status=0
deallocate(work) 1140 1148 deallocate(work)
deallocate(vr) 1141 1149 deallocate(vr)
deallocate(wi) 1142 1150 deallocate(wi)
deallocate(wr) 1143 1151 deallocate(wr)
deallocate(wc_a) 1144 1152 deallocate(wc_a)
return 1145 1153 return
end if 1146 1154 end if
1147 1155
! now fill in the results 1148 1156 ! now fill in the results
i=1 1149 1157 i=1
do while(i<=d) 1150 1158 do while(i<=d)
evala(i)=dcmplx(wr(i),wi(i)) 1151 1159 evala(i)=dcmplx(wr(i),wi(i))
if (wi(i) == 0.) then ! eigenvalue is real 1152 1160 if (wi(i) == 0.) then ! eigenvalue is real
eveca(:,i)=dcmplx(vr(:,i),0.) 1153 1161 eveca(:,i)=dcmplx(vr(:,i),0.)
else ! eigenvalue is complex 1154 1162 else ! eigenvalue is complex
evala(i+1)=dcmplx(wr(i+1),wi(i+1)) 1155 1163 evala(i+1)=dcmplx(wr(i+1),wi(i+1))
eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1)) 1156 1164 eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))
eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1)) 1157 1165 eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))
i=i+1 1158 1166 i=i+1
end if 1159 1167 end if
i=i+1 1160 1168 i=i+1
enddo 1161 1169 enddo
1162 1170
deallocate(work) 1163 1171 deallocate(work)
deallocate(vr) 1164 1172 deallocate(vr)
deallocate(wi) 1165 1173 deallocate(wi)
deallocate(wr) 1166 1174 deallocate(wr)
deallocate(wc_a) 1167 1175 deallocate(wc_a)
1168 1176
end subroutine 1169 1177 end subroutine
1170 1178
subroutine fvn_c_matev(d,a,evala,eveca,status) 1171 1179 subroutine fvn_c_matev(d,a,evala,eveca,status)
! 1172 1180 !
! integer d (in) : matrice rank 1173 1181 ! integer d (in) : matrice rank
! complex a(d,d) (in) : The Matrix 1174 1182 ! complex a(d,d) (in) : The Matrix
! complex evala(d) (out) : eigenvalues 1175 1183 ! complex evala(d) (out) : eigenvalues
! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1176 1184 ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 1177 1185 ! integer (out) : status =0 if something went wrong
! 1178 1186 !
! interfacing Lapack routine CGEEV 1179 1187 ! interfacing Lapack routine CGEEV
implicit none 1180 1188 implicit none
integer, intent(in) :: d 1181 1189 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 1182 1190 complex, intent(in) :: a(d,d)
complex, intent(out) :: evala(d) 1183 1191 complex, intent(out) :: evala(d)
complex, intent(out) :: eveca(d,d) 1184 1192 complex, intent(out) :: eveca(d,d)
integer, intent(out), optional :: status 1185 1193 integer, intent(out), optional :: status
1186 1194
complex, allocatable :: wc_a(:,:) ! a working copy of a 1187 1195 complex, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 1188 1196 integer :: info
integer :: lwork 1189 1197 integer :: lwork
complex, allocatable :: work(:) 1190 1198 complex, allocatable :: work(:)
complex :: twork(1) 1191 1199 complex :: twork(1)
real, allocatable :: rwork(:) 1192 1200 real, allocatable :: rwork(:)
complex :: vl ! unused but necessary for the call 1193 1201 complex :: vl ! unused but necessary for the call
1194 1202
if (present(status)) status=1 1195 1203 if (present(status)) status=1
1196 1204
! making a working copy of a 1197 1205 ! making a working copy of a
allocate(wc_a(d,d)) 1198 1206 allocate(wc_a(d,d))
!call ccopy(d*d,a,1,wc_a,1) 1199 1207 !call ccopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1200 1208 wc_a(:,:)=a(:,:)
1201 1209
1202 1210
! query optimal work size 1203 1211 ! query optimal work size
call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) 1204 1212 call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
lwork=int(twork(1)) 1205 1213 lwork=int(twork(1))
allocate(work(lwork)) 1206 1214 allocate(work(lwork))
allocate(rwork(2*d)) 1207 1215 allocate(rwork(2*d))
call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) 1208 1216 call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
1209 1217
if (info /= 0) then 1210 1218 if (info /= 0) then
if (present(status)) status=0 1211 1219 if (present(status)) status=0
end if 1212 1220 end if
deallocate(rwork) 1213 1221 deallocate(rwork)
deallocate(work) 1214 1222 deallocate(work)
deallocate(wc_a) 1215 1223 deallocate(wc_a)
1216 1224
end subroutine 1217 1225 end subroutine
1218 1226
subroutine fvn_z_matev(d,a,evala,eveca,status) 1219 1227 subroutine fvn_z_matev(d,a,evala,eveca,status)
! 1220 1228 !
! integer d (in) : matrice rank 1221 1229 ! integer d (in) : matrice rank
! double complex a(d,d) (in) : The Matrix 1222 1230 ! double complex a(d,d) (in) : The Matrix
! double complex evala(d) (out) : eigenvalues 1223 1231 ! double complex evala(d) (out) : eigenvalues
! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1224 1232 ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 1225 1233 ! integer (out) : status =0 if something went wrong
! 1226 1234 !
! interfacing Lapack routine ZGEEV 1227 1235 ! interfacing Lapack routine ZGEEV
implicit none 1228 1236 implicit none
integer, intent(in) :: d 1229 1237 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 1230 1238 double complex, intent(in) :: a(d,d)
double complex, intent(out) :: evala(d) 1231 1239 double complex, intent(out) :: evala(d)
double complex, intent(out) :: eveca(d,d) 1232 1240 double complex, intent(out) :: eveca(d,d)
integer, intent(out), optional :: status 1233 1241 integer, intent(out), optional :: status
1234 1242
double complex, allocatable :: wc_a(:,:) ! a working copy of a 1235 1243 double complex, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 1236 1244 integer :: info
integer :: lwork 1237 1245 integer :: lwork
double complex, allocatable :: work(:) 1238 1246 double complex, allocatable :: work(:)
double complex :: twork(1) 1239 1247 double complex :: twork(1)
double precision, allocatable :: rwork(:) 1240 1248 double precision, allocatable :: rwork(:)
double complex :: vl ! unused but necessary for the call 1241 1249 double complex :: vl ! unused but necessary for the call
1242 1250
if (present(status)) status=1 1243 1251 if (present(status)) status=1
1244 1252
! making a working copy of a 1245 1253 ! making a working copy of a
allocate(wc_a(d,d)) 1246 1254 allocate(wc_a(d,d))
!call zcopy(d*d,a,1,wc_a,1) 1247 1255 !call zcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1248 1256 wc_a(:,:)=a(:,:)
1249 1257
! query optimal work size 1250 1258 ! query optimal work size
call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) 1251 1259 call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
lwork=int(twork(1)) 1252 1260 lwork=int(twork(1))
allocate(work(lwork)) 1253 1261 allocate(work(lwork))
allocate(rwork(2*d)) 1254 1262 allocate(rwork(2*d))
call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) 1255 1263 call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
1256 1264
if (info /= 0) then 1257 1265 if (info /= 0) then
if (present(status)) status=0 1258 1266 if (present(status)) status=0
end if 1259 1267 end if
deallocate(rwork) 1260 1268 deallocate(rwork)
deallocate(work) 1261 1269 deallocate(work)
deallocate(wc_a) 1262 1270 deallocate(wc_a)
1263 1271
end subroutine 1264 1272 end subroutine
1265 1273
1266 1274
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1267 1275 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1268 1276 !
! Least square problem 1269 1277 ! Least square problem
! 1270 1278 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1271 1279 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1272 1280 !
! 1273 1281 !
1274 1282
1275 1283
1276 1284
1277 1285
subroutine fvn_d_lspoly(np,x,y,deg,coeff,status) 1278 1286 subroutine fvn_d_lspoly(np,x,y,deg,coeff,status)
! 1279 1287 !
! Least square polynomial fitting 1280 1288 ! Least square polynomial fitting
! 1281 1289 !
! Find the coefficients of the least square polynomial of a given degree 1282 1290 ! Find the coefficients of the least square polynomial of a given degree
! for a set of coordinates. 1283 1291 ! for a set of coordinates.
! 1284 1292 !
! The degree must be lower than the number of points 1285 1293 ! The degree must be lower than the number of points
! 1286 1294 !
! np (in) : number of points 1287 1295 ! np (in) : number of points
! x(np) (in) : x data 1288 1296 ! x(np) (in) : x data
! y(np) (in) : y data 1289 1297 ! y(np) (in) : y data
! deg (in) : polynomial's degree 1290 1298 ! deg (in) : polynomial's degree
! coeff(deg+1) (out) : polynomial coefficients 1291 1299 ! coeff(deg+1) (out) : polynomial coefficients
! status (out) : =0 if a problem occurs 1292 1300 ! status (out) : =0 if a problem occurs
! 1293 1301 !
implicit none 1294 1302 implicit none
1295 1303
integer, intent(in) :: np,deg 1296 1304 integer, intent(in) :: np,deg
real(kind=8), intent(in), dimension(np) :: x,y 1297 1305 real(kind=8), intent(in), dimension(np) :: x,y
real(kind=8), intent(out), dimension(deg+1) :: coeff 1298 1306 real(kind=8), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status 1299 1307 integer, intent(out), optional :: status
1300 1308
real(kind=8), allocatable, dimension(:,:) :: mat,bmat 1301 1309 real(kind=8), allocatable, dimension(:,:) :: mat,bmat
real(kind=8),dimension(:),allocatable :: work 1302 1310 real(kind=8),dimension(:),allocatable :: work