lxsd / fvn

Download as
Browse Code ยป

Commit 68642e5d21b0cb88b82aa6c7a681e3b448197363

Authored by William Daniau
1 parent 204ded3d3d
Exists in geevx_2020

Using [SDCZ]GEEVX instead of [SDCZ]GEEV to avoid use of 

[SDCZ]GEBAL which cause NaN in some circumstances

Showing 1 changed file with 49 additions and 9 deletions Inline Diff

fvn_linear/fvn_linear.f90
Diff comments   View file @ 68642e5
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(kind=sp_kind), dimension(:,:),intent(in) :: a,b 105 105 real(kind=sp_kind), dimension(:,:),intent(in) :: a,b
real(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_s_matmul 106 106 real(kind=sp_kind), 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(kind=dp_kind), dimension(:,:),intent(in) :: a,b 111 111 real(kind=dp_kind), dimension(:,:),intent(in) :: a,b
real(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_d_matmul 112 112 real(kind=dp_kind), 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(kind=sp_kind), dimension(:,:),intent(in) :: a,b 117 117 complex(kind=sp_kind), dimension(:,:),intent(in) :: a,b
complex(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_c_matmul 118 118 complex(kind=sp_kind), 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(kind=dp_kind), dimension(:,:), intent(in) :: a,b 123 123 complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_z_matmul 124 124 complex(kind=dp_kind), 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(kind=sp_kind), dimension(:,:), intent(in) :: a,b 133 133 real(kind=sp_kind), dimension(:,:), intent(in) :: a,b
real(kind=sp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_s_tx 134 134 real(kind=sp_kind), 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(kind=dp_kind), dimension(:,:), intent(in) :: a,b 139 139 real(kind=dp_kind), dimension(:,:), intent(in) :: a,b
real(kind=dp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_d_tx 140 140 real(kind=dp_kind), 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(kind=sp_kind), dimension(:,:), intent(in) :: a,b 145 145 complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=sp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_c_tx 146 146 complex(kind=sp_kind), 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(kind=dp_kind), dimension(:,:), intent(in) :: a,b 151 151 complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=dp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_z_tx 152 152 complex(kind=dp_kind), 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(kind=sp_kind), dimension(:,:), intent(in) :: a,b 162 162 real(kind=sp_kind), dimension(:,:), intent(in) :: a,b
real(kind=sp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_s_xt 163 163 real(kind=sp_kind), 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(kind=dp_kind), dimension(:,:), intent(in) :: a,b 168 168 real(kind=dp_kind), dimension(:,:), intent(in) :: a,b
real(kind=dp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_d_xt 169 169 real(kind=dp_kind), 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(kind=sp_kind), dimension(:,:), intent(in) :: a,b 174 174 complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=sp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_c_xt 175 175 complex(kind=sp_kind), 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(kind=dp_kind), dimension(:,:), intent(in) :: a,b 180 180 complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=dp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_z_xt 181 181 complex(kind=dp_kind), 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(kind=sp_kind),dimension(:,:),intent(in) :: a 190 190 real(kind=sp_kind),dimension(:,:),intent(in) :: a
real(kind=sp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_s_transpose 191 191 real(kind=sp_kind),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(kind=dp_kind),dimension(:,:),intent(in) :: a 196 196 real(kind=dp_kind),dimension(:,:),intent(in) :: a
real(kind=dp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_d_transpose 197 197 real(kind=dp_kind),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(kind=sp_kind),dimension(:,:),intent(in) :: a 202 202 complex(kind=sp_kind),dimension(:,:),intent(in) :: a
complex(kind=sp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_c_transpose 203 203 complex(kind=sp_kind),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(kind=dp_kind),dimension(:,:),intent(in) :: a 208 208 complex(kind=dp_kind),dimension(:,:),intent(in) :: a
complex(kind=dp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_z_transpose 209 209 complex(kind=dp_kind),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 !
! It seems that there's a problem with automatic arrays with gfortran 216 216 ! It seems that there's a problem with automatic arrays with gfortran
! in some circumstances. To allow compilation with gfortran we use here a temporary array 217 217 ! in some circumstances. To allow compilation with gfortran we use here a temporary array
! for the call. Without that there's a warning at compile time and a segmentation fault 218 218 ! for the call. Without that there's a warning at compile time and a segmentation fault
! during execution. This is odd as we double memory use. 219 219 ! during execution. This is odd as we double memory use.
function fvn_op_s_matinv(a) 220 220 function fvn_op_s_matinv(a)
implicit none 221 221 implicit none
real(kind=sp_kind),dimension(:,:),intent(in) :: a 222 222 real(kind=sp_kind),dimension(:,:),intent(in) :: a
real(kind=sp_kind),dimension(size(a,1),size(a,1)) :: fvn_op_s_matinv,tmp_array 223 223 real(kind=sp_kind),dimension(size(a,1),size(a,1)) :: fvn_op_s_matinv,tmp_array
call fvn_s_matinv(size(a,1),a,tmp_array,fvn_status) 224 224 call fvn_s_matinv(size(a,1),a,tmp_array,fvn_status)
fvn_op_s_matinv=tmp_array 225 225 fvn_op_s_matinv=tmp_array
end function 226 226 end function
function fvn_op_d_matinv(a) 227 227 function fvn_op_d_matinv(a)
implicit none 228 228 implicit none
real(kind=dp_kind),dimension(:,:),intent(in) :: a 229 229 real(kind=dp_kind),dimension(:,:),intent(in) :: a
real(kind=dp_kind),dimension(size(a,1),size(a,1)) :: fvn_op_d_matinv,tmp_array 230 230 real(kind=dp_kind),dimension(size(a,1),size(a,1)) :: fvn_op_d_matinv,tmp_array
call fvn_d_matinv(size(a,1),a,tmp_array,fvn_status) 231 231 call fvn_d_matinv(size(a,1),a,tmp_array,fvn_status)
fvn_op_d_matinv=tmp_array 232 232 fvn_op_d_matinv=tmp_array
end function 233 233 end function
function fvn_op_c_matinv(a) 234 234 function fvn_op_c_matinv(a)
implicit none 235 235 implicit none
complex(kind=sp_kind),dimension(:,:),intent(in) :: a 236 236 complex(kind=sp_kind),dimension(:,:),intent(in) :: a
complex(kind=sp_kind),dimension(size(a,1),size(a,1)) :: fvn_op_c_matinv,tmp_array 237 237 complex(kind=sp_kind),dimension(size(a,1),size(a,1)) :: fvn_op_c_matinv,tmp_array
call fvn_c_matinv(size(a,1),a,tmp_array,fvn_status) 238 238 call fvn_c_matinv(size(a,1),a,tmp_array,fvn_status)
fvn_op_c_matinv=tmp_array 239 239 fvn_op_c_matinv=tmp_array
end function 240 240 end function
function fvn_op_z_matinv(a) 241 241 function fvn_op_z_matinv(a)
implicit none 242 242 implicit none
complex(kind=dp_kind),dimension(:,:),intent(in) :: a 243 243 complex(kind=dp_kind),dimension(:,:),intent(in) :: a
complex(kind=dp_kind),dimension(size(a,1),size(a,1)) :: fvn_op_z_matinv,tmp_array 244 244 complex(kind=dp_kind),dimension(size(a,1),size(a,1)) :: fvn_op_z_matinv,tmp_array
call fvn_z_matinv(size(a,1),a,tmp_array,fvn_status) 245 245 call fvn_z_matinv(size(a,1),a,tmp_array,fvn_status)
fvn_op_z_matinv=tmp_array 246 246 fvn_op_z_matinv=tmp_array
end function 247 247 end function
248 248
! 249 249 !
! .ix. 250 250 ! .ix.
! 251 251 !
function fvn_op_s_ix(a,b) 252 252 function fvn_op_s_ix(a,b)
implicit none 253 253 implicit none
real(kind=sp_kind), dimension(:,:), intent(in) :: a,b 254 254 real(kind=sp_kind), dimension(:,:), intent(in) :: a,b
real(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_s_ix 255 255 real(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_s_ix
fvn_op_s_ix=matmul(fvn_op_s_matinv(a),b) 256 256 fvn_op_s_ix=matmul(fvn_op_s_matinv(a),b)
end function 257 257 end function
function fvn_op_d_ix(a,b) 258 258 function fvn_op_d_ix(a,b)
implicit none 259 259 implicit none
real(kind=dp_kind), dimension(:,:), intent(in) :: a,b 260 260 real(kind=dp_kind), dimension(:,:), intent(in) :: a,b
real(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_d_ix 261 261 real(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_d_ix
fvn_op_d_ix=matmul(fvn_op_d_matinv(a),b) 262 262 fvn_op_d_ix=matmul(fvn_op_d_matinv(a),b)
end function 263 263 end function
function fvn_op_c_ix(a,b) 264 264 function fvn_op_c_ix(a,b)
implicit none 265 265 implicit none
complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b 266 266 complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_c_ix 267 267 complex(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_c_ix
fvn_op_c_ix=matmul(fvn_op_c_matinv(a),b) 268 268 fvn_op_c_ix=matmul(fvn_op_c_matinv(a),b)
end function 269 269 end function
function fvn_op_z_ix(a,b) 270 270 function fvn_op_z_ix(a,b)
implicit none 271 271 implicit none
complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b 272 272 complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_z_ix 273 273 complex(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_z_ix
fvn_op_z_ix=matmul(fvn_op_z_matinv(a),b) 274 274 fvn_op_z_ix=matmul(fvn_op_z_matinv(a),b)
end function 275 275 end function
276 276
! 277 277 !
! .xi. 278 278 ! .xi.
! 279 279 !
function fvn_op_s_xi(a,b) 280 280 function fvn_op_s_xi(a,b)
implicit none 281 281 implicit none
real(kind=sp_kind), dimension(:,:), intent(in) :: a,b 282 282 real(kind=sp_kind), dimension(:,:), intent(in) :: a,b
real(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_s_xi 283 283 real(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_s_xi
fvn_op_s_xi=matmul(a,fvn_op_s_matinv(b)) 284 284 fvn_op_s_xi=matmul(a,fvn_op_s_matinv(b))
end function 285 285 end function
function fvn_op_d_xi(a,b) 286 286 function fvn_op_d_xi(a,b)
implicit none 287 287 implicit none
real(kind=dp_kind), dimension(:,:), intent(in) :: a,b 288 288 real(kind=dp_kind), dimension(:,:), intent(in) :: a,b
real(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_d_xi 289 289 real(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_d_xi
fvn_op_d_xi=matmul(a,fvn_op_d_matinv(b)) 290 290 fvn_op_d_xi=matmul(a,fvn_op_d_matinv(b))
end function 291 291 end function
function fvn_op_c_xi(a,b) 292 292 function fvn_op_c_xi(a,b)
implicit none 293 293 implicit none
complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b 294 294 complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_c_xi 295 295 complex(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_c_xi
fvn_op_c_xi=matmul(a,fvn_op_c_matinv(b)) 296 296 fvn_op_c_xi=matmul(a,fvn_op_c_matinv(b))
end function 297 297 end function
function fvn_op_z_xi(a,b) 298 298 function fvn_op_z_xi(a,b)
implicit none 299 299 implicit none
complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b 300 300 complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_z_xi 301 301 complex(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_z_xi
fvn_op_z_xi=matmul(a,fvn_op_z_matinv(b)) 302 302 fvn_op_z_xi=matmul(a,fvn_op_z_matinv(b))
end function 303 303 end function
304 304
! 305 305 !
! .h. 306 306 ! .h.
! 307 307 !
function fvn_op_c_conj_transpose(a) 308 308 function fvn_op_c_conj_transpose(a)
implicit none 309 309 implicit none
complex(kind=sp_kind),dimension(:,:),intent(in) :: a 310 310 complex(kind=sp_kind),dimension(:,:),intent(in) :: a
complex(kind=sp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_c_conj_transpose 311 311 complex(kind=sp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_c_conj_transpose
fvn_op_c_conj_transpose=transpose(conjg(a)) 312 312 fvn_op_c_conj_transpose=transpose(conjg(a))
end function 313 313 end function
function fvn_op_z_conj_transpose(a) 314 314 function fvn_op_z_conj_transpose(a)
implicit none 315 315 implicit none
complex(kind=dp_kind),dimension(:,:),intent(in) :: a 316 316 complex(kind=dp_kind),dimension(:,:),intent(in) :: a
complex(kind=dp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_z_conj_transpose 317 317 complex(kind=dp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_z_conj_transpose
fvn_op_z_conj_transpose=transpose(conjg(a)) 318 318 fvn_op_z_conj_transpose=transpose(conjg(a))
end function 319 319 end function
320 320
321 321
! 322 322 !
! .hx. 323 323 ! .hx.
! 324 324 !
function fvn_op_c_hx(a,b) 325 325 function fvn_op_c_hx(a,b)
implicit none 326 326 implicit none
complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b 327 327 complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=sp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_c_hx 328 328 complex(kind=sp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_c_hx
fvn_op_c_hx=matmul(transpose(conjg(a)),b) 329 329 fvn_op_c_hx=matmul(transpose(conjg(a)),b)
end function 330 330 end function
function fvn_op_z_hx(a,b) 331 331 function fvn_op_z_hx(a,b)
implicit none 332 332 implicit none
complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b 333 333 complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=dp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_z_hx 334 334 complex(kind=dp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_z_hx
fvn_op_z_hx=matmul(transpose(conjg(a)),b) 335 335 fvn_op_z_hx=matmul(transpose(conjg(a)),b)
end function 336 336 end function
337 337
338 338
! 339 339 !
! .xh. 340 340 ! .xh.
! 341 341 !
function fvn_op_c_xh(a,b) 342 342 function fvn_op_c_xh(a,b)
implicit none 343 343 implicit none
complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b 344 344 complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=sp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_c_xh 345 345 complex(kind=sp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_c_xh
fvn_op_c_xh=matmul(a,transpose(conjg(b))) 346 346 fvn_op_c_xh=matmul(a,transpose(conjg(b)))
end function 347 347 end function
function fvn_op_z_xh(a,b) 348 348 function fvn_op_z_xh(a,b)
implicit none 349 349 implicit none
complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b 350 350 complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
complex(kind=dp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_z_xh 351 351 complex(kind=dp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_z_xh
fvn_op_z_xh=matmul(a,transpose(conjg(b))) 352 352 fvn_op_z_xh=matmul(a,transpose(conjg(b)))
end function 353 353 end function
354 354
355 355
356 356
357 357
358 358
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 359 359 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 360 360 !
! Identity Matrix 361 361 ! Identity Matrix
! 362 362 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 363 363 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function fvn_d_ident(n) 364 364 function fvn_d_ident(n)
implicit none 365 365 implicit none
integer(kind=ip_kind) :: n 366 366 integer(kind=ip_kind) :: n
real(kind=dp_kind), dimension(n,n) :: fvn_d_ident 367 367 real(kind=dp_kind), dimension(n,n) :: fvn_d_ident
368 368
real(kind=dp_kind),dimension(n*n) :: vect 369 369 real(kind=dp_kind),dimension(n*n) :: vect
integer(kind=ip_kind) :: i 370 370 integer(kind=ip_kind) :: i
371 371
vect=0._dp_kind 372 372 vect=0._dp_kind
vect(1:n*n:n+1) = 1._dp_kind 373 373 vect(1:n*n:n+1) = 1._dp_kind
fvn_d_ident=reshape(vect, shape = (/ n,n /)) 374 374 fvn_d_ident=reshape(vect, shape = (/ n,n /))
end function 375 375 end function
376 376
function fvn_s_ident(n) 377 377 function fvn_s_ident(n)
implicit none 378 378 implicit none
integer(kind=ip_kind) :: n 379 379 integer(kind=ip_kind) :: n
real(kind=sp_kind), dimension(n,n) :: fvn_s_ident 380 380 real(kind=sp_kind), dimension(n,n) :: fvn_s_ident
381 381
real(kind=sp_kind),dimension(n*n) :: vect 382 382 real(kind=sp_kind),dimension(n*n) :: vect
integer(kind=ip_kind) :: i 383 383 integer(kind=ip_kind) :: i
384 384
vect=0._sp_kind 385 385 vect=0._sp_kind
vect(1:n*n:n+1) = 1._sp_kind 386 386 vect(1:n*n:n+1) = 1._sp_kind
fvn_s_ident=reshape(vect, shape = (/ n,n /)) 387 387 fvn_s_ident=reshape(vect, shape = (/ n,n /))
end function 388 388 end function
389 389
function fvn_c_ident(n) 390 390 function fvn_c_ident(n)
implicit none 391 391 implicit none
integer(kind=ip_kind) :: n 392 392 integer(kind=ip_kind) :: n
complex(kind=sp_kind), dimension(n,n) :: fvn_c_ident 393 393 complex(kind=sp_kind), dimension(n,n) :: fvn_c_ident
394 394
complex(kind=sp_kind),dimension(n*n) :: vect 395 395 complex(kind=sp_kind),dimension(n*n) :: vect
integer(kind=ip_kind) :: i 396 396 integer(kind=ip_kind) :: i
397 397
vect=(0._sp_kind,0._sp_kind) 398 398 vect=(0._sp_kind,0._sp_kind)
vect(1:n*n:n+1) = (1._sp_kind,0._sp_kind) 399 399 vect(1:n*n:n+1) = (1._sp_kind,0._sp_kind)
fvn_c_ident=reshape(vect, shape = (/ n,n /)) 400 400 fvn_c_ident=reshape(vect, shape = (/ n,n /))
end function 401 401 end function
402 402
function fvn_z_ident(n) 403 403 function fvn_z_ident(n)
implicit none 404 404 implicit none
integer(kind=ip_kind) :: n 405 405 integer(kind=ip_kind) :: n
complex(kind=dp_kind), dimension(n,n) :: fvn_z_ident 406 406 complex(kind=dp_kind), dimension(n,n) :: fvn_z_ident
407 407
complex(kind=dp_kind),dimension(n*n) :: vect 408 408 complex(kind=dp_kind),dimension(n*n) :: vect
integer(kind=ip_kind) :: i 409 409 integer(kind=ip_kind) :: i
410 410
vect=(0._dp_kind,0._dp_kind) 411 411 vect=(0._dp_kind,0._dp_kind)
vect(1:n*n:n+1) = (1._dp_kind,0._dp_kind) 412 412 vect(1:n*n:n+1) = (1._dp_kind,0._dp_kind)
fvn_z_ident=reshape(vect, shape = (/ n,n /)) 413 413 fvn_z_ident=reshape(vect, shape = (/ n,n /))
end function 414 414 end function
415 415
416 416
417 417
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 418 418 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 419 419 !
! Matrix inversion subroutines 420 420 ! Matrix inversion subroutines
! 421 421 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 422 422 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
423 423
subroutine fvn_s_matinv(d,a,inva,status) 424 424 subroutine fvn_s_matinv(d,a,inva,status)
! 425 425 !
! Matrix inversion of a real(kind=sp_kind) matrix using BLAS and LAPACK 426 426 ! Matrix inversion of a real(kind=sp_kind) matrix using BLAS and LAPACK
! 427 427 !
! d (in) : matrix rank 428 428 ! d (in) : matrix rank
! a (in) : input matrix 429 429 ! a (in) : input matrix
! inva (out) : inversed matrix 430 430 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 431 431 ! status (ou) : =0 if something failed
! 432 432 !
implicit none 433 433 implicit none
integer(kind=ip_kind), intent(in) :: d 434 434 integer(kind=ip_kind), intent(in) :: d
real(kind=sp_kind), intent(in) :: a(d,d) 435 435 real(kind=sp_kind), intent(in) :: a(d,d)
real(kind=sp_kind), intent(out) :: inva(d,d) 436 436 real(kind=sp_kind), intent(out) :: inva(d,d)
integer(kind=ip_kind), intent(out),optional :: status 437 437 integer(kind=ip_kind), intent(out),optional :: status
438 438
integer(kind=ip_kind), allocatable :: ipiv(:) 439 439 integer(kind=ip_kind), allocatable :: ipiv(:)
real(kind=sp_kind), allocatable :: work(:) 440 440 real(kind=sp_kind), allocatable :: work(:)
real(kind=sp_kind) twork(1) 441 441 real(kind=sp_kind) twork(1)
integer(kind=ip_kind) :: info 442 442 integer(kind=ip_kind) :: info
integer(kind=ip_kind) :: lwork 443 443 integer(kind=ip_kind) :: lwork
444 444
if (present(status)) status=1 445 445 if (present(status)) status=1
446 446
allocate(ipiv(d)) 447 447 allocate(ipiv(d))
! copy a into inva using BLAS 448 448 ! copy a into inva using BLAS
!call scopy(d*d,a,1,inva,1) 449 449 !call scopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 450 450 inva(:,:)=a(:,:)
! LU factorization using LAPACK 451 451 ! LU factorization using LAPACK
call sgetrf(d,d,inva,d,ipiv,info) 452 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 453 453 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 454 454 if (info /= 0) then
if (present(status)) status=0 455 455 if (present(status)) status=0
deallocate(ipiv) 456 456 deallocate(ipiv)
return 457 457 return
end if 458 458 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 459 459 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call sgetri(d,inva,d,ipiv,twork,-1,info) 460 460 call sgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 461 461 lwork=int(twork(1))
allocate(work(lwork)) 462 462 allocate(work(lwork))
! Matrix inversion using LAPACK 463 463 ! Matrix inversion using LAPACK
call sgetri(d,inva,d,ipiv,work,lwork,info) 464 464 call sgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 465 465 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 466 466 if (info /= 0) then
if (present(status)) status=0 467 467 if (present(status)) status=0
end if 468 468 end if
deallocate(work) 469 469 deallocate(work)
deallocate(ipiv) 470 470 deallocate(ipiv)
end subroutine 471 471 end subroutine
472 472
subroutine fvn_d_matinv(d,a,inva,status) 473 473 subroutine fvn_d_matinv(d,a,inva,status)
! 474 474 !
! Matrix inversion of a real(kind=dp_kind) matrix using BLAS and LAPACK 475 475 ! Matrix inversion of a real(kind=dp_kind) matrix using BLAS and LAPACK
! 476 476 !
! d (in) : matrix rank 477 477 ! d (in) : matrix rank
! a (in) : input matrix 478 478 ! a (in) : input matrix
! inva (out) : inversed matrix 479 479 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 480 480 ! status (ou) : =0 if something failed
! 481 481 !
implicit none 482 482 implicit none
integer(kind=ip_kind), intent(in) :: d 483 483 integer(kind=ip_kind), intent(in) :: d
real(kind=dp_kind), intent(in) :: a(d,d) 484 484 real(kind=dp_kind), intent(in) :: a(d,d)
real(kind=dp_kind), intent(out) :: inva(d,d) 485 485 real(kind=dp_kind), intent(out) :: inva(d,d)
integer(kind=ip_kind), intent(out),optional :: status 486 486 integer(kind=ip_kind), intent(out),optional :: status
487 487
integer(kind=ip_kind), allocatable :: ipiv(:) 488 488 integer(kind=ip_kind), allocatable :: ipiv(:)
real(kind=dp_kind), allocatable :: work(:) 489 489 real(kind=dp_kind), allocatable :: work(:)
real(kind=dp_kind) :: twork(1) 490 490 real(kind=dp_kind) :: twork(1)
integer(kind=ip_kind) :: info 491 491 integer(kind=ip_kind) :: info
integer(kind=ip_kind) :: lwork 492 492 integer(kind=ip_kind) :: lwork
493 493
if (present(status)) status=1 494 494 if (present(status)) status=1
495 495
allocate(ipiv(d)) 496 496 allocate(ipiv(d))
! copy a into inva using BLAS 497 497 ! copy a into inva using BLAS
!call dcopy(d*d,a,1,inva,1) 498 498 !call dcopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 499 499 inva(:,:)=a(:,:)
! LU factorization using LAPACK 500 500 ! LU factorization using LAPACK
call dgetrf(d,d,inva,d,ipiv,info) 501 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 502 502 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 503 503 if (info /= 0) then
if (present(status)) status=0 504 504 if (present(status)) status=0
deallocate(ipiv) 505 505 deallocate(ipiv)
return 506 506 return
end if 507 507 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 508 508 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call dgetri(d,inva,d,ipiv,twork,-1,info) 509 509 call dgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 510 510 lwork=int(twork(1))
allocate(work(lwork)) 511 511 allocate(work(lwork))
! Matrix inversion using LAPACK 512 512 ! Matrix inversion using LAPACK
call dgetri(d,inva,d,ipiv,work,lwork,info) 513 513 call dgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 514 514 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 515 515 if (info /= 0) then
if (present(status)) status=0 516 516 if (present(status)) status=0
end if 517 517 end if
deallocate(work) 518 518 deallocate(work)
deallocate(ipiv) 519 519 deallocate(ipiv)
end subroutine 520 520 end subroutine
521 521
subroutine fvn_c_matinv(d,a,inva,status) 522 522 subroutine fvn_c_matinv(d,a,inva,status)
! 523 523 !
! Matrix inversion of a complex(kind=sp_kind) matrix using BLAS and LAPACK 524 524 ! Matrix inversion of a complex(kind=sp_kind) matrix using BLAS and LAPACK
! 525 525 !
! d (in) : matrix rank 526 526 ! d (in) : matrix rank
! a (in) : input matrix 527 527 ! a (in) : input matrix
! inva (out) : inversed matrix 528 528 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 529 529 ! status (ou) : =0 if something failed
! 530 530 !
implicit none 531 531 implicit none
integer(kind=ip_kind), intent(in) :: d 532 532 integer(kind=ip_kind), intent(in) :: d
complex(kind=sp_kind), intent(in) :: a(d,d) 533 533 complex(kind=sp_kind), intent(in) :: a(d,d)
complex(kind=sp_kind), intent(out) :: inva(d,d) 534 534 complex(kind=sp_kind), intent(out) :: inva(d,d)
integer(kind=ip_kind), intent(out),optional :: status 535 535 integer(kind=ip_kind), intent(out),optional :: status
536 536
integer(kind=ip_kind), allocatable :: ipiv(:) 537 537 integer(kind=ip_kind), allocatable :: ipiv(:)
complex(kind=sp_kind), allocatable :: work(:) 538 538 complex(kind=sp_kind), allocatable :: work(:)
complex(kind=sp_kind) :: twork(1) 539 539 complex(kind=sp_kind) :: twork(1)
integer(kind=ip_kind) :: info 540 540 integer(kind=ip_kind) :: info
integer(kind=ip_kind) :: lwork 541 541 integer(kind=ip_kind) :: lwork
542 542
if (present(status)) status=1 543 543 if (present(status)) status=1
544 544
allocate(ipiv(d)) 545 545 allocate(ipiv(d))
! copy a into inva using BLAS 546 546 ! copy a into inva using BLAS
!call ccopy(d*d,a,1,inva,1) 547 547 !call ccopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 548 548 inva(:,:)=a(:,:)
549 549
! LU factorization using LAPACK 550 550 ! LU factorization using LAPACK
call cgetrf(d,d,inva,d,ipiv,info) 551 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 552 552 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 553 553 if (info /= 0) then
if (present(status)) status=0 554 554 if (present(status)) status=0
deallocate(ipiv) 555 555 deallocate(ipiv)
return 556 556 return
end if 557 557 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 558 558 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call cgetri(d,inva,d,ipiv,twork,-1,info) 559 559 call cgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 560 560 lwork=int(twork(1))
allocate(work(lwork)) 561 561 allocate(work(lwork))
! Matrix inversion using LAPACK 562 562 ! Matrix inversion using LAPACK
call cgetri(d,inva,d,ipiv,work,lwork,info) 563 563 call cgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 564 564 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 565 565 if (info /= 0) then
if (present(status)) status=0 566 566 if (present(status)) status=0
end if 567 567 end if
deallocate(work) 568 568 deallocate(work)
deallocate(ipiv) 569 569 deallocate(ipiv)
end subroutine 570 570 end subroutine
571 571
subroutine fvn_z_matinv(d,a,inva,status) 572 572 subroutine fvn_z_matinv(d,a,inva,status)
! 573 573 !
! Matrix inversion of a complex(kind=dp_kind)(kind=sp_kind) matrix using BLAS and LAPACK 574 574 ! Matrix inversion of a complex(kind=dp_kind)(kind=sp_kind) matrix using BLAS and LAPACK
! 575 575 !
! d (in) : matrix rank 576 576 ! d (in) : matrix rank
! a (in) : input matrix 577 577 ! a (in) : input matrix
! inva (out) : inversed matrix 578 578 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 579 579 ! status (ou) : =0 if something failed
! 580 580 !
implicit none 581 581 implicit none
integer(kind=ip_kind), intent(in) :: d 582 582 integer(kind=ip_kind), intent(in) :: d
complex(kind=dp_kind), intent(in) :: a(d,d) 583 583 complex(kind=dp_kind), intent(in) :: a(d,d)
complex(kind=dp_kind), intent(out) :: inva(d,d) 584 584 complex(kind=dp_kind), intent(out) :: inva(d,d)
integer(kind=ip_kind), intent(out),optional :: status 585 585 integer(kind=ip_kind), intent(out),optional :: status
586 586
integer(kind=ip_kind), allocatable :: ipiv(:) 587 587 integer(kind=ip_kind), allocatable :: ipiv(:)
complex(kind=dp_kind), allocatable :: work(:) 588 588 complex(kind=dp_kind), allocatable :: work(:)
complex(kind=dp_kind) :: twork(1) 589 589 complex(kind=dp_kind) :: twork(1)
integer(kind=ip_kind) :: info 590 590 integer(kind=ip_kind) :: info
integer(kind=ip_kind) :: lwork 591 591 integer(kind=ip_kind) :: lwork
592 592
if (present(status)) status=1 593 593 if (present(status)) status=1
594 594
allocate(ipiv(d)) 595 595 allocate(ipiv(d))
! copy a into inva using BLAS 596 596 ! copy a into inva using BLAS
!call zcopy(d*d,a,1,inva,1) 597 597 !call zcopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 598 598 inva(:,:)=a(:,:)
599 599
! LU factorization using LAPACK 600 600 ! LU factorization using LAPACK
call zgetrf(d,d,inva,d,ipiv,info) 601 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 602 602 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 603 603 if (info /= 0) then
if (present(status)) status=0 604 604 if (present(status)) status=0
deallocate(ipiv) 605 605 deallocate(ipiv)
return 606 606 return
end if 607 607 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 608 608 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call zgetri(d,inva,d,ipiv,twork,-1,info) 609 609 call zgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 610 610 lwork=int(twork(1))
allocate(work(lwork)) 611 611 allocate(work(lwork))
! Matrix inversion using LAPACK 612 612 ! Matrix inversion using LAPACK
call zgetri(d,inva,d,ipiv,work,lwork,info) 613 613 call zgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 614 614 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 615 615 if (info /= 0) then
if (present(status)) status=0 616 616 if (present(status)) status=0
end if 617 617 end if
deallocate(work) 618 618 deallocate(work)
deallocate(ipiv) 619 619 deallocate(ipiv)
end subroutine 620 620 end subroutine
621 621
622 622
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 623 623 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 624 624 !
! Determinants 625 625 ! Determinants
! 626 626 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 627 627 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function fvn_s_det(d,a,status) 628 628 function fvn_s_det(d,a,status)
! 629 629 !
! Evaluate the determinant of a square matrix using lapack LU factorization 630 630 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 631 631 !
! d (in) : matrix rank 632 632 ! d (in) : matrix rank
! a (in) : The Matrix 633 633 ! a (in) : The Matrix
! status (out) : =0 if LU factorization failed 634 634 ! status (out) : =0 if LU factorization failed
! 635 635 !
implicit none 636 636 implicit none
integer(kind=ip_kind), intent(in) :: d 637 637 integer(kind=ip_kind), intent(in) :: d
real(kind=sp_kind), intent(in) :: a(d,d) 638 638 real(kind=sp_kind), intent(in) :: a(d,d)
integer(kind=ip_kind), intent(out), optional :: status 639 639 integer(kind=ip_kind), intent(out), optional :: status
real(kind=sp_kind) :: fvn_s_det 640 640 real(kind=sp_kind) :: fvn_s_det
641 641
real(kind=sp_kind), allocatable :: wc_a(:,:) 642 642 real(kind=sp_kind), allocatable :: wc_a(:,:)
integer(kind=ip_kind), allocatable :: ipiv(:) 643 643 integer(kind=ip_kind), allocatable :: ipiv(:)
integer(kind=ip_kind) :: info,i 644 644 integer(kind=ip_kind) :: info,i
645 645
if (present(status)) status=1 646 646 if (present(status)) status=1
allocate(wc_a(d,d)) 647 647 allocate(wc_a(d,d))
allocate(ipiv(d)) 648 648 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 649 649 wc_a(:,:)=a(:,:)
call sgetrf(d,d,wc_a,d,ipiv,info) 650 650 call sgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 651 651 if (info/= 0) then
if (present(status)) status=0 652 652 if (present(status)) status=0
fvn_s_det=0.e0 653 653 fvn_s_det=0.e0
deallocate(ipiv) 654 654 deallocate(ipiv)
deallocate(wc_a) 655 655 deallocate(wc_a)
return 656 656 return
end if 657 657 end if
fvn_s_det=1.e0 658 658 fvn_s_det=1.e0
do i=1,d 659 659 do i=1,d
if (ipiv(i)==i) then 660 660 if (ipiv(i)==i) then
fvn_s_det=fvn_s_det*wc_a(i,i) 661 661 fvn_s_det=fvn_s_det*wc_a(i,i)
else 662 662 else
fvn_s_det=-fvn_s_det*wc_a(i,i) 663 663 fvn_s_det=-fvn_s_det*wc_a(i,i)
end if 664 664 end if
end do 665 665 end do
deallocate(ipiv) 666 666 deallocate(ipiv)
deallocate(wc_a) 667 667 deallocate(wc_a)
668 668
end function 669 669 end function
670 670
function fvn_d_det(d,a,status) 671 671 function fvn_d_det(d,a,status)
! 672 672 !
! Evaluate the determinant of a square matrix using lapack LU factorization 673 673 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 674 674 !
! d (in) : matrix rank 675 675 ! d (in) : matrix rank
! a (in) : The Matrix 676 676 ! a (in) : The Matrix
! status (out) : =0 if LU factorization failed 677 677 ! status (out) : =0 if LU factorization failed
! 678 678 !
implicit none 679 679 implicit none
integer(kind=ip_kind), intent(in) :: d 680 680 integer(kind=ip_kind), intent(in) :: d
real(kind=dp_kind), intent(in) :: a(d,d) 681 681 real(kind=dp_kind), intent(in) :: a(d,d)
integer(kind=ip_kind), intent(out), optional :: status 682 682 integer(kind=ip_kind), intent(out), optional :: status
real(kind=dp_kind) :: fvn_d_det 683 683 real(kind=dp_kind) :: fvn_d_det
684 684
real(kind=dp_kind), allocatable :: wc_a(:,:) 685 685 real(kind=dp_kind), allocatable :: wc_a(:,:)
integer(kind=ip_kind), allocatable :: ipiv(:) 686 686 integer(kind=ip_kind), allocatable :: ipiv(:)
integer(kind=ip_kind) :: info,i 687 687 integer(kind=ip_kind) :: info,i
688 688
if (present(status)) status=1 689 689 if (present(status)) status=1
allocate(wc_a(d,d)) 690 690 allocate(wc_a(d,d))
allocate(ipiv(d)) 691 691 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 692 692 wc_a(:,:)=a(:,:)
call dgetrf(d,d,wc_a,d,ipiv,info) 693 693 call dgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 694 694 if (info/= 0) then
if (present(status)) status=0 695 695 if (present(status)) status=0
fvn_d_det=0.d0 696 696 fvn_d_det=0.d0
deallocate(ipiv) 697 697 deallocate(ipiv)
deallocate(wc_a) 698 698 deallocate(wc_a)
return 699 699 return
end if 700 700 end if
fvn_d_det=1.d0 701 701 fvn_d_det=1.d0
do i=1,d 702 702 do i=1,d
if (ipiv(i)==i) then 703 703 if (ipiv(i)==i) then
fvn_d_det=fvn_d_det*wc_a(i,i) 704 704 fvn_d_det=fvn_d_det*wc_a(i,i)
else 705 705 else
fvn_d_det=-fvn_d_det*wc_a(i,i) 706 706 fvn_d_det=-fvn_d_det*wc_a(i,i)
end if 707 707 end if
end do 708 708 end do
deallocate(ipiv) 709 709 deallocate(ipiv)
deallocate(wc_a) 710 710 deallocate(wc_a)
711 711
end function 712 712 end function
713 713
function fvn_c_det(d,a,status) ! 714 714 function fvn_c_det(d,a,status) !
! Evaluate the determinant of a square matrix using lapack LU factorization 715 715 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 716 716 !
! d (in) : matrix rank 717 717 ! d (in) : matrix rank
! a (in) : The Matrix 718 718 ! a (in) : The Matrix
! status (out) : =0 if LU factorization failed 719 719 ! status (out) : =0 if LU factorization failed
! 720 720 !
implicit none 721 721 implicit none
integer(kind=ip_kind), intent(in) :: d 722 722 integer(kind=ip_kind), intent(in) :: d
complex(kind=sp_kind), intent(in) :: a(d,d) 723 723 complex(kind=sp_kind), intent(in) :: a(d,d)
integer(kind=ip_kind), intent(out), optional :: status 724 724 integer(kind=ip_kind), intent(out), optional :: status
complex(kind=sp_kind) :: fvn_c_det 725 725 complex(kind=sp_kind) :: fvn_c_det
726 726
complex(kind=sp_kind), allocatable :: wc_a(:,:) 727 727 complex(kind=sp_kind), allocatable :: wc_a(:,:)
integer(kind=ip_kind), allocatable :: ipiv(:) 728 728 integer(kind=ip_kind), allocatable :: ipiv(:)
integer(kind=ip_kind) :: info,i 729 729 integer(kind=ip_kind) :: info,i
730 730
if (present(status)) status=1 731 731 if (present(status)) status=1
allocate(wc_a(d,d)) 732 732 allocate(wc_a(d,d))
allocate(ipiv(d)) 733 733 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 734 734 wc_a(:,:)=a(:,:)
call cgetrf(d,d,wc_a,d,ipiv,info) 735 735 call cgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 736 736 if (info/= 0) then
if (present(status)) status=0 737 737 if (present(status)) status=0
fvn_c_det=(0.e0,0.e0) 738 738 fvn_c_det=(0.e0,0.e0)
deallocate(ipiv) 739 739 deallocate(ipiv)
deallocate(wc_a) 740 740 deallocate(wc_a)
return 741 741 return
end if 742 742 end if
fvn_c_det=(1.e0,0.e0) 743 743 fvn_c_det=(1.e0,0.e0)
do i=1,d 744 744 do i=1,d
if (ipiv(i)==i) then 745 745 if (ipiv(i)==i) then
fvn_c_det=fvn_c_det*wc_a(i,i) 746 746 fvn_c_det=fvn_c_det*wc_a(i,i)
else 747 747 else
fvn_c_det=-fvn_c_det*wc_a(i,i) 748 748 fvn_c_det=-fvn_c_det*wc_a(i,i)
end if 749 749 end if
end do 750 750 end do
deallocate(ipiv) 751 751 deallocate(ipiv)
deallocate(wc_a) 752 752 deallocate(wc_a)
753 753
end function 754 754 end function
755 755
function fvn_z_det(d,a,status) 756 756 function fvn_z_det(d,a,status)
! 757 757 !
! Evaluate the determinant of a square matrix using lapack LU factorization 758 758 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 759 759 !
! d (in) : matrix rank 760 760 ! d (in) : matrix rank
! a (in) : The Matrix 761 761 ! a (in) : The Matrix
! det (out) : determinant 762 762 ! det (out) : determinant
! status (out) : =0 if LU factorization failed 763 763 ! status (out) : =0 if LU factorization failed
! 764 764 !
implicit none 765 765 implicit none
integer(kind=ip_kind), intent(in) :: d 766 766 integer(kind=ip_kind), intent(in) :: d
complex(kind=dp_kind), intent(in) :: a(d,d) 767 767 complex(kind=dp_kind), intent(in) :: a(d,d)
integer(kind=ip_kind), intent(out), optional :: status 768 768 integer(kind=ip_kind), intent(out), optional :: status
complex(kind=dp_kind) :: fvn_z_det 769 769 complex(kind=dp_kind) :: fvn_z_det
770 770
complex(kind=dp_kind), allocatable :: wc_a(:,:) 771 771 complex(kind=dp_kind), allocatable :: wc_a(:,:)
integer(kind=ip_kind), allocatable :: ipiv(:) 772 772 integer(kind=ip_kind), allocatable :: ipiv(:)
integer(kind=ip_kind) :: info,i 773 773 integer(kind=ip_kind) :: info,i
774 774
if (present(status)) status=1 775 775 if (present(status)) status=1
allocate(wc_a(d,d)) 776 776 allocate(wc_a(d,d))
allocate(ipiv(d)) 777 777 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 778 778 wc_a(:,:)=a(:,:)
call zgetrf(d,d,wc_a,d,ipiv,info) 779 779 call zgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 780 780 if (info/= 0) then
if (present(status)) status=0 781 781 if (present(status)) status=0
fvn_z_det=(0.d0,0.d0) 782 782 fvn_z_det=(0.d0,0.d0)
deallocate(ipiv) 783 783 deallocate(ipiv)
deallocate(wc_a) 784 784 deallocate(wc_a)
return 785 785 return
end if 786 786 end if
fvn_z_det=(1.d0,0.d0) 787 787 fvn_z_det=(1.d0,0.d0)
do i=1,d 788 788 do i=1,d
if (ipiv(i)==i) then 789 789 if (ipiv(i)==i) then
fvn_z_det=fvn_z_det*wc_a(i,i) 790 790 fvn_z_det=fvn_z_det*wc_a(i,i)
else 791 791 else
fvn_z_det=-fvn_z_det*wc_a(i,i) 792 792 fvn_z_det=-fvn_z_det*wc_a(i,i)
end if 793 793 end if
end do 794 794 end do
deallocate(ipiv) 795 795 deallocate(ipiv)
deallocate(wc_a) 796 796 deallocate(wc_a)
797 797
end function 798 798 end function
799 799
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 800 800 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 801 801 !
! Condition test 802 802 ! Condition test
! 803 803 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 804 804 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1-norm 805 805 ! 1-norm
! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm 806 806 ! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm
! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond 807 807 ! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond
! 808 808 !
subroutine fvn_s_matcon(d,a,rcond,status) 809 809 subroutine fvn_s_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 810 810 ! Matrix condition (reciprocal of condition number)
! 811 811 !
! d (in) : matrix rank 812 812 ! d (in) : matrix rank
! a (in) : The Matrix 813 813 ! a (in) : The Matrix
! rcond (out) : guess what 814 814 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 815 815 ! status (out) : =0 if something went wrong
! 816 816 !
implicit none 817 817 implicit none
integer(kind=ip_kind), intent(in) :: d 818 818 integer(kind=ip_kind), intent(in) :: d
real(kind=sp_kind), intent(in) :: a(d,d) 819 819 real(kind=sp_kind), intent(in) :: a(d,d)
real(kind=sp_kind), intent(out) :: rcond 820 820 real(kind=sp_kind), intent(out) :: rcond
integer(kind=ip_kind), intent(out), optional :: status 821 821 integer(kind=ip_kind), intent(out), optional :: status
822 822
real(kind=sp_kind), allocatable :: work(:) 823 823 real(kind=sp_kind), allocatable :: work(:)
integer(kind=ip_kind), allocatable :: iwork(:) 824 824 integer(kind=ip_kind), allocatable :: iwork(:)
real(kind=sp_kind) :: anorm 825 825 real(kind=sp_kind) :: anorm
real(kind=sp_kind), allocatable :: wc_a(:,:) ! working copy of a 826 826 real(kind=sp_kind), allocatable :: wc_a(:,:) ! working copy of a
integer(kind=ip_kind) :: info 827 827 integer(kind=ip_kind) :: info
integer(kind=ip_kind), allocatable :: ipiv(:) 828 828 integer(kind=ip_kind), allocatable :: ipiv(:)
829 829
real(kind=sp_kind), external :: slange 830 830 real(kind=sp_kind), external :: slange
831 831
832 832
if (present(status)) status=1 833 833 if (present(status)) status=1
834 834
anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm 835 835 anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
836 836
allocate(wc_a(d,d)) 837 837 allocate(wc_a(d,d))
!call scopy(d*d,a,1,wc_a,1) 838 838 !call scopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 839 839 wc_a(:,:)=a(:,:)
840 840
allocate(ipiv(d)) 841 841 allocate(ipiv(d))
call sgetrf(d,d,wc_a,d,ipiv,info) 842 842 call sgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 843 843 if (info /= 0) then
if (present(status)) status=0 844 844 if (present(status)) status=0
deallocate(ipiv) 845 845 deallocate(ipiv)
deallocate(wc_a) 846 846 deallocate(wc_a)
return 847 847 return
end if 848 848 end if
allocate(work(4*d)) 849 849 allocate(work(4*d))
allocate(iwork(d)) 850 850 allocate(iwork(d))
call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info) 851 851 call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
if (info /= 0) then 852 852 if (info /= 0) then
if (present(status)) status=0 853 853 if (present(status)) status=0
end if 854 854 end if
deallocate(iwork) 855 855 deallocate(iwork)
deallocate(work) 856 856 deallocate(work)
deallocate(ipiv) 857 857 deallocate(ipiv)
deallocate(wc_a) 858 858 deallocate(wc_a)
859 859
end subroutine 860 860 end subroutine
861 861
subroutine fvn_d_matcon(d,a,rcond,status) 862 862 subroutine fvn_d_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 863 863 ! Matrix condition (reciprocal of condition number)
! 864 864 !
! d (in) : matrix rank 865 865 ! d (in) : matrix rank
! a (in) : The Matrix 866 866 ! a (in) : The Matrix
! rcond (out) : guess what 867 867 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 868 868 ! status (out) : =0 if something went wrong
! 869 869 !
implicit none 870 870 implicit none
integer(kind=ip_kind), intent(in) :: d 871 871 integer(kind=ip_kind), intent(in) :: d
real(kind=dp_kind), intent(in) :: a(d,d) 872 872 real(kind=dp_kind), intent(in) :: a(d,d)
real(kind=dp_kind), intent(out) :: rcond 873 873 real(kind=dp_kind), intent(out) :: rcond
integer(kind=ip_kind), intent(out), optional :: status 874 874 integer(kind=ip_kind), intent(out), optional :: status
875 875
real(kind=dp_kind), allocatable :: work(:) 876 876 real(kind=dp_kind), allocatable :: work(:)
integer(kind=ip_kind), allocatable :: iwork(:) 877 877 integer(kind=ip_kind), allocatable :: iwork(:)
real(kind=dp_kind) :: anorm 878 878 real(kind=dp_kind) :: anorm
real(kind=dp_kind), allocatable :: wc_a(:,:) ! working copy of a 879 879 real(kind=dp_kind), allocatable :: wc_a(:,:) ! working copy of a
integer(kind=ip_kind) :: info 880 880 integer(kind=ip_kind) :: info
integer(kind=ip_kind), allocatable :: ipiv(:) 881 881 integer(kind=ip_kind), allocatable :: ipiv(:)
882 882
real(kind=dp_kind), external :: dlange 883 883 real(kind=dp_kind), external :: dlange
884 884
885 885
if (present(status)) status=1 886 886 if (present(status)) status=1
887 887
anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm 888 888 anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
889 889
allocate(wc_a(d,d)) 890 890 allocate(wc_a(d,d))
!call dcopy(d*d,a,1,wc_a,1) 891 891 !call dcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 892 892 wc_a(:,:)=a(:,:)
893 893
allocate(ipiv(d)) 894 894 allocate(ipiv(d))
call dgetrf(d,d,wc_a,d,ipiv,info) 895 895 call dgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 896 896 if (info /= 0) then
if (present(status)) status=0 897 897 if (present(status)) status=0
deallocate(ipiv) 898 898 deallocate(ipiv)
deallocate(wc_a) 899 899 deallocate(wc_a)
return 900 900 return
end if 901 901 end if
902 902
allocate(work(4*d)) 903 903 allocate(work(4*d))
allocate(iwork(d)) 904 904 allocate(iwork(d))
call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info) 905 905 call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
if (info /= 0) then 906 906 if (info /= 0) then
if (present(status)) status=0 907 907 if (present(status)) status=0
end if 908 908 end if
deallocate(iwork) 909 909 deallocate(iwork)
deallocate(work) 910 910 deallocate(work)
deallocate(ipiv) 911 911 deallocate(ipiv)
deallocate(wc_a) 912 912 deallocate(wc_a)
913 913
end subroutine 914 914 end subroutine
915 915
subroutine fvn_c_matcon(d,a,rcond,status) 916 916 subroutine fvn_c_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 917 917 ! Matrix condition (reciprocal of condition number)
! 918 918 !
! d (in) : matrix rank 919 919 ! d (in) : matrix rank
! a (in) : The Matrix 920 920 ! a (in) : The Matrix
! rcond (out) : guess what 921 921 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 922 922 ! status (out) : =0 if something went wrong
! 923 923 !
implicit none 924 924 implicit none
integer(kind=ip_kind), intent(in) :: d 925 925 integer(kind=ip_kind), intent(in) :: d
complex(kind=sp_kind), intent(in) :: a(d,d) 926 926 complex(kind=sp_kind), intent(in) :: a(d,d)
real(kind=sp_kind), intent(out) :: rcond 927 927 real(kind=sp_kind), intent(out) :: rcond
integer(kind=ip_kind), intent(out), optional :: status 928 928 integer(kind=ip_kind), intent(out), optional :: status
929 929
real(kind=sp_kind), allocatable :: rwork(:) 930 930 real(kind=sp_kind), allocatable :: rwork(:)
complex(kind=sp_kind), allocatable :: work(:) 931 931 complex(kind=sp_kind), allocatable :: work(:)
integer(kind=ip_kind), allocatable :: iwork(:) 932 932 integer(kind=ip_kind), allocatable :: iwork(:)
real(kind=sp_kind) :: anorm 933 933 real(kind=sp_kind) :: anorm
complex(kind=sp_kind), allocatable :: wc_a(:,:) ! working copy of a 934 934 complex(kind=sp_kind), allocatable :: wc_a(:,:) ! working copy of a
integer(kind=ip_kind) :: info 935 935 integer(kind=ip_kind) :: info
integer(kind=ip_kind), allocatable :: ipiv(:) 936 936 integer(kind=ip_kind), allocatable :: ipiv(:)
937 937
real(kind=sp_kind), external :: clange 938 938 real(kind=sp_kind), external :: clange
939 939
940 940
if (present(status)) status=1 941 941 if (present(status)) status=1
942 942
anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm 943 943 anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
944 944
allocate(wc_a(d,d)) 945 945 allocate(wc_a(d,d))
!call ccopy(d*d,a,1,wc_a,1) 946 946 !call ccopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 947 947 wc_a(:,:)=a(:,:)
948 948
allocate(ipiv(d)) 949 949 allocate(ipiv(d))
call cgetrf(d,d,wc_a,d,ipiv,info) 950 950 call cgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 951 951 if (info /= 0) then
if (present(status)) status=0 952 952 if (present(status)) status=0
deallocate(ipiv) 953 953 deallocate(ipiv)
deallocate(wc_a) 954 954 deallocate(wc_a)
return 955 955 return
end if 956 956 end if
allocate(work(2*d)) 957 957 allocate(work(2*d))
allocate(rwork(2*d)) 958 958 allocate(rwork(2*d))
call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info) 959 959 call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
if (info /= 0) then 960 960 if (info /= 0) then
if (present(status)) status=0 961 961 if (present(status)) status=0
end if 962 962 end if
deallocate(rwork) 963 963 deallocate(rwork)
deallocate(work) 964 964 deallocate(work)
deallocate(ipiv) 965 965 deallocate(ipiv)
deallocate(wc_a) 966 966 deallocate(wc_a)
end subroutine 967 967 end subroutine
968 968
subroutine fvn_z_matcon(d,a,rcond,status) 969 969 subroutine fvn_z_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 970 970 ! Matrix condition (reciprocal of condition number)
! 971 971 !
! d (in) : matrix rank 972 972 ! d (in) : matrix rank
! a (in) : The Matrix 973 973 ! a (in) : The Matrix
! rcond (out) : guess what 974 974 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 975 975 ! status (out) : =0 if something went wrong
! 976 976 !
implicit none 977 977 implicit none
integer(kind=ip_kind), intent(in) :: d 978 978 integer(kind=ip_kind), intent(in) :: d
complex(kind=dp_kind), intent(in) :: a(d,d) 979 979 complex(kind=dp_kind), intent(in) :: a(d,d)
real(kind=dp_kind), intent(out) :: rcond 980 980 real(kind=dp_kind), intent(out) :: rcond
integer(kind=ip_kind), intent(out), optional :: status 981 981 integer(kind=ip_kind), intent(out), optional :: status
982 982
complex(kind=dp_kind), allocatable :: work(:) 983 983 complex(kind=dp_kind), allocatable :: work(:)
real(kind=dp_kind), allocatable :: rwork(:) 984 984 real(kind=dp_kind), allocatable :: rwork(:)
real(kind=dp_kind) :: anorm 985 985 real(kind=dp_kind) :: anorm
complex(kind=dp_kind), allocatable :: wc_a(:,:) ! working copy of a 986 986 complex(kind=dp_kind), allocatable :: wc_a(:,:) ! working copy of a
integer(kind=ip_kind) :: info 987 987 integer(kind=ip_kind) :: info
integer(kind=ip_kind), allocatable :: ipiv(:) 988 988 integer(kind=ip_kind), allocatable :: ipiv(:)
989 989
real(kind=dp_kind), external :: zlange 990 990 real(kind=dp_kind), external :: zlange
991 991
992 992
if (present(status)) status=1 993 993 if (present(status)) status=1
994 994
anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm 995 995 anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
996 996
allocate(wc_a(d,d)) 997 997 allocate(wc_a(d,d))
!call zcopy(d*d,a,1,wc_a,1) 998 998 !call zcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 999 999 wc_a(:,:)=a(:,:)
1000 1000
allocate(ipiv(d)) 1001 1001 allocate(ipiv(d))
call zgetrf(d,d,wc_a,d,ipiv,info) 1002 1002 call zgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 1003 1003 if (info /= 0) then
if (present(status)) status=0 1004 1004 if (present(status)) status=0
deallocate(ipiv) 1005 1005 deallocate(ipiv)
deallocate(wc_a) 1006 1006 deallocate(wc_a)
return 1007 1007 return
end if 1008 1008 end if
1009 1009
allocate(work(2*d)) 1010 1010 allocate(work(2*d))
allocate(rwork(2*d)) 1011 1011 allocate(rwork(2*d))
call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info) 1012 1012 call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
if (info /= 0) then 1013 1013 if (info /= 0) then
if (present(status)) status=0 1014 1014 if (present(status)) status=0
end if 1015 1015 end if
deallocate(rwork) 1016 1016 deallocate(rwork)
deallocate(work) 1017 1017 deallocate(work)
deallocate(ipiv) 1018 1018 deallocate(ipiv)
deallocate(wc_a) 1019 1019 deallocate(wc_a)
end subroutine 1020 1020 end subroutine
1021 1021
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1022 1022 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1023 1023 !
! Valeurs propres/ Vecteurs propre 1024 1024 ! Valeurs propres/ Vecteurs propre
! 1025 1025 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1026 1026 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1027 1027
! August 2009 1028 1028 ! August 2009
! William Daniau 1029 1029 ! William Daniau
! Adding sorting of eigenvalues and vectors 1030 1030 ! Adding sorting of eigenvalues and vectors
subroutine fvn_z_sort_eigen(d,v,vec) 1031 1031 subroutine fvn_z_sort_eigen(d,v,vec)
! this routine takes in input : 1032 1032 ! this routine takes in input :
! v : a vector containing unsorted eigenvalues 1033 1033 ! v : a vector containing unsorted eigenvalues
! vec : a matrix where vec(:,j) is the eigenvector corresponding to v(j) 1034 1034 ! vec : a matrix where vec(:,j) is the eigenvector corresponding to v(j)
! 1035 1035 !
! At the end of subroutine the eigenvalues are sorted by decreasing order of 1036 1036 ! At the end of subroutine the eigenvalues are sorted by decreasing order of
! modulus so as the vectors. 1037 1037 ! modulus so as the vectors.
implicit none 1038 1038 implicit none
integer(kind=ip_kind) :: d,i,j 1039 1039 integer(kind=ip_kind) :: d,i,j
complex(kind=dp_kind),dimension(d) :: v 1040 1040 complex(kind=dp_kind),dimension(d) :: v
complex(kind=dp_kind),dimension(d,d) :: vec 1041 1041 complex(kind=dp_kind),dimension(d,d) :: vec
complex(kind=dp_kind) :: cur_v 1042 1042 complex(kind=dp_kind) :: cur_v
complex(kind=dp_kind), dimension(d) :: cur_vec 1043 1043 complex(kind=dp_kind), dimension(d) :: cur_vec
1044 1044
do i=2,d 1045 1045 do i=2,d
cur_v=v(i) 1046 1046 cur_v=v(i)
cur_vec=vec(:,i) 1047 1047 cur_vec=vec(:,i)
j=i-1 1048 1048 j=i-1
do while( (j>=1) .and. (abs(v(j)) < abs(cur_v)) ) 1049 1049 do while( (j>=1) .and. (abs(v(j)) < abs(cur_v)) )
v(j+1)=v(j) 1050 1050 v(j+1)=v(j)
vec(:,j+1)=vec(:,j) 1051 1051 vec(:,j+1)=vec(:,j)
j=j-1 1052 1052 j=j-1
end do 1053 1053 end do
v(j+1)=cur_v 1054 1054 v(j+1)=cur_v
vec(:,j+1)=cur_vec 1055 1055 vec(:,j+1)=cur_vec
end do 1056 1056 end do
1057 1057
end subroutine 1058 1058 end subroutine
1059 1059
subroutine fvn_c_sort_eigen(d,v,vec) 1060 1060 subroutine fvn_c_sort_eigen(d,v,vec)
! this routine takes in input : 1061 1061 ! this routine takes in input :
! v : a vector containing unsorted eigenvalues 1062 1062 ! v : a vector containing unsorted eigenvalues
! vec : a matrix where vec(:,j) is the eigenvector corresponding to v(j) 1063 1063 ! vec : a matrix where vec(:,j) is the eigenvector corresponding to v(j)
! 1064 1064 !
! At the end of subroutine the eigenvalues are sorted by decreasing order of 1065 1065 ! At the end of subroutine the eigenvalues are sorted by decreasing order of
! modulus so as the vectors. 1066 1066 ! modulus so as the vectors.
implicit none 1067 1067 implicit none
integer(kind=ip_kind) :: d,i,j 1068 1068 integer(kind=ip_kind) :: d,i,j
complex(kind=sp_kind),dimension(d) :: v 1069 1069 complex(kind=sp_kind),dimension(d) :: v
complex(kind=sp_kind),dimension(d,d) :: vec 1070 1070 complex(kind=sp_kind),dimension(d,d) :: vec
complex(kind=sp_kind) :: cur_v 1071 1071 complex(kind=sp_kind) :: cur_v
complex(kind=sp_kind), dimension(d) :: cur_vec 1072 1072 complex(kind=sp_kind), dimension(d) :: cur_vec
1073 1073
do i=2,d 1074 1074 do i=2,d
cur_v=v(i) 1075 1075 cur_v=v(i)
cur_vec=vec(:,i) 1076 1076 cur_vec=vec(:,i)
j=i-1 1077 1077 j=i-1
do while( (j>=1) .and. (abs(v(j)) < abs(cur_v)) ) 1078 1078 do while( (j>=1) .and. (abs(v(j)) < abs(cur_v)) )
v(j+1)=v(j) 1079 1079 v(j+1)=v(j)
vec(:,j+1)=vec(:,j) 1080 1080 vec(:,j+1)=vec(:,j)
j=j-1 1081 1081 j=j-1
end do 1082 1082 end do
v(j+1)=cur_v 1083 1083 v(j+1)=cur_v
vec(:,j+1)=cur_vec 1084 1084 vec(:,j+1)=cur_vec
end do 1085 1085 end do
1086 1086
end subroutine 1087 1087 end subroutine
1088 1088
subroutine fvn_s_matev(d,a,evala,eveca,status,sortval) 1089 1089 subroutine fvn_s_matev(d,a,evala,eveca,status,sortval)
! 1090 1090 !
! integer(kind=ip_kind) d (in) : matrice rank 1091 1091 ! integer(kind=ip_kind) d (in) : matrice rank
! real(kind=sp_kind) a(d,d) (in) : The Matrix 1092 1092 ! real(kind=sp_kind) a(d,d) (in) : The Matrix
! complex(kind=sp_kind) evala(d) (out) : eigenvalues 1093 1093 ! complex(kind=sp_kind) evala(d) (out) : eigenvalues
! complex(kind=sp_kind) eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1094 1094 ! complex(kind=sp_kind) eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer(kind=ip_kind) (out) : status =0 if something went wrong 1095 1095 ! integer(kind=ip_kind) (out) : status =0 if something went wrong
! 1096 1096 !
! interfacing Lapack routine SGEEV 1097 1097 ! interfacing Lapack routine SGEEV
1098 ! 2020 11 03 : SGEEVX with BALANC='N'
implicit none 1098 1099 implicit none
integer(kind=ip_kind), intent(in) :: d 1099 1100 integer(kind=ip_kind), intent(in) :: d
real(kind=sp_kind), intent(in) :: a(d,d) 1100 1101 real(kind=sp_kind), intent(in) :: a(d,d)
complex(kind=sp_kind), intent(out) :: evala(d) 1101 1102 complex(kind=sp_kind), intent(out) :: evala(d)
complex(kind=sp_kind), intent(out) :: eveca(d,d) 1102 1103 complex(kind=sp_kind), intent(out) :: eveca(d,d)
integer(kind=ip_kind), intent(out), optional :: status 1103 1104 integer(kind=ip_kind), intent(out), optional :: status
logical(kind=ip_kind), intent(in), optional :: sortval 1104 1105 logical(kind=ip_kind), intent(in), optional :: sortval
1105 1106
real(kind=sp_kind), allocatable :: wc_a(:,:) ! a working copy of a 1106 1107 real(kind=sp_kind), allocatable :: wc_a(:,:) ! a working copy of a
integer(kind=ip_kind) :: info 1107 1108 integer(kind=ip_kind) :: info
integer(kind=ip_kind) :: lwork 1108 1109 integer(kind=ip_kind) :: lwork
real(kind=sp_kind), allocatable :: wr(:),wi(:) 1109 1110 real(kind=sp_kind), allocatable :: wr(:),wi(:)
real(kind=sp_kind) :: vl ! unused but necessary for the call 1110 1111 real(kind=sp_kind) :: vl ! unused but necessary for the call
real(kind=sp_kind), allocatable :: vr(:,:) 1111 1112 real(kind=sp_kind), allocatable :: vr(:,:)
real(kind=sp_kind), allocatable :: work(:) 1112 1113 real(kind=sp_kind), allocatable :: work(:)
real(kind=sp_kind) :: twork(1) 1113 1114 real(kind=sp_kind) :: twork(1)
integer(kind=ip_kind) i 1114 1115 integer(kind=ip_kind) i
integer(kind=ip_kind) j 1115 1116 integer(kind=ip_kind) j
1116 1117
1118 integer(kind=ip_kind) :: ilo,ihi, iwork
1119 real(kind=sp_kind), allocatable, dimension(:) :: scal, rconde, rcondv
1120 real(kind=sp_kind) :: abnrm
1121
if (present(status)) status=1 1117 1122 if (present(status)) status=1
1118 1123
! making a working copy of a 1119 1124 ! making a working copy of a
allocate(wc_a(d,d)) 1120 1125 allocate(wc_a(d,d))
!call scopy(d*d,a,1,wc_a,1) 1121 1126 !call scopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1122 1127 wc_a(:,:)=a(:,:)
1123 1128
allocate(wr(d)) 1124 1129 allocate(wr(d))
allocate(wi(d)) 1125 1130 allocate(wi(d))
allocate(vr(d,d)) 1126 1131 allocate(vr(d,d))
1132 allocate(scal(d),rconde(d),rcondv(d))
! query optimal work size 1127 1133 ! query optimal work size
call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) 1128 1134 ! call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
1135 call sgeevx('N','N','V','N',d,wc_a,d,wr,wi,vl,1,vr,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,twork,-1,IWORK,info)
lwork=int(twork(1)) 1129 1136 lwork=int(twork(1))
allocate(work(lwork)) 1130 1137 allocate(work(lwork))
call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) 1131 1138 ! call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
1139 call sgeevx('N','N','V','N',d,wc_a,d,wr,wi,vl,1,vr,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,work,lwork,IWORK,info)
1132 1140
1141
if (info /= 0) then 1133 1142 if (info /= 0) then
if (present(status)) status=0 1134 1143 if (present(status)) status=0
deallocate(work) 1135 1144 deallocate(work)
deallocate(vr) 1136 1145 deallocate(vr)
deallocate(wi) 1137 1146 deallocate(wi)
deallocate(wr) 1138 1147 deallocate(wr)
deallocate(wc_a) 1139 1148 deallocate(wc_a)
1149 deallocate(scal,rconde,rcondv)
return 1140 1150 return
end if 1141 1151 end if
1142 1152
! now fill in the results 1143 1153 ! now fill in the results
i=1 1144 1154 i=1
do while(i<=d) 1145 1155 do while(i<=d)
evala(i)=cmplx(wr(i),wi(i)) 1146 1156 evala(i)=cmplx(wr(i),wi(i))
if (wi(i) == 0.) then ! eigenvalue is real 1147 1157 if (wi(i) == 0.) then ! eigenvalue is real
eveca(:,i)=cmplx(vr(:,i),0.) 1148 1158 eveca(:,i)=cmplx(vr(:,i),0.)
else ! eigenvalue is complex 1149 1159 else ! eigenvalue is complex
evala(i+1)=cmplx(wr(i+1),wi(i+1)) 1150 1160 evala(i+1)=cmplx(wr(i+1),wi(i+1))
eveca(:,i)=cmplx(vr(:,i),vr(:,i+1)) 1151 1161 eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))
eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1)) 1152 1162 eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))
i=i+1 1153 1163 i=i+1
end if 1154 1164 end if
i=i+1 1155 1165 i=i+1
enddo 1156 1166 enddo
deallocate(work) 1157 1167 deallocate(work)
deallocate(vr) 1158 1168 deallocate(vr)
deallocate(wi) 1159 1169 deallocate(wi)
deallocate(wr) 1160 1170 deallocate(wr)
deallocate(wc_a) 1161 1171 deallocate(wc_a)
1172 deallocate(scal,rconde,rcondv)
1162 1173
! sorting 1163 1174 ! sorting
if (present(sortval) .and. sortval) then 1164 1175 if (present(sortval) .and. sortval) then
call fvn_c_sort_eigen(d,evala,eveca) 1165 1176 call fvn_c_sort_eigen(d,evala,eveca)
end if 1166 1177 end if
1167 1178
end subroutine 1168 1179 end subroutine
1169 1180
subroutine fvn_d_matev(d,a,evala,eveca,status,sortval) 1170 1181 subroutine fvn_d_matev(d,a,evala,eveca,status,sortval)
! 1171 1182 !
! integer(kind=ip_kind) d (in) : matrice rank 1172 1183 ! integer(kind=ip_kind) d (in) : matrice rank
! real(kind=dp_kind) a(d,d) (in) : The Matrix 1173 1184 ! real(kind=dp_kind) a(d,d) (in) : The Matrix
! complex(kind=dp_kind)(kind=sp_kind) evala(d) (out) : eigenvalues 1174 1185 ! complex(kind=dp_kind)(kind=sp_kind) evala(d) (out) : eigenvalues
! complex(kind=dp_kind)(kind=sp_kind) eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1175 1186 ! complex(kind=dp_kind)(kind=sp_kind) eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer(kind=ip_kind) (out) : status =0 if something went wrong 1176 1187 ! integer(kind=ip_kind) (out) : status =0 if something went wrong
! 1177 1188 !
! interfacing Lapack routine DGEEV 1178 1189 ! interfacing Lapack routine DGEEV
1190 ! 2020 11 03 : DGEEVX with BALANC='N'
implicit none 1179 1191 implicit none
integer(kind=ip_kind), intent(in) :: d 1180 1192 integer(kind=ip_kind), intent(in) :: d
real(kind=dp_kind), intent(in) :: a(d,d) 1181 1193 real(kind=dp_kind), intent(in) :: a(d,d)
complex(kind=dp_kind), intent(out) :: evala(d) 1182 1194 complex(kind=dp_kind), intent(out) :: evala(d)
complex(kind=dp_kind), intent(out) :: eveca(d,d) 1183 1195 complex(kind=dp_kind), intent(out) :: eveca(d,d)
integer(kind=ip_kind), intent(out), optional :: status 1184 1196 integer(kind=ip_kind), intent(out), optional :: status
logical(kind=ip_kind), intent(in), optional :: sortval 1185 1197 logical(kind=ip_kind), intent(in), optional :: sortval
1186 1198
real(kind=dp_kind), allocatable :: wc_a(:,:) ! a working copy of a 1187 1199 real(kind=dp_kind), allocatable :: wc_a(:,:) ! a working copy of a
integer(kind=ip_kind) :: info 1188 1200 integer(kind=ip_kind) :: info
integer(kind=ip_kind) :: lwork 1189 1201 integer(kind=ip_kind) :: lwork
real(kind=dp_kind), allocatable :: wr(:),wi(:) 1190 1202 real(kind=dp_kind), allocatable :: wr(:),wi(:)
real(kind=dp_kind) :: vl ! unused but necessary for the call 1191 1203 real(kind=dp_kind) :: vl ! unused but necessary for the call
real(kind=dp_kind), allocatable :: vr(:,:) 1192 1204 real(kind=dp_kind), allocatable :: vr(:,:)
real(kind=dp_kind), allocatable :: work(:) 1193 1205 real(kind=dp_kind), allocatable :: work(:)
real(kind=dp_kind) :: twork(1) 1194 1206 real(kind=dp_kind) :: twork(1)
integer(kind=ip_kind) i 1195 1207 integer(kind=ip_kind) i
integer(kind=ip_kind) j 1196 1208 integer(kind=ip_kind) j
1197 1209
1210 integer(kind=ip_kind) :: ilo,ihi, iwork
1211 real(kind=dp_kind), allocatable, dimension(:) :: scal, rconde, rcondv
1212 real(kind=dp_kind) :: abnrm
1213
if (present(status)) status=1 1198 1214 if (present(status)) status=1
1199 1215
! making a working copy of a 1200 1216 ! making a working copy of a
allocate(wc_a(d,d)) 1201 1217 allocate(wc_a(d,d))
!call dcopy(d*d,a,1,wc_a,1) 1202 1218 !call dcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1203 1219 wc_a(:,:)=a(:,:)
1204 1220
allocate(wr(d)) 1205 1221 allocate(wr(d))
allocate(wi(d)) 1206 1222 allocate(wi(d))
allocate(vr(d,d)) 1207 1223 allocate(vr(d,d))
1224 allocate(scal(d),rconde(d),rcondv(d))
! query optimal work size 1208 1225 ! query optimal work size
call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) 1209 1226 ! call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
1227 call dgeevx('N','N','V','N',d,wc_a,d,wr,wi,vl,1,vr,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,twork,-1,IWORK,info)
lwork=int(twork(1)) 1210 1228 lwork=int(twork(1))
allocate(work(lwork)) 1211 1229 allocate(work(lwork))
call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) 1212 1230 ! call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
1231 call dgeevx('N','N','V','N',d,wc_a,d,wr,wi,vl,1,vr,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,work,lwork,IWORK,info)
1213 1232
if (info /= 0) then 1214 1233 if (info /= 0) then
if (present(status)) status=0 1215 1234 if (present(status)) status=0
deallocate(work) 1216 1235 deallocate(work)
deallocate(vr) 1217 1236 deallocate(vr)
deallocate(wi) 1218 1237 deallocate(wi)
deallocate(wr) 1219 1238 deallocate(wr)
deallocate(wc_a) 1220 1239 deallocate(wc_a)
1240 deallocate(scal,rconde,rcondv)
return 1221 1241 return
end if 1222 1242 end if
1223 1243
! now fill in the results 1224 1244 ! now fill in the results
i=1 1225 1245 i=1
do while(i<=d) 1226 1246 do while(i<=d)
evala(i)=cmplx(wr(i),wi(i),dp_kind) 1227 1247 evala(i)=cmplx(wr(i),wi(i),dp_kind)
if (wi(i) == 0.) then ! eigenvalue is real 1228 1248 if (wi(i) == 0.) then ! eigenvalue is real
eveca(:,i)=cmplx(vr(:,i),0._dp_kind,dp_kind) 1229 1249 eveca(:,i)=cmplx(vr(:,i),0._dp_kind,dp_kind)
else ! eigenvalue is complex 1230 1250 else ! eigenvalue is complex
evala(i+1)=cmplx(wr(i+1),wi(i+1),dp_kind) 1231 1251 evala(i+1)=cmplx(wr(i+1),wi(i+1),dp_kind)
eveca(:,i)=cmplx(vr(:,i),vr(:,i+1),dp_kind) 1232 1252 eveca(:,i)=cmplx(vr(:,i),vr(:,i+1),dp_kind)
eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1),dp_kind) 1233 1253 eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1),dp_kind)
i=i+1 1234 1254 i=i+1
end if 1235 1255 end if
i=i+1 1236 1256 i=i+1
enddo 1237 1257 enddo
1238 1258
deallocate(work) 1239 1259 deallocate(work)
deallocate(vr) 1240 1260 deallocate(vr)
deallocate(wi) 1241 1261 deallocate(wi)
deallocate(wr) 1242 1262 deallocate(wr)
deallocate(wc_a) 1243 1263 deallocate(wc_a)
1264 deallocate(scal,rconde,rcondv)
1244 1265
! sorting 1245 1266 ! sorting
if (present(sortval) .and. sortval) then 1246 1267 if (present(sortval) .and. sortval) then
call fvn_z_sort_eigen(d,evala,eveca) 1247 1268 call fvn_z_sort_eigen(d,evala,eveca)
end if 1248 1269 end if
1249 1270
end subroutine 1250 1271 end subroutine
1251 1272
subroutine fvn_c_matev(d,a,evala,eveca,status,sortval) 1252 1273 subroutine fvn_c_matev(d,a,evala,eveca,status,sortval)
! 1253 1274 !
! integer(kind=ip_kind) d (in) : matrice rank 1254 1275 ! integer(kind=ip_kind) d (in) : matrice rank
! complex(kind=sp_kind) a(d,d) (in) : The Matrix 1255 1276 ! complex(kind=sp_kind) a(d,d) (in) : The Matrix
! complex(kind=sp_kind) evala(d) (out) : eigenvalues 1256 1277 ! complex(kind=sp_kind) evala(d) (out) : eigenvalues
! complex(kind=sp_kind) eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1257 1278 ! complex(kind=sp_kind) eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer(kind=ip_kind) (out) : status =0 if something went wrong 1258 1279 ! integer(kind=ip_kind) (out) : status =0 if something went wrong
! 1259 1280 !
! interfacing Lapack routine CGEEV 1260 1281 ! interfacing Lapack routine CGEEV
1282 ! 2020 11 03 : CGEEVX with BALANC='N'
implicit none 1261 1283 implicit none
integer(kind=ip_kind), intent(in) :: d 1262 1284 integer(kind=ip_kind), intent(in) :: d
complex(kind=sp_kind), intent(in) :: a(d,d) 1263 1285 complex(kind=sp_kind), intent(in) :: a(d,d)
complex(kind=sp_kind), intent(out) :: evala(d) 1264 1286 complex(kind=sp_kind), intent(out) :: evala(d)
complex(kind=sp_kind), intent(out) :: eveca(d,d) 1265 1287 complex(kind=sp_kind), intent(out) :: eveca(d,d)
integer(kind=ip_kind), intent(out), optional :: status 1266 1288 integer(kind=ip_kind), intent(out), optional :: status
logical(kind=ip_kind), intent(in), optional :: sortval 1267 1289 logical(kind=ip_kind), intent(in), optional :: sortval
1268 1290
complex(kind=sp_kind), allocatable :: wc_a(:,:) ! a working copy of a 1269 1291 complex(kind=sp_kind), allocatable :: wc_a(:,:) ! a working copy of a
integer(kind=ip_kind) :: info 1270 1292 integer(kind=ip_kind) :: info
integer(kind=ip_kind) :: lwork 1271 1293 integer(kind=ip_kind) :: lwork
complex(kind=sp_kind), allocatable :: work(:) 1272 1294 complex(kind=sp_kind), allocatable :: work(:)
complex(kind=sp_kind) :: twork(1) 1273 1295 complex(kind=sp_kind) :: twork(1)
real(kind=sp_kind), allocatable :: rwork(:) 1274 1296 real(kind=sp_kind), allocatable :: rwork(:)
complex(kind=sp_kind) :: vl ! unused but necessary for the call 1275 1297 complex(kind=sp_kind) :: vl ! unused but necessary for the call
1276 1298
1299 integer(kind=ip_kind) :: ilo,ihi
1300 real(kind=sp_kind), allocatable, dimension(:) :: scal, rconde, rcondv
1301 real(kind=sp_kind) :: abnrm
1302
if (present(status)) status=1 1277 1303 if (present(status)) status=1
1278 1304
! making a working copy of a 1279 1305 ! making a working copy of a
allocate(wc_a(d,d)) 1280 1306 allocate(wc_a(d,d))
!call ccopy(d*d,a,1,wc_a,1) 1281 1307 !call ccopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1282 1308 wc_a(:,:)=a(:,:)
1283 1309
! rwork must be allocated before query 1284 1310 ! rwork must be allocated before query
allocate(rwork(2*d)) 1285 1311 allocate(rwork(2*d))
1312 allocate(scal(d),rconde(d),rcondv(d))
! query optimal work size 1286 1313 ! query optimal work size
call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) 1287 1314 ! call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
1315 call cgeevx('N','N','V','N',d,wc_a,d,evala,vl,1,eveca,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,twork,-1,rwork,info)
lwork=int(twork(1)) 1288 1316 lwork=int(twork(1))
allocate(work(lwork)) 1289 1317 allocate(work(lwork))
call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) 1290 1318 ! call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
1319 call cgeevx('N','N','V','N',d,wc_a,d,evala,vl,1,eveca,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,work,lwork,rwork,info)
1291 1320
if (info /= 0) then 1292 1321 if (info /= 0) then
if (present(status)) status=0 1293 1322 if (present(status)) status=0
end if 1294 1323 end if
deallocate(rwork) 1295 1324 deallocate(rwork)
deallocate(work) 1296 1325 deallocate(work)
deallocate(wc_a) 1297 1326 deallocate(wc_a)
1327 deallocate(scal,rconde,rcondv)
1298 1328
! sorting 1299 1329 ! sorting
if (present(sortval) .and. sortval) then 1300 1330 if (present(sortval) .and. sortval) then
call fvn_c_sort_eigen(d,evala,eveca) 1301 1331 call fvn_c_sort_eigen(d,evala,eveca)
end if 1302 1332 end if
1303 1333
end subroutine 1304 1334 end subroutine
1305 1335
subroutine fvn_z_matev(d,a,evala,eveca,status,sortval) 1306 1336 subroutine fvn_z_matev(d,a,evala,eveca,status,sortval)
! 1307 1337 !
! integer(kind=ip_kind) d (in) : matrice rank 1308 1338 ! integer(kind=ip_kind) d (in) : matrice rank
! complex(kind=dp_kind)(kind=sp_kind) a(d,d) (in) : The Matrix 1309 1339 ! complex(kind=dp_kind)(kind=sp_kind) a(d,d) (in) : The Matrix
! complex(kind=dp_kind)(kind=sp_kind) evala(d) (out) : eigenvalues 1310 1340 ! complex(kind=dp_kind)(kind=sp_kind) evala(d) (out) : eigenvalues
! complex(kind=dp_kind)(kind=sp_kind) eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1311 1341 ! complex(kind=dp_kind)(kind=sp_kind) eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer(kind=ip_kind) (out) : status =0 if something went wrong 1312 1342 ! integer(kind=ip_kind) (out) : status =0 if something went wrong
! 1313 1343 !
! interfacing Lapack routine ZGEEV 1314 1344 ! interfacing Lapack routine ZGEEV
1345 ! 2020 11 03 : ZGEEVX with BALANC='N'
implicit none 1315 1346 implicit none
integer(kind=ip_kind), intent(in) :: d 1316 1347 integer(kind=ip_kind), intent(in) :: d
complex(kind=dp_kind), intent(in) :: a(d,d) 1317 1348 complex(kind=dp_kind), intent(in) :: a(d,d)
complex(kind=dp_kind), intent(out) :: evala(d) 1318 1349 complex(kind=dp_kind), intent(out) :: evala(d)
complex(kind=dp_kind), intent(out) :: eveca(d,d) 1319 1350 complex(kind=dp_kind), intent(out) :: eveca(d,d)
integer(kind=ip_kind), intent(out), optional :: status 1320 1351 integer(kind=ip_kind), intent(out), optional :: status
logical(kind=ip_kind), intent(in), optional :: sortval 1321 1352 logical(kind=ip_kind), intent(in), optional :: sortval
1322 1353
complex(kind=dp_kind), allocatable :: wc_a(:,:) ! a working copy of a 1323 1354 complex(kind=dp_kind), allocatable :: wc_a(:,:) ! a working copy of a
integer(kind=ip_kind) :: info 1324 1355 integer(kind=ip_kind) :: info
integer(kind=ip_kind) :: lwork 1325 1356 integer(kind=ip_kind) :: lwork
complex(kind=dp_kind), allocatable :: work(:) 1326 1357 complex(kind=dp_kind), allocatable :: work(:)
complex(kind=dp_kind) :: twork(1) 1327 1358 complex(kind=dp_kind) :: twork(1)
real(kind=dp_kind), allocatable :: rwork(:) 1328 1359 real(kind=dp_kind), allocatable :: rwork(:)
complex(kind=dp_kind) :: vl ! unused but necessary for the call 1329 1360 complex(kind=dp_kind) :: vl ! unused but necessary for the call
1330 1361
1362 integer(kind=ip_kind) :: ilo,ihi
1363 real(kind=dp_kind), allocatable, dimension(:) :: scal, rconde, rcondv
1364 real(kind=dp_kind) :: abnrm
1365
1366
if (present(status)) status=1 1331 1367 if (present(status)) status=1
1332 1368
! making a working copy of a 1333 1369 ! making a working copy of a
allocate(wc_a(d,d)) 1334 1370 allocate(wc_a(d,d))
!call zcopy(d*d,a,1,wc_a,1) 1335 1371 !call zcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1336 1372 wc_a(:,:)=a(:,:)
1337 1373
! rwork must be allocated before query 1338 1374 ! rwork must be allocated before query
allocate(rwork(2*d)) 1339 1375 allocate(rwork(2*d))
1376 allocate(scal(d),rconde(d),rcondv(d))
! query optimal work size 1340 1377 ! query optimal work size
call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) 1341 1378 ! call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
1379 call zgeevx('N','N','V','N',d,wc_a,d,evala,vl,1,eveca,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,twork,-1,rwork,info)
lwork=int(twork(1)) 1342 1380 lwork=int(twork(1))
allocate(work(lwork)) 1343 1381 allocate(work(lwork))
call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) 1344 1382 ! call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
1383 call zgeevx('N','N','V','N',d,wc_a,d,evala,vl,1,eveca,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,work,lwork,rwork,info)
1345 1384
if (info /= 0) then 1346 1385 if (info /= 0) then
if (present(status)) status=0 1347 1386 if (present(status)) status=0
end if 1348 1387 end if
deallocate(rwork) 1349 1388 deallocate(rwork)
deallocate(work) 1350 1389 deallocate(work)
deallocate(wc_a) 1351 1390 deallocate(wc_a)
1391 deallocate(scal,rconde,rcondv)
1352 1392
! sorting 1353 1393 ! sorting
if (present(sortval) .and. sortval) then 1354 1394 if (present(sortval) .and. sortval) then
call fvn_z_sort_eigen(d,evala,eveca) 1355 1395 call fvn_z_sort_eigen(d,evala,eveca)
end if 1356 1396 end if
end subroutine 1357 1397 end subroutine
1358 1398
1359 1399
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1360 1400 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1361 1401 !
! Least square problem 1362 1402 ! Least square problem
! 1363 1403 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1364 1404 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1365 1405 !
! 1366 1406 !
1367 1407
1368 1408
1369 1409
1370 1410
subroutine fvn_d_lspoly(np,x,y,deg,coeff,status) 1371 1411 subroutine fvn_d_lspoly(np,x,y,deg,coeff,status)
! 1372 1412 !
! Least square polynomial fitting 1373 1413 ! Least square polynomial fitting
! 1374 1414 !
! Find the coefficients of the least square polynomial of a given degree 1375 1415 ! Find the coefficients of the least square polynomial of a given degree
! for a set of coordinates. 1376 1416 ! for a set of coordinates.
! 1377 1417 !
! The degree must be lower than the number of points 1378 1418 ! The degree must be lower than the number of points
! 1379 1419 !
! np (in) : number of points 1380 1420 ! np (in) : number of points
! x(np) (in) : x data 1381 1421 ! x(np) (in) : x data
! y(np) (in) : y data 1382 1422 ! y(np) (in) : y data
! deg (in) : polynomial's degree 1383 1423 ! deg (in) : polynomial's degree
! coeff(deg+1) (out) : polynomial coefficients 1384 1424 ! coeff(deg+1) (out) : polynomial coefficients
! status (out) : =0 if a problem occurs 1385 1425 ! status (out) : =0 if a problem occurs
! 1386 1426 !
implicit none 1387 1427 implicit none