lxsd / fvn

Download as
Browse Code ยป

Commit b85eed0171def3a9c57818323d3494afee5fc0ed

Authored by wdaniau
1 parent 19e954a45e
Exists in master and in 3 other branches geevx, geevx_2020, multi

Added sorting of eigenvalues and vectors using an insertion sort algo to better 

match imsl behavior

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

Showing 1 changed file with 73 additions and 0 deletions Inline Diff

fvn_linear/fvn_linear.f90
Diff comments   View file @ b85eed0
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 !
! 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(4),dimension(:,:),intent(in) :: a 222 222 real(4),dimension(:,:),intent(in) :: a
real(4),dimension(size(a,1),size(a,1)) :: fvn_op_s_matinv,tmp_array 223 223 real(4),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(8),dimension(:,:),intent(in) :: a 229 229 real(8),dimension(:,:),intent(in) :: a
real(8),dimension(size(a,1),size(a,1)) :: fvn_op_d_matinv,tmp_array 230 230 real(8),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(4),dimension(:,:),intent(in) :: a 236 236 complex(4),dimension(:,:),intent(in) :: a
complex(4),dimension(size(a,1),size(a,1)) :: fvn_op_c_matinv,tmp_array 237 237 complex(4),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(8),dimension(:,:),intent(in) :: a 243 243 complex(8),dimension(:,:),intent(in) :: a
complex(8),dimension(size(a,1),size(a,1)) :: fvn_op_z_matinv,tmp_array 244 244 complex(8),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(4), dimension(:,:), intent(in) :: a,b 254 254 real(4), dimension(:,:), intent(in) :: a,b
real(4), dimension(size(a,1),size(b,2)) :: fvn_op_s_ix 255 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) 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(8), dimension(:,:), intent(in) :: a,b 260 260 real(8), dimension(:,:), intent(in) :: a,b
real(8), dimension(size(a,1),size(b,2)) :: fvn_op_d_ix 261 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) 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(4), dimension(:,:), intent(in) :: a,b 266 266 complex(4), dimension(:,:), intent(in) :: a,b
complex(4), dimension(size(a,1),size(b,2)) :: fvn_op_c_ix 267 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) 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(8), dimension(:,:), intent(in) :: a,b 272 272 complex(8), dimension(:,:), intent(in) :: a,b
complex(8), dimension(size(a,1),size(b,2)) :: fvn_op_z_ix 273 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) 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(4), dimension(:,:), intent(in) :: a,b 282 282 real(4), dimension(:,:), intent(in) :: a,b
real(4), dimension(size(a,1),size(b,2)) :: fvn_op_s_xi 283 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)) 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(8), dimension(:,:), intent(in) :: a,b 288 288 real(8), dimension(:,:), intent(in) :: a,b
real(8), dimension(size(a,1),size(b,2)) :: fvn_op_d_xi 289 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)) 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(4), dimension(:,:), intent(in) :: a,b 294 294 complex(4), dimension(:,:), intent(in) :: a,b
complex(4), dimension(size(a,1),size(b,2)) :: fvn_op_c_xi 295 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)) 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(8), dimension(:,:), intent(in) :: a,b 300 300 complex(8), dimension(:,:), intent(in) :: a,b
complex(8), dimension(size(a,1),size(b,2)) :: fvn_op_z_xi 301 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)) 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(4),dimension(:,:),intent(in) :: a 310 310 complex(4),dimension(:,:),intent(in) :: a
complex(4),dimension(size(a,2),size(a,1)) :: fvn_op_c_conj_transpose 311 311 complex(4),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(8),dimension(:,:),intent(in) :: a 316 316 complex(8),dimension(:,:),intent(in) :: a
complex(8),dimension(size(a,2),size(a,1)) :: fvn_op_z_conj_transpose 317 317 complex(8),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(4), dimension(:,:), intent(in) :: a,b 327 327 complex(4), dimension(:,:), intent(in) :: a,b
complex(4), dimension(size(a,2),size(b,2)) :: fvn_op_c_hx 328 328 complex(4), 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(8), dimension(:,:), intent(in) :: a,b 333 333 complex(8), dimension(:,:), intent(in) :: a,b
complex(8), dimension(size(a,2),size(b,2)) :: fvn_op_z_hx 334 334 complex(8), 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(4), dimension(:,:), intent(in) :: a,b 344 344 complex(4), dimension(:,:), intent(in) :: a,b
complex(4), dimension(size(a,1),size(b,1)) :: fvn_op_c_xh 345 345 complex(4), 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(8), dimension(:,:), intent(in) :: a,b 350 350 complex(8), dimension(:,:), intent(in) :: a,b
complex(8), dimension(size(a,1),size(b,1)) :: fvn_op_z_xh 351 351 complex(8), 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=4) :: n 366 366 integer(kind=4) :: n
real(kind=8), dimension(n,n) :: fvn_d_ident 367 367 real(kind=8), dimension(n,n) :: fvn_d_ident
368 368
real(kind=8),dimension(n*n) :: vect 369 369 real(kind=8),dimension(n*n) :: vect
integer(kind=4) :: i 370 370 integer(kind=4) :: i
371 371
vect=0._8 372 372 vect=0._8
vect(1:n*n:n+1) = 1._8 373 373 vect(1:n*n:n+1) = 1._8
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=4) :: n 379 379 integer(kind=4) :: n
real(kind=4), dimension(n,n) :: fvn_s_ident 380 380 real(kind=4), dimension(n,n) :: fvn_s_ident
381 381
real(kind=4),dimension(n*n) :: vect 382 382 real(kind=4),dimension(n*n) :: vect
integer(kind=4) :: i 383 383 integer(kind=4) :: i
384 384
vect=0._4 385 385 vect=0._4
vect(1:n*n:n+1) = 1._4 386 386 vect(1:n*n:n+1) = 1._4
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=4) :: n 392 392 integer(kind=4) :: n
complex(kind=4), dimension(n,n) :: fvn_c_ident 393 393 complex(kind=4), dimension(n,n) :: fvn_c_ident
394 394
complex(kind=4),dimension(n*n) :: vect 395 395 complex(kind=4),dimension(n*n) :: vect
integer(kind=4) :: i 396 396 integer(kind=4) :: i
397 397
vect=(0._4,0._4) 398 398 vect=(0._4,0._4)
vect(1:n*n:n+1) = (1._4,0._4) 399 399 vect(1:n*n:n+1) = (1._4,0._4)
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=4) :: n 405 405 integer(kind=4) :: n
complex(kind=8), dimension(n,n) :: fvn_z_ident 406 406 complex(kind=8), dimension(n,n) :: fvn_z_ident
407 407
complex(kind=8),dimension(n*n) :: vect 408 408 complex(kind=8),dimension(n*n) :: vect
integer(kind=4) :: i 409 409 integer(kind=4) :: i
410 410
vect=(0._8,0._8) 411 411 vect=(0._8,0._8)
vect(1:n*n:n+1) = (1._8,0._8) 412 412 vect(1:n*n:n+1) = (1._8,0._8)
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 matrix using BLAS and LAPACK 426 426 ! Matrix inversion of a real 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, intent(in) :: d 434 434 integer, intent(in) :: d
real, intent(in) :: a(d,d) 435 435 real, intent(in) :: a(d,d)
real, intent(out) :: inva(d,d) 436 436 real, intent(out) :: inva(d,d)
integer, intent(out),optional :: status 437 437 integer, intent(out),optional :: status
438 438
integer, allocatable :: ipiv(:) 439 439 integer, allocatable :: ipiv(:)
real, allocatable :: work(:) 440 440 real, allocatable :: work(:)
real twork(1) 441 441 real twork(1)
integer :: info 442 442 integer :: info
integer :: lwork 443 443 integer :: 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 double precision matrix using BLAS and LAPACK 475 475 ! Matrix inversion of a double precision 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, intent(in) :: d 483 483 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 484 484 double precision, intent(in) :: a(d,d)
double precision, intent(out) :: inva(d,d) 485 485 double precision, intent(out) :: inva(d,d)
integer, intent(out),optional :: status 486 486 integer, intent(out),optional :: status
487 487
integer, allocatable :: ipiv(:) 488 488 integer, allocatable :: ipiv(:)
double precision, allocatable :: work(:) 489 489 double precision, allocatable :: work(:)
double precision :: twork(1) 490 490 double precision :: twork(1)
integer :: info 491 491 integer :: info
integer :: lwork 492 492 integer :: 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 matrix using BLAS and LAPACK 524 524 ! Matrix inversion of a complex 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, intent(in) :: d 532 532 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 533 533 complex, intent(in) :: a(d,d)
complex, intent(out) :: inva(d,d) 534 534 complex, intent(out) :: inva(d,d)
integer, intent(out),optional :: status 535 535 integer, intent(out),optional :: status
536 536
integer, allocatable :: ipiv(:) 537 537 integer, allocatable :: ipiv(:)
complex, allocatable :: work(:) 538 538 complex, allocatable :: work(:)
complex :: twork(1) 539 539 complex :: twork(1)
integer :: info 540 540 integer :: info
integer :: lwork 541 541 integer :: 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 double complex matrix using BLAS and LAPACK 574 574 ! Matrix inversion of a double complex 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, intent(in) :: d 582 582 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 583 583 double complex, intent(in) :: a(d,d)
double complex, intent(out) :: inva(d,d) 584 584 double complex, intent(out) :: inva(d,d)
integer, intent(out),optional :: status 585 585 integer, intent(out),optional :: status
586 586
integer, allocatable :: ipiv(:) 587 587 integer, allocatable :: ipiv(:)
double complex, allocatable :: work(:) 588 588 double complex, allocatable :: work(:)
double complex :: twork(1) 589 589 double complex :: twork(1)
integer :: info 590 590 integer :: info
integer :: lwork 591 591 integer :: 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, intent(in) :: d 637 637 integer, intent(in) :: d
real, intent(in) :: a(d,d) 638 638 real, intent(in) :: a(d,d)
integer, intent(out), optional :: status 639 639 integer, intent(out), optional :: status
real :: fvn_s_det 640 640 real :: fvn_s_det
641 641
real, allocatable :: wc_a(:,:) 642 642 real, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 643 643 integer, allocatable :: ipiv(:)
integer :: info,i 644 644 integer :: 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, intent(in) :: d 680 680 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 681 681 double precision, intent(in) :: a(d,d)
integer, intent(out), optional :: status 682 682 integer, intent(out), optional :: status
double precision :: fvn_d_det 683 683 double precision :: fvn_d_det
684 684
double precision, allocatable :: wc_a(:,:) 685 685 double precision, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 686 686 integer, allocatable :: ipiv(:)
integer :: info,i 687 687 integer :: 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, intent(in) :: d 722 722 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 723 723 complex, intent(in) :: a(d,d)
integer, intent(out), optional :: status 724 724 integer, intent(out), optional :: status
complex :: fvn_c_det 725 725 complex :: fvn_c_det
726 726
complex, allocatable :: wc_a(:,:) 727 727 complex, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 728 728 integer, allocatable :: ipiv(:)
integer :: info,i 729 729 integer :: 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, intent(in) :: d 766 766 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 767 767 double complex, intent(in) :: a(d,d)
integer, intent(out), optional :: status 768 768 integer, intent(out), optional :: status
double complex :: fvn_z_det 769 769 double complex :: fvn_z_det
770 770
double complex, allocatable :: wc_a(:,:) 771 771 double complex, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 772 772 integer, allocatable :: ipiv(:)
integer :: info,i 773 773 integer :: 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, intent(in) :: d 818 818 integer, intent(in) :: d
real, intent(in) :: a(d,d) 819 819 real, intent(in) :: a(d,d)
real, intent(out) :: rcond 820 820 real, intent(out) :: rcond
integer, intent(out), optional :: status 821 821 integer, intent(out), optional :: status
822 822
real, allocatable :: work(:) 823 823 real, allocatable :: work(:)
integer, allocatable :: iwork(:) 824 824 integer, allocatable :: iwork(:)
real :: anorm 825 825 real :: anorm
real, allocatable :: wc_a(:,:) ! working copy of a 826 826 real, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 827 827 integer :: info
integer, allocatable :: ipiv(:) 828 828 integer, allocatable :: ipiv(:)
829 829
real, external :: slange 830 830 real, 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, intent(in) :: d 871 871 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 872 872 double precision, intent(in) :: a(d,d)
double precision, intent(out) :: rcond 873 873 double precision, intent(out) :: rcond
integer, intent(out), optional :: status 874 874 integer, intent(out), optional :: status
875 875
double precision, allocatable :: work(:) 876 876 double precision, allocatable :: work(:)
integer, allocatable :: iwork(:) 877 877 integer, allocatable :: iwork(:)
double precision :: anorm 878 878 double precision :: anorm
double precision, allocatable :: wc_a(:,:) ! working copy of a 879 879 double precision, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 880 880 integer :: info
integer, allocatable :: ipiv(:) 881 881 integer, allocatable :: ipiv(:)
882 882
double precision, external :: dlange 883 883 double precision, 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, intent(in) :: d 925 925 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 926 926 complex, intent(in) :: a(d,d)
real, intent(out) :: rcond 927 927 real, intent(out) :: rcond
integer, intent(out), optional :: status 928 928 integer, intent(out), optional :: status
929 929
real, allocatable :: rwork(:) 930 930 real, allocatable :: rwork(:)
complex, allocatable :: work(:) 931 931 complex, allocatable :: work(:)
integer, allocatable :: iwork(:) 932 932 integer, allocatable :: iwork(:)
real :: anorm 933 933 real :: anorm
complex, allocatable :: wc_a(:,:) ! working copy of a 934 934 complex, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 935 935 integer :: info
integer, allocatable :: ipiv(:) 936 936 integer, allocatable :: ipiv(:)
937 937
real, external :: clange 938 938 real, 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, intent(in) :: d 978 978 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 979 979 double complex, intent(in) :: a(d,d)
double precision, intent(out) :: rcond 980 980 double precision, intent(out) :: rcond
integer, intent(out), optional :: status 981 981 integer, intent(out), optional :: status
982 982
double complex, allocatable :: work(:) 983 983 double complex, allocatable :: work(:)
double precision, allocatable :: rwork(:) 984 984 double precision, allocatable :: rwork(:)
double precision :: anorm 985 985 double precision :: anorm
double complex, allocatable :: wc_a(:,:) ! working copy of a 986 986 double complex, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 987 987 integer :: info
integer, allocatable :: ipiv(:) 988 988 integer, allocatable :: ipiv(:)
989 989
double precision, external :: zlange 990 990 double precision, 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
1028 ! August 2009
1029 ! William Daniau
1030 ! Adding sorting of eigenvalues and vectors
1031 subroutine fvn_z_sort_eigen(d,v,vec)
1032 ! this routine takes in input :
1033 ! v : a vector containing unsorted eigenvalues
1034 ! vec : a matrix where vec(:,j) is the eigenvector corresponding to v(j)
1035 !
1036 ! At the end of subroutine the eigenvalues are sorted by decreasing order of
1037 ! modulus so as the vectors.
1038 implicit none
1039 integer :: d,i,j
1040 complex(8),dimension(d) :: v
1041 complex(8),dimension(d,d) :: vec
1042 complex(8) :: cur_v
1043 complex(8), dimension(d) :: cur_vec
1044
1045 do i=2,d
1046 cur_v=v(i)
1047 cur_vec=vec(:,i)
1048 j=i-1
1049 do while( (j>=1) .and. (abs(v(j)) < abs(cur_v)) )
1050 v(j+1)=v(j)
1051 vec(:,j+1)=vec(:,j)
1052 j=j-1
1053 end do
1054 v(j+1)=cur_v
1055 vec(:,j+1)=cur_vec
1056 end do
1057
1058 end subroutine
1059
1060 subroutine fvn_c_sort_eigen(d,v,vec)
1061 ! this routine takes in input :
1062 ! v : a vector containing unsorted eigenvalues
1063 ! vec : a matrix where vec(:,j) is the eigenvector corresponding to v(j)
1064 !
1065 ! At the end of subroutine the eigenvalues are sorted by decreasing order of
1066 ! modulus so as the vectors.
1067 implicit none
1068 integer :: d,i,j
1069 complex(4),dimension(d) :: v
1070 complex(4),dimension(d,d) :: vec
1071 complex(4) :: cur_v
1072 complex(4), dimension(d) :: cur_vec
1073
1074 do i=2,d
1075 cur_v=v(i)
1076 cur_vec=vec(:,i)
1077 j=i-1
1078 do while( (j>=1) .and. (abs(v(j)) < abs(cur_v)) )
1079 v(j+1)=v(j)
1080 vec(:,j+1)=vec(:,j)
1081 j=j-1
1082 end do
1083 v(j+1)=cur_v
1084 vec(:,j+1)=cur_vec
1085 end do
1086
1087 end subroutine
1088
subroutine fvn_s_matev(d,a,evala,eveca,status) 1028 1089 subroutine fvn_s_matev(d,a,evala,eveca,status)
! 1029 1090 !
! integer d (in) : matrice rank 1030 1091 ! integer d (in) : matrice rank
! real a(d,d) (in) : The Matrix 1031 1092 ! real a(d,d) (in) : The Matrix
! complex evala(d) (out) : eigenvalues 1032 1093 ! complex evala(d) (out) : eigenvalues
! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1033 1094 ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 1034 1095 ! integer (out) : status =0 if something went wrong
! 1035 1096 !
! interfacing Lapack routine SGEEV 1036 1097 ! interfacing Lapack routine SGEEV
implicit none 1037 1098 implicit none
integer, intent(in) :: d 1038 1099 integer, intent(in) :: d
real, intent(in) :: a(d,d) 1039 1100 real, intent(in) :: a(d,d)
complex, intent(out) :: evala(d) 1040 1101 complex, intent(out) :: evala(d)
complex, intent(out) :: eveca(d,d) 1041 1102 complex, intent(out) :: eveca(d,d)
integer, intent(out), optional :: status 1042 1103 integer, intent(out), optional :: status
1043 1104
real, allocatable :: wc_a(:,:) ! a working copy of a 1044 1105 real, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 1045 1106 integer :: info
integer :: lwork 1046 1107 integer :: lwork
real, allocatable :: wr(:),wi(:) 1047 1108 real, allocatable :: wr(:),wi(:)
real :: vl ! unused but necessary for the call 1048 1109 real :: vl ! unused but necessary for the call
real, allocatable :: vr(:,:) 1049 1110 real, allocatable :: vr(:,:)
real, allocatable :: work(:) 1050 1111 real, allocatable :: work(:)
real :: twork(1) 1051 1112 real :: twork(1)
integer i 1052 1113 integer i
integer j 1053 1114 integer j
1054 1115
if (present(status)) status=1 1055 1116 if (present(status)) status=1
1056 1117
! making a working copy of a 1057 1118 ! making a working copy of a
allocate(wc_a(d,d)) 1058 1119 allocate(wc_a(d,d))
!call scopy(d*d,a,1,wc_a,1) 1059 1120 !call scopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1060 1121 wc_a(:,:)=a(:,:)
1061 1122
allocate(wr(d)) 1062 1123 allocate(wr(d))
allocate(wi(d)) 1063 1124 allocate(wi(d))
allocate(vr(d,d)) 1064 1125 allocate(vr(d,d))
! query optimal work size 1065 1126 ! query optimal work size
call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) 1066 1127 call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
lwork=int(twork(1)) 1067 1128 lwork=int(twork(1))
allocate(work(lwork)) 1068 1129 allocate(work(lwork))
call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) 1069 1130 call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
1070 1131
if (info /= 0) then 1071 1132 if (info /= 0) then
if (present(status)) status=0 1072 1133 if (present(status)) status=0
deallocate(work) 1073 1134 deallocate(work)
deallocate(vr) 1074 1135 deallocate(vr)
deallocate(wi) 1075 1136 deallocate(wi)
deallocate(wr) 1076 1137 deallocate(wr)
deallocate(wc_a) 1077 1138 deallocate(wc_a)
return 1078 1139 return
end if 1079 1140 end if
1080 1141
! now fill in the results 1081 1142 ! now fill in the results
i=1 1082 1143 i=1
do while(i<=d) 1083 1144 do while(i<=d)
evala(i)=cmplx(wr(i),wi(i)) 1084 1145 evala(i)=cmplx(wr(i),wi(i))
if (wi(i) == 0.) then ! eigenvalue is real 1085 1146 if (wi(i) == 0.) then ! eigenvalue is real
eveca(:,i)=cmplx(vr(:,i),0.) 1086 1147 eveca(:,i)=cmplx(vr(:,i),0.)
else ! eigenvalue is complex 1087 1148 else ! eigenvalue is complex
evala(i+1)=cmplx(wr(i+1),wi(i+1)) 1088 1149 evala(i+1)=cmplx(wr(i+1),wi(i+1))
eveca(:,i)=cmplx(vr(:,i),vr(:,i+1)) 1089 1150 eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))
eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1)) 1090 1151 eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))
i=i+1 1091 1152 i=i+1
end if 1092 1153 end if
i=i+1 1093 1154 i=i+1
enddo 1094 1155 enddo
deallocate(work) 1095 1156 deallocate(work)
deallocate(vr) 1096 1157 deallocate(vr)
deallocate(wi) 1097 1158 deallocate(wi)
deallocate(wr) 1098 1159 deallocate(wr)
deallocate(wc_a) 1099 1160 deallocate(wc_a)
1100 1161
1162 ! sorting
1163 call fvn_c_sort_eigen(d,evala,eveca)
1164
end subroutine 1101 1165 end subroutine
1102 1166
subroutine fvn_d_matev(d,a,evala,eveca,status) 1103 1167 subroutine fvn_d_matev(d,a,evala,eveca,status)
! 1104 1168 !
! integer d (in) : matrice rank 1105 1169 ! integer d (in) : matrice rank
! double precision a(d,d) (in) : The Matrix 1106 1170 ! double precision a(d,d) (in) : The Matrix
! double complex evala(d) (out) : eigenvalues 1107 1171 ! double complex evala(d) (out) : eigenvalues
! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1108 1172 ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 1109 1173 ! integer (out) : status =0 if something went wrong
! 1110 1174 !
! interfacing Lapack routine DGEEV 1111 1175 ! interfacing Lapack routine DGEEV
implicit none 1112 1176 implicit none
integer, intent(in) :: d 1113 1177 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 1114 1178 double precision, intent(in) :: a(d,d)
double complex, intent(out) :: evala(d) 1115 1179 double complex, intent(out) :: evala(d)
double complex, intent(out) :: eveca(d,d) 1116 1180 double complex, intent(out) :: eveca(d,d)
integer, intent(out), optional :: status 1117 1181 integer, intent(out), optional :: status
1118 1182
double precision, allocatable :: wc_a(:,:) ! a working copy of a 1119 1183 double precision, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 1120 1184 integer :: info
integer :: lwork 1121 1185 integer :: lwork
double precision, allocatable :: wr(:),wi(:) 1122 1186 double precision, allocatable :: wr(:),wi(:)
double precision :: vl ! unused but necessary for the call 1123 1187 double precision :: vl ! unused but necessary for the call
double precision, allocatable :: vr(:,:) 1124 1188 double precision, allocatable :: vr(:,:)
double precision, allocatable :: work(:) 1125 1189 double precision, allocatable :: work(:)
double precision :: twork(1) 1126 1190 double precision :: twork(1)
integer i 1127 1191 integer i
integer j 1128 1192 integer j
1129 1193
if (present(status)) status=1 1130 1194 if (present(status)) status=1
1131 1195
! making a working copy of a 1132 1196 ! making a working copy of a
allocate(wc_a(d,d)) 1133 1197 allocate(wc_a(d,d))
!call dcopy(d*d,a,1,wc_a,1) 1134 1198 !call dcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1135 1199 wc_a(:,:)=a(:,:)
1136 1200
allocate(wr(d)) 1137 1201 allocate(wr(d))
allocate(wi(d)) 1138 1202 allocate(wi(d))
allocate(vr(d,d)) 1139 1203 allocate(vr(d,d))
! query optimal work size 1140 1204 ! query optimal work size
call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) 1141 1205 call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
lwork=int(twork(1)) 1142 1206 lwork=int(twork(1))
allocate(work(lwork)) 1143 1207 allocate(work(lwork))
call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) 1144 1208 call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
1145 1209
if (info /= 0) then 1146 1210 if (info /= 0) then
if (present(status)) status=0 1147 1211 if (present(status)) status=0
deallocate(work) 1148 1212 deallocate(work)
deallocate(vr) 1149 1213 deallocate(vr)
deallocate(wi) 1150 1214 deallocate(wi)
deallocate(wr) 1151 1215 deallocate(wr)
deallocate(wc_a) 1152 1216 deallocate(wc_a)
return 1153 1217 return
end if 1154 1218 end if
1155 1219
! now fill in the results 1156 1220 ! now fill in the results
i=1 1157 1221 i=1
do while(i<=d) 1158 1222 do while(i<=d)
evala(i)=dcmplx(wr(i),wi(i)) 1159 1223 evala(i)=dcmplx(wr(i),wi(i))
if (wi(i) == 0.) then ! eigenvalue is real 1160 1224 if (wi(i) == 0.) then ! eigenvalue is real
eveca(:,i)=dcmplx(vr(:,i),0.) 1161 1225 eveca(:,i)=dcmplx(vr(:,i),0.)
else ! eigenvalue is complex 1162 1226 else ! eigenvalue is complex
evala(i+1)=dcmplx(wr(i+1),wi(i+1)) 1163 1227 evala(i+1)=dcmplx(wr(i+1),wi(i+1))
eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1)) 1164 1228 eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))
eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1)) 1165 1229 eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))
i=i+1 1166 1230 i=i+1
end if 1167 1231 end if
i=i+1 1168 1232 i=i+1
enddo 1169 1233 enddo
1170 1234
deallocate(work) 1171 1235 deallocate(work)
deallocate(vr) 1172 1236 deallocate(vr)
deallocate(wi) 1173 1237 deallocate(wi)
deallocate(wr) 1174 1238 deallocate(wr)
deallocate(wc_a) 1175 1239 deallocate(wc_a)
1176 1240
1241 ! sorting
1242 call fvn_z_sort_eigen(d,evala,eveca)
1243
end subroutine 1177 1244 end subroutine
1178 1245
subroutine fvn_c_matev(d,a,evala,eveca,status) 1179 1246 subroutine fvn_c_matev(d,a,evala,eveca,status)
! 1180 1247 !
! integer d (in) : matrice rank 1181 1248 ! integer d (in) : matrice rank
! complex a(d,d) (in) : The Matrix 1182 1249 ! complex a(d,d) (in) : The Matrix
! complex evala(d) (out) : eigenvalues 1183 1250 ! complex evala(d) (out) : eigenvalues
! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1184 1251 ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 1185 1252 ! integer (out) : status =0 if something went wrong
! 1186 1253 !
! interfacing Lapack routine CGEEV 1187 1254 ! interfacing Lapack routine CGEEV
implicit none 1188 1255 implicit none
integer, intent(in) :: d 1189 1256 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 1190 1257 complex, intent(in) :: a(d,d)
complex, intent(out) :: evala(d) 1191 1258 complex, intent(out) :: evala(d)
complex, intent(out) :: eveca(d,d) 1192 1259 complex, intent(out) :: eveca(d,d)
integer, intent(out), optional :: status 1193 1260 integer, intent(out), optional :: status
1194 1261
complex, allocatable :: wc_a(:,:) ! a working copy of a 1195 1262 complex, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 1196 1263 integer :: info
integer :: lwork 1197 1264 integer :: lwork
complex, allocatable :: work(:) 1198 1265 complex, allocatable :: work(:)
complex :: twork(1) 1199 1266 complex :: twork(1)
real, allocatable :: rwork(:) 1200 1267 real, allocatable :: rwork(:)
complex :: vl ! unused but necessary for the call 1201 1268 complex :: vl ! unused but necessary for the call
1202 1269
if (present(status)) status=1 1203 1270 if (present(status)) status=1
1204 1271
! making a working copy of a 1205 1272 ! making a working copy of a
allocate(wc_a(d,d)) 1206 1273 allocate(wc_a(d,d))
!call ccopy(d*d,a,1,wc_a,1) 1207 1274 !call ccopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1208 1275 wc_a(:,:)=a(:,:)
1209 1276
1210 1277
! query optimal work size 1211 1278 ! query optimal work size
call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) 1212 1279 call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
lwork=int(twork(1)) 1213 1280 lwork=int(twork(1))
allocate(work(lwork)) 1214 1281 allocate(work(lwork))
allocate(rwork(2*d)) 1215 1282 allocate(rwork(2*d))
call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) 1216 1283 call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
1217 1284
if (info /= 0) then 1218 1285 if (info /= 0) then
if (present(status)) status=0 1219 1286 if (present(status)) status=0
end if 1220 1287 end if
deallocate(rwork) 1221 1288 deallocate(rwork)
deallocate(work) 1222 1289 deallocate(work)
deallocate(wc_a) 1223 1290 deallocate(wc_a)
1224 1291
1292 ! sorting
1293 call fvn_c_sort_eigen(d,evala,eveca)
1294
end subroutine 1225 1295 end subroutine
1226 1296
subroutine fvn_z_matev(d,a,evala,eveca,status) 1227 1297 subroutine fvn_z_matev(d,a,evala,eveca,status)
! 1228 1298 !
! integer d (in) : matrice rank 1229 1299 ! integer d (in) : matrice rank
! double complex a(d,d) (in) : The Matrix 1230 1300 ! double complex a(d,d) (in) : The Matrix
! double complex evala(d) (out) : eigenvalues 1231 1301 ! double complex evala(d) (out) : eigenvalues
! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1232 1302 ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 1233 1303 ! integer (out) : status =0 if something went wrong
! 1234 1304 !
! interfacing Lapack routine ZGEEV 1235 1305 ! interfacing Lapack routine ZGEEV
implicit none 1236 1306 implicit none
integer, intent(in) :: d 1237 1307 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 1238 1308 double complex, intent(in) :: a(d,d)
double complex, intent(out) :: evala(d) 1239 1309 double complex, intent(out) :: evala(d)
double complex, intent(out) :: eveca(d,d) 1240 1310 double complex, intent(out) :: eveca(d,d)
integer, intent(out), optional :: status 1241 1311 integer, intent(out), optional :: status
1242 1312
double complex, allocatable :: wc_a(:,:) ! a working copy of a 1243 1313 double complex, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 1244 1314 integer :: info
integer :: lwork 1245 1315 integer :: lwork
double complex, allocatable :: work(:) 1246 1316 double complex, allocatable :: work(:)
double complex :: twork(1) 1247 1317 double complex :: twork(1)
double precision, allocatable :: rwork(:) 1248 1318 double precision, allocatable :: rwork(:)
double complex :: vl ! unused but necessary for the call 1249 1319 double complex :: vl ! unused but necessary for the call
1250 1320
if (present(status)) status=1 1251 1321 if (present(status)) status=1
1252 1322
! making a working copy of a 1253 1323 ! making a working copy of a
allocate(wc_a(d,d)) 1254 1324 allocate(wc_a(d,d))
!call zcopy(d*d,a,1,wc_a,1) 1255 1325 !call zcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1256 1326 wc_a(:,:)=a(:,:)
1257 1327
! query optimal work size 1258 1328 ! query optimal work size
call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) 1259 1329 call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
lwork=int(twork(1)) 1260 1330 lwork=int(twork(1))
allocate(work(lwork)) 1261 1331 allocate(work(lwork))
allocate(rwork(2*d)) 1262 1332 allocate(rwork(2*d))
call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) 1263 1333 call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
1264 1334
if (info /= 0) then 1265 1335 if (info /= 0) then
if (present(status)) status=0 1266 1336 if (present(status)) status=0
end if 1267 1337 end if
deallocate(rwork) 1268 1338 deallocate(rwork)
deallocate(work) 1269 1339 deallocate(work)
deallocate(wc_a) 1270 1340 deallocate(wc_a)
1341
1342 ! sorting
1343 call fvn_z_sort_eigen(d,evala,eveca)
1271 1344
end subroutine 1272 1345 end subroutine
1273 1346
1274 1347
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1275 1348 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1276 1349 !
! Least square problem 1277 1350 ! Least square problem
! 1278 1351 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1279 1352 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1280 1353 !
! 1281 1354 !
1282 1355
1283 1356
1284 1357
1285 1358
subroutine fvn_d_lspoly(np,x,y,deg,coeff,status) 1286 1359 subroutine fvn_d_lspoly(np,x,y,deg,coeff,status)
! 1287 1360 !
! Least square polynomial fitting 1288 1361 ! Least square polynomial fitting
! 1289 1362 !
! Find the coefficients of the least square polynomial of a given degree 1290 1363 ! Find the coefficients of the least square polynomial of a given degree
! for a set of coordinates. 1291 1364 ! for a set of coordinates.
! 1292 1365 !
! The degree must be lower than the number of points 1293 1366 ! The degree must be lower than the number of points
! 1294 1367 !
! np (in) : number of points 1295 1368 ! np (in) : number of points
! x(np) (in) : x data 1296 1369 ! x(np) (in) : x data
! y(np) (in) : y data 1297 1370 ! y(np) (in) : y data
! deg (in) : polynomial's degree 1298 1371 ! deg (in) : polynomial's degree
! coeff(deg+1) (out) : polynomial coefficients 1299 1372 ! coeff(deg+1) (out) : polynomial coefficients
! status (out) : =0 if a problem occurs 1300 1373 ! status (out) : =0 if a problem occurs
! 1301 1374 !
implicit none 1302 1375 implicit none