lxsd / fvn

1

module fvn

2

module fvn

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

3

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

4

!

! fvn : a f95 module replacement for some imsl routines

5

! fvn : a f95 module replacement for some imsl routines

! it uses lapack for linear algebra

6

! it uses lapack for linear algebra

! it uses modified quadpack for integration

7

! it uses modified quadpack for integration

!

8

!

! William Daniau 2007->today

9

! William Daniau 2007->today

! william.daniau@femto-st.fr

10

! william.daniau@femto-st.fr

!

11

!

! Routines naming scheme :

12

! Routines naming scheme :

!

13

!

! fvn_x_name

14

! fvn_x_name

! where x can be s : real

15

! where x can be s : real

! d : real double precision

16

! d : real double precision

! c : complex

17

! c : complex

! z : double complex

18

! z : double complex

!

19

!

20

!

! This piece of code is totally free! Do whatever you want with it. However

21

! This piece of code is totally free! Do whatever you want with it. However

! if you find it usefull it would be kind to give credits ;-)

22

! if you find it usefull it would be kind to give credits ;-)

!

23

!

! svn version

24

! svn version

! February 2008 : added fnlib to repository so there's no use to have

25

! February 2008 : added fnlib to repository so there's no use to have

! special functions and trigonometry here. Some functions

26

! special functions and trigonometry here. Some functions

! are then removed

27

! are then removed

! January 2008 : added quadratic interpolation, gamma/factorial function,

28

! January 2008 : added quadratic interpolation, gamma/factorial function,

! a function which return identity matrix,

29

! a function which return identity matrix,

! evaluation of nterm chebyshev series

30

! evaluation of nterm chebyshev series

! September 2007 : added sparse system solving by interfacing umfpack

31

! September 2007 : added sparse system solving by interfacing umfpack

! June 2007 : added some complex trigonometric functions

32

! June 2007 : added some complex trigonometric functions

!

33

!

! TO DO LIST :

34

! TO DO LIST :

! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm

35

! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm

! eigenvalues are given with no particular order.

36

! eigenvalues are given with no particular order.

! + Generic interface for fvn_x_name family -> fvn_name

37

! + Generic interface for fvn_x_name family -> fvn_name

! + Make some parameters optional, status for example

38

! + Make some parameters optional, status for example

! + use f95 kinds "double complex" -> complex(kind=8)

39

! + use f95 kinds "double complex" -> complex(kind=8)

! + unify quadpack routines

40

! + unify quadpack routines

! + ...

41

! + ...

!

42

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

43

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

44

implicit none

45

implicit none

! We define pi and i for the module

46

! We define pi and i for the module

real(kind=8),parameter :: fvn_pi = 3.141592653589793_8

47

real(kind=8),parameter :: fvn_pi = 3.141592653589793_8

complex(kind=8),parameter :: fvn_i = (0._8,1._8)

48

complex(kind=8),parameter :: fvn_i = (0._8,1._8)

49

50

!ajout par CW 29/10/07

51

integer :: fvn_status

50

52

! All quadpack routines are private to the module

51

53

! All quadpack routines are private to the module

private :: d1mach,dqag,dqag_2d_inner,dqag_2d_outer,dqage,dqage_2d_inner, &

52

54

private :: d1mach,dqag,dqag_2d_inner,dqag_2d_outer,dqage,dqage_2d_inner, &

dqage_2d_outer,dqk15,dqk15_2d_inner,dqk15_2d_outer,dqk21,dqk21_2d_inner,dqk21_2d_outer, &

53

55

dqage_2d_outer,dqk15,dqk15_2d_inner,dqk15_2d_outer,dqk21,dqk21_2d_inner,dqk21_2d_outer, &

dqk31,dqk31_2d_inner,dqk31_2d_outer,dqk41,dqk41_2d_inner,dqk41_2d_outer, &

54

56

dqk31,dqk31_2d_inner,dqk31_2d_outer,dqk41,dqk41_2d_inner,dqk41_2d_outer, &

dqk51,dqk51_2d_inner,dqk51_2d_outer,dqk61,dqk61_2d_inner,dqk61_2d_outer,dqpsrt

55

57

dqk51,dqk51_2d_inner,dqk51_2d_outer,dqk61,dqk61_2d_inner,dqk61_2d_outer,dqpsrt

56

58

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

57

59

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

58

60

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! Generic interface Definition

59

61

! Generic interface Definition

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

60

62

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

61

63

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

62

64

63

65

! Matrix inversion

64

66

! Matrix inversion

interface fvn_matinv

65

67

interface fvn_matinv

module procedure fvn_s_matinv,fvn_d_matinv,fvn_c_matinv,fvn_z_matinv

66

68

module procedure fvn_s_matinv,fvn_d_matinv,fvn_c_matinv,fvn_z_matinv

end interface fvn_matinv

67

69

end interface fvn_matinv

68

70

! Determinant

69

71

! Determinant

interface fvn_det

70

72

interface fvn_det

module procedure fvn_s_det,fvn_d_det,fvn_c_det,fvn_z_det

71

73

module procedure fvn_s_det,fvn_d_det,fvn_c_det,fvn_z_det

end interface fvn_det

72

74

end interface fvn_det

73

75

! Condition

74

76

! Condition

interface fvn_matcon

75

77

interface fvn_matcon

module procedure fvn_s_matcon,fvn_d_matcon,fvn_c_matcon,fvn_z_matcon

76

78

module procedure fvn_s_matcon,fvn_d_matcon,fvn_c_matcon,fvn_z_matcon

end interface fvn_matcon

77

79

end interface fvn_matcon

78

80

! Eigen

79

81

! Eigen

interface fvn_matev

80

82

interface fvn_matev

module procedure fvn_s_matev,fvn_d_matev,fvn_c_matev,fvn_z_matev

81

83

module procedure fvn_s_matev,fvn_d_matev,fvn_c_matev,fvn_z_matev

end interface fvn_matev

82

84

end interface fvn_matev

83

85

! Utility procedure find interval

84

86

! Utility procedure find interval

interface fvn_find_interval

85

87

interface fvn_find_interval

module procedure fvn_s_find_interval,fvn_d_find_interval

86

88

module procedure fvn_s_find_interval,fvn_d_find_interval

end interface fvn_find_interval

87

89

end interface fvn_find_interval

88

90

! Quadratic 1D interpolation

89

91

! Quadratic 1D interpolation

interface fvn_quad_interpol

90

92

interface fvn_quad_interpol

module procedure fvn_s_quad_interpol,fvn_d_quad_interpol

91

93

module procedure fvn_s_quad_interpol,fvn_d_quad_interpol

end interface fvn_quad_interpol

92

94

end interface fvn_quad_interpol

93

95

! Quadratic 2D interpolation

94

96

! Quadratic 2D interpolation

interface fvn_quad_2d_interpol

95

97

interface fvn_quad_2d_interpol

module procedure fvn_s_quad_2d_interpol,fvn_d_quad_2d_interpol

96

98

module procedure fvn_s_quad_2d_interpol,fvn_d_quad_2d_interpol

end interface fvn_quad_2d_interpol

97

99

end interface fvn_quad_2d_interpol

98

100

! Quadratic 3D interpolation

99

101

! Quadratic 3D interpolation

interface fvn_quad_3d_interpol

100

102

interface fvn_quad_3d_interpol

module procedure fvn_s_quad_3d_interpol,fvn_d_quad_3d_interpol

101

103

module procedure fvn_s_quad_3d_interpol,fvn_d_quad_3d_interpol

end interface fvn_quad_3d_interpol

102

104

end interface fvn_quad_3d_interpol

103

105

! Akima interpolation

104

106

! Akima interpolation

interface fvn_akima

105

107

interface fvn_akima

module procedure fvn_s_akima,fvn_d_akima

106

108

module procedure fvn_s_akima,fvn_d_akima

end interface fvn_akima

107

109

end interface fvn_akima

108

110

! Akima evaluation

109

111

! Akima evaluation

interface fvn_spline_eval

110

112

interface fvn_spline_eval

module procedure fvn_s_spline_eval,fvn_d_spline_eval

111

113

module procedure fvn_s_spline_eval,fvn_d_spline_eval

end interface fvn_spline_eval

112

114

end interface fvn_spline_eval

113

115

! Least square polynomial

114

116

! Least square polynomial

interface fvn_lspoly

115

117

interface fvn_lspoly

module procedure fvn_s_lspoly,fvn_d_lspoly

116

118

module procedure fvn_s_lspoly,fvn_d_lspoly

end interface fvn_lspoly

117

119

end interface fvn_lspoly

118

120

! Muller

119

121

! Muller

interface fvn_muller

120

122

interface fvn_muller

module procedure fvn_z_muller

121

123

module procedure fvn_z_muller

end interface fvn_muller

122

124

end interface fvn_muller

123

125

! Gauss legendre

124

126

! Gauss legendre

interface fvn_gauss_legendre

125

127

interface fvn_gauss_legendre

module procedure fvn_d_gauss_legendre

126

128

module procedure fvn_d_gauss_legendre

end interface fvn_gauss_legendre

127

129

end interface fvn_gauss_legendre

128

130

! Simple Gauss Legendre integration

129

131

! Simple Gauss Legendre integration

interface fvn_gl_integ

130

132

interface fvn_gl_integ

module procedure fvn_d_gl_integ

131

133

module procedure fvn_d_gl_integ

end interface fvn_gl_integ

132

134

end interface fvn_gl_integ

133

135

! Adaptative Gauss Kronrod integration f(x)

134

136

! Adaptative Gauss Kronrod integration f(x)

interface fvn_integ_1_gk

135

137

interface fvn_integ_1_gk

module procedure fvn_d_integ_1_gk

136

138

module procedure fvn_d_integ_1_gk

end interface fvn_integ_1_gk

137

139

end interface fvn_integ_1_gk

138

140

! Adaptative Gauss Kronrod integration f(x,y)

139

141

! Adaptative Gauss Kronrod integration f(x,y)

interface fvn_integ_2_gk

140

142

interface fvn_integ_2_gk

module procedure fvn_d_integ_2_gk

141

143

module procedure fvn_d_integ_2_gk

end interface fvn_integ_2_gk

142

144

end interface fvn_integ_2_gk

143

145

! Sparse solving

144

146

! Sparse solving

interface fvn_sparse_solve

145

147

interface fvn_sparse_solve

module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve

146

148

module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve

end interface fvn_sparse_solve

147

149

end interface fvn_sparse_solve

148

150

149

151

contains

150

152

contains

151

153

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

152

154

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

153

155

!

! Identity Matrix

154

156

! Identity Matrix

!

155

157

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

156

158

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function fvn_d_ident(n)

157

159

function fvn_d_ident(n)

implicit none

158

160

implicit none

integer(kind=4) :: n

159

161

integer(kind=4) :: n

real(kind=8), dimension(n,n) :: fvn_d_ident

160

162

real(kind=8), dimension(n,n) :: fvn_d_ident

161

163

real(kind=8),dimension(n*n) :: vect

162

164

real(kind=8),dimension(n*n) :: vect

integer(kind=4) :: i

163

165

integer(kind=4) :: i

164

166

vect=0._8

165

167

vect=0._8

vect(1:n*n:n+1) = 1._8

166

168

vect(1:n*n:n+1) = 1._8

fvn_d_ident=reshape(vect, shape = (/ n,n /))

167

169

fvn_d_ident=reshape(vect, shape = (/ n,n /))

end function

168

170

end function

169

171

function fvn_s_ident(n)

170

172

function fvn_s_ident(n)

implicit none

171

173

implicit none

integer(kind=4) :: n

172

174

integer(kind=4) :: n

real(kind=4), dimension(n,n) :: fvn_s_ident

173

175

real(kind=4), dimension(n,n) :: fvn_s_ident

174

176

real(kind=4),dimension(n*n) :: vect

175

177

real(kind=4),dimension(n*n) :: vect

integer(kind=4) :: i

176

178

integer(kind=4) :: i

177

179

vect=0._4

178

180

vect=0._4

vect(1:n*n:n+1) = 1._4

179

181

vect(1:n*n:n+1) = 1._4

fvn_s_ident=reshape(vect, shape = (/ n,n /))

180

182

fvn_s_ident=reshape(vect, shape = (/ n,n /))

end function

181

183

end function

182

184

function fvn_c_ident(n)

183

185

function fvn_c_ident(n)

implicit none

184

186

implicit none

integer(kind=4) :: n

185

187

integer(kind=4) :: n

complex(kind=4), dimension(n,n) :: fvn_c_ident

186

188

complex(kind=4), dimension(n,n) :: fvn_c_ident

187

189

complex(kind=4),dimension(n*n) :: vect

188

190

complex(kind=4),dimension(n*n) :: vect

integer(kind=4) :: i

189

191

integer(kind=4) :: i

190

192

vect=(0._4,0._4)

191

193

vect=(0._4,0._4)

vect(1:n*n:n+1) = (1._4,0._4)

192

194

vect(1:n*n:n+1) = (1._4,0._4)

fvn_c_ident=reshape(vect, shape = (/ n,n /))

193

195

fvn_c_ident=reshape(vect, shape = (/ n,n /))

end function

194

196

end function

195

197

function fvn_z_ident(n)

196

198

function fvn_z_ident(n)

implicit none

197

199

implicit none

integer(kind=4) :: n

198

200

integer(kind=4) :: n

complex(kind=8), dimension(n,n) :: fvn_z_ident

199

201

complex(kind=8), dimension(n,n) :: fvn_z_ident

200

202

complex(kind=8),dimension(n*n) :: vect

201

203

complex(kind=8),dimension(n*n) :: vect

integer(kind=4) :: i

202

204

integer(kind=4) :: i

203

205

vect=(0._8,0._8)

204

206

vect=(0._8,0._8)

vect(1:n*n:n+1) = (1._8,0._8)

205

207

vect(1:n*n:n+1) = (1._8,0._8)

fvn_z_ident=reshape(vect, shape = (/ n,n /))

206

208

fvn_z_ident=reshape(vect, shape = (/ n,n /))

end function

207

209

end function

208

210

209

211

210

212

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

211

213

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

212

214

!

! Matrix inversion subroutines

213

215

! Matrix inversion subroutines

!

214

216

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

215

217

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

216

218

subroutine fvn_s_matinv(d,a,inva,status)

217

219

subroutine fvn_s_matinv(d,a,inva,status)

!

218

220

!

! Matrix inversion of a real matrix using BLAS and LAPACK

219

221

! Matrix inversion of a real matrix using BLAS and LAPACK

!

220

222

!

! d (in) : matrix rank

221

223

! d (in) : matrix rank

! a (in) : input matrix

222

224

! a (in) : input matrix

! inva (out) : inversed matrix

223

225

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

224

226

! status (ou) : =0 if something failed

!

225

227

!

implicit none

226

228

implicit none

integer, intent(in) :: d

227

229

integer, intent(in) :: d

real, intent(in) :: a(d,d)

228

230

real, intent(in) :: a(d,d)

real, intent(out) :: inva(d,d)

229

231

real, intent(out) :: inva(d,d)

integer, intent(out),optional :: status

230

232

integer, intent(out),optional :: status

231

233

integer, allocatable :: ipiv(:)

232

234

integer, allocatable :: ipiv(:)

real, allocatable :: work(:)

233

235

real, allocatable :: work(:)

real twork(1)

234

236

real twork(1)

integer :: info

235

237

integer :: info

integer :: lwork

236

238

integer :: lwork

237

239

if (present(status)) status=1

238

240

if (present(status)) status=1

239

241

allocate(ipiv(d))

240

242

allocate(ipiv(d))

! copy a into inva using BLAS

241

243

! copy a into inva using BLAS

!call scopy(d*d,a,1,inva,1)

242

244

!call scopy(d*d,a,1,inva,1)

inva(:,:)=a(:,:)

243

245

inva(:,:)=a(:,:)

! LU factorization using LAPACK

244

246

! LU factorization using LAPACK

call sgetrf(d,d,inva,d,ipiv,info)

245

247

call sgetrf(d,d,inva,d,ipiv,info)

! if info is not equal to 0, something went wrong we exit setting status to 0

246

248

! if info is not equal to 0, something went wrong we exit setting status to 0

if (info /= 0) then

247

249

if (info /= 0) then

if (present(status)) status=0

248

250

if (present(status)) status=0

deallocate(ipiv)

249

251

deallocate(ipiv)

return

250

252

return

end if

251

253

end if

! we use the query fonction of xxxtri to obtain the optimal workspace size

252

254

! we use the query fonction of xxxtri to obtain the optimal workspace size

call sgetri(d,inva,d,ipiv,twork,-1,info)

253

255

call sgetri(d,inva,d,ipiv,twork,-1,info)

lwork=int(twork(1))

254

256

lwork=int(twork(1))

allocate(work(lwork))

255

257

allocate(work(lwork))

! Matrix inversion using LAPACK

256

258

! Matrix inversion using LAPACK

call sgetri(d,inva,d,ipiv,work,lwork,info)

257

259

call sgetri(d,inva,d,ipiv,work,lwork,info)

! again if info is not equal to 0, we exit setting status to 0

258

260

! again if info is not equal to 0, we exit setting status to 0

if (info /= 0) then

259

261

if (info /= 0) then

if (present(status)) status=0

260

262

if (present(status)) status=0

end if

261

263

end if

deallocate(work)

262

264

deallocate(work)

deallocate(ipiv)

263

265

deallocate(ipiv)

end subroutine

264

266

end subroutine

265

267

subroutine fvn_d_matinv(d,a,inva,status)

266

268

subroutine fvn_d_matinv(d,a,inva,status)

!

267

269

!

! Matrix inversion of a double precision matrix using BLAS and LAPACK

268

270

! Matrix inversion of a double precision matrix using BLAS and LAPACK

!

269

271

!

! d (in) : matrix rank

270

272

! d (in) : matrix rank

! a (in) : input matrix

271

273

! a (in) : input matrix

! inva (out) : inversed matrix

272

274

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

273

275

! status (ou) : =0 if something failed

!

274

276

!

implicit none

275

277

implicit none

integer, intent(in) :: d

276

278

integer, intent(in) :: d

double precision, intent(in) :: a(d,d)

277

279

double precision, intent(in) :: a(d,d)

double precision, intent(out) :: inva(d,d)

278

280

double precision, intent(out) :: inva(d,d)

integer, intent(out),optional :: status

279

281

integer, intent(out),optional :: status

280

282

integer, allocatable :: ipiv(:)

281

283

integer, allocatable :: ipiv(:)

double precision, allocatable :: work(:)

282

284

double precision, allocatable :: work(:)

double precision :: twork(1)

283

285

double precision :: twork(1)

integer :: info

284

286

integer :: info

integer :: lwork

285

287

integer :: lwork

286

288

if (present(status)) status=1

287

289

if (present(status)) status=1

288

290

allocate(ipiv(d))

289

291

allocate(ipiv(d))

! copy a into inva using BLAS

290

292

! copy a into inva using BLAS

!call dcopy(d*d,a,1,inva,1)

291

293

!call dcopy(d*d,a,1,inva,1)

inva(:,:)=a(:,:)

292

294

inva(:,:)=a(:,:)

! LU factorization using LAPACK

293

295

! LU factorization using LAPACK

call dgetrf(d,d,inva,d,ipiv,info)

294

296

call dgetrf(d,d,inva,d,ipiv,info)

! if info is not equal to 0, something went wrong we exit setting status to 0

295

297

! if info is not equal to 0, something went wrong we exit setting status to 0

if (info /= 0) then

296

298

if (info /= 0) then

if (present(status)) status=0

297

299

if (present(status)) status=0

deallocate(ipiv)

298

300

deallocate(ipiv)

return

299

301

return

end if

300

302

end if

! we use the query fonction of xxxtri to obtain the optimal workspace size

301

303

! we use the query fonction of xxxtri to obtain the optimal workspace size

call dgetri(d,inva,d,ipiv,twork,-1,info)

302

304

call dgetri(d,inva,d,ipiv,twork,-1,info)

lwork=int(twork(1))

303

305

lwork=int(twork(1))

allocate(work(lwork))

304

306

allocate(work(lwork))

! Matrix inversion using LAPACK

305

307

! Matrix inversion using LAPACK

call dgetri(d,inva,d,ipiv,work,lwork,info)

306

308

call dgetri(d,inva,d,ipiv,work,lwork,info)

! again if info is not equal to 0, we exit setting status to 0

307

309

! again if info is not equal to 0, we exit setting status to 0

if (info /= 0) then

308

310

if (info /= 0) then

if (present(status)) status=0

309

311

if (present(status)) status=0

end if

310

312

end if

deallocate(work)

311

313

deallocate(work)

deallocate(ipiv)

312

314

deallocate(ipiv)

end subroutine

313

315

end subroutine

314

316

subroutine fvn_c_matinv(d,a,inva,status)

315

317

subroutine fvn_c_matinv(d,a,inva,status)

!

316

318

!

! Matrix inversion of a complex matrix using BLAS and LAPACK

317

319

! Matrix inversion of a complex matrix using BLAS and LAPACK

!

318

320

!

! d (in) : matrix rank

319

321

! d (in) : matrix rank

! a (in) : input matrix

320

322

! a (in) : input matrix

! inva (out) : inversed matrix

321

323

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

322

324

! status (ou) : =0 if something failed

!

323

325

!

implicit none

324

326

implicit none

integer, intent(in) :: d

325

327

integer, intent(in) :: d

complex, intent(in) :: a(d,d)

326

328

complex, intent(in) :: a(d,d)

complex, intent(out) :: inva(d,d)

327

329

complex, intent(out) :: inva(d,d)

integer, intent(out),optional :: status

328

330

integer, intent(out),optional :: status

329

331

integer, allocatable :: ipiv(:)

330

332

integer, allocatable :: ipiv(:)

complex, allocatable :: work(:)

331

333

complex, allocatable :: work(:)

complex :: twork(1)

332

334

complex :: twork(1)

integer :: info

333

335

integer :: info

integer :: lwork

334

336

integer :: lwork

335

337

if (present(status)) status=1

336

338

if (present(status)) status=1

337

339

allocate(ipiv(d))

338

340

allocate(ipiv(d))

! copy a into inva using BLAS

339

341

! copy a into inva using BLAS

!call ccopy(d*d,a,1,inva,1)

340

342

!call ccopy(d*d,a,1,inva,1)

inva(:,:)=a(:,:)

341

343

inva(:,:)=a(:,:)

342

344

! LU factorization using LAPACK

343

345

! LU factorization using LAPACK

call cgetrf(d,d,inva,d,ipiv,info)

344

346

call cgetrf(d,d,inva,d,ipiv,info)

! if info is not equal to 0, something went wrong we exit setting status to 0

345

347

! if info is not equal to 0, something went wrong we exit setting status to 0

if (info /= 0) then

346

348

if (info /= 0) then

if (present(status)) status=0

347

349

if (present(status)) status=0

deallocate(ipiv)

348

350

deallocate(ipiv)

return

349

351

return

end if

350

352

end if

! we use the query fonction of xxxtri to obtain the optimal workspace size

351

353

! we use the query fonction of xxxtri to obtain the optimal workspace size

call cgetri(d,inva,d,ipiv,twork,-1,info)

352

354

call cgetri(d,inva,d,ipiv,twork,-1,info)

lwork=int(twork(1))

353

355

lwork=int(twork(1))

allocate(work(lwork))

354

356

allocate(work(lwork))

! Matrix inversion using LAPACK

355

357

! Matrix inversion using LAPACK

call cgetri(d,inva,d,ipiv,work,lwork,info)

356

358

call cgetri(d,inva,d,ipiv,work,lwork,info)

! again if info is not equal to 0, we exit setting status to 0

357

359

! again if info is not equal to 0, we exit setting status to 0

if (info /= 0) then

358

360

if (info /= 0) then

if (present(status)) status=0

359

361

if (present(status)) status=0

end if

360

362

end if

deallocate(work)

361

363

deallocate(work)

deallocate(ipiv)

362

364

deallocate(ipiv)

end subroutine

363

365

end subroutine

364

366

subroutine fvn_z_matinv(d,a,inva,status)

365

367

subroutine fvn_z_matinv(d,a,inva,status)

!

366

368

!

! Matrix inversion of a double complex matrix using BLAS and LAPACK

367

369

! Matrix inversion of a double complex matrix using BLAS and LAPACK

!

368

370

!

! d (in) : matrix rank

369

371

! d (in) : matrix rank

! a (in) : input matrix

370

372

! a (in) : input matrix

! inva (out) : inversed matrix

371

373

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

372

374

! status (ou) : =0 if something failed

!

373

375

!

implicit none

374

376

implicit none

integer, intent(in) :: d

375

377

integer, intent(in) :: d

double complex, intent(in) :: a(d,d)

376

378

double complex, intent(in) :: a(d,d)

double complex, intent(out) :: inva(d,d)

377

379

double complex, intent(out) :: inva(d,d)

integer, intent(out),optional :: status

378

380

integer, intent(out),optional :: status

379

381

integer, allocatable :: ipiv(:)

380

382

integer, allocatable :: ipiv(:)

double complex, allocatable :: work(:)

381

383

double complex, allocatable :: work(:)

double complex :: twork(1)

382

384

double complex :: twork(1)

integer :: info

383

385

integer :: info

integer :: lwork

384

386

integer :: lwork

385

387

if (present(status)) status=1

386

388

if (present(status)) status=1

387

389

allocate(ipiv(d))

388

390

allocate(ipiv(d))

! copy a into inva using BLAS

389

391

! copy a into inva using BLAS

!call zcopy(d*d,a,1,inva,1)

390

392

!call zcopy(d*d,a,1,inva,1)

inva(:,:)=a(:,:)

391

393

inva(:,:)=a(:,:)

392

394

! LU factorization using LAPACK

393

395

! LU factorization using LAPACK

call zgetrf(d,d,inva,d,ipiv,info)

394

396

call zgetrf(d,d,inva,d,ipiv,info)

! if info is not equal to 0, something went wrong we exit setting status to 0

395

397

! if info is not equal to 0, something went wrong we exit setting status to 0

if (info /= 0) then

396

398

if (info /= 0) then

if (present(status)) status=0

397

399

if (present(status)) status=0

deallocate(ipiv)

398

400

deallocate(ipiv)

return

399

401

return

end if

400

402

end if

! we use the query fonction of xxxtri to obtain the optimal workspace size

401

403

! we use the query fonction of xxxtri to obtain the optimal workspace size

call zgetri(d,inva,d,ipiv,twork,-1,info)

402

404

call zgetri(d,inva,d,ipiv,twork,-1,info)

lwork=int(twork(1))

403

405

lwork=int(twork(1))

allocate(work(lwork))

404

406

allocate(work(lwork))

! Matrix inversion using LAPACK

405

407

! Matrix inversion using LAPACK

call zgetri(d,inva,d,ipiv,work,lwork,info)

406

408

call zgetri(d,inva,d,ipiv,work,lwork,info)

! again if info is not equal to 0, we exit setting status to 0

407

409

! again if info is not equal to 0, we exit setting status to 0

if (info /= 0) then

408

410

if (info /= 0) then

if (present(status)) status=0

409

411

if (present(status)) status=0

end if

410

412

end if

deallocate(work)

411

413

deallocate(work)

deallocate(ipiv)

412

414

deallocate(ipiv)

end subroutine

413

415

end subroutine

414

416

415

417

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

416

418

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

417

419

!

! Determinants

418

420

! Determinants

!

419

421

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

420

422

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function fvn_s_det(d,a,status)

421

423

function fvn_s_det(d,a,status)

!

422

424

!

! Evaluate the determinant of a square matrix using lapack LU factorization

423

425

! Evaluate the determinant of a square matrix using lapack LU factorization

!

424

426

!

! d (in) : matrix rank

425

427

! d (in) : matrix rank

! a (in) : The Matrix

426

428

! a (in) : The Matrix

! status (out) : =0 if LU factorization failed

427

429

! status (out) : =0 if LU factorization failed

!

428

430

!

implicit none

429

431

implicit none

integer, intent(in) :: d

430

432

integer, intent(in) :: d

real, intent(in) :: a(d,d)

431

433

real, intent(in) :: a(d,d)

integer, intent(out), optional :: status

432

434

integer, intent(out), optional :: status

real :: fvn_s_det

433

435

real :: fvn_s_det

434

436

real, allocatable :: wc_a(:,:)

435

437

real, allocatable :: wc_a(:,:)

integer, allocatable :: ipiv(:)

436

438

integer, allocatable :: ipiv(:)

integer :: info,i

437

439

integer :: info,i

438

440

if (present(status)) status=1

439

441

if (present(status)) status=1

allocate(wc_a(d,d))

440

442

allocate(wc_a(d,d))

allocate(ipiv(d))

441

443

allocate(ipiv(d))

wc_a(:,:)=a(:,:)

442

444

wc_a(:,:)=a(:,:)

call sgetrf(d,d,wc_a,d,ipiv,info)

443

445

call sgetrf(d,d,wc_a,d,ipiv,info)

if (info/= 0) then

444

446

if (info/= 0) then

if (present(status)) status=0

445

447

if (present(status)) status=0

fvn_s_det=0.e0

446

448

fvn_s_det=0.e0

deallocate(ipiv)

447

449

deallocate(ipiv)

deallocate(wc_a)

448

450

deallocate(wc_a)

return

449

451

return

end if

450

452

end if

fvn_s_det=1.e0

451

453

fvn_s_det=1.e0

do i=1,d

452

454

do i=1,d

if (ipiv(i)==i) then

453

455

if (ipiv(i)==i) then

fvn_s_det=fvn_s_det*wc_a(i,i)

454

456

fvn_s_det=fvn_s_det*wc_a(i,i)

else

455

457

else

fvn_s_det=-fvn_s_det*wc_a(i,i)

456

458

fvn_s_det=-fvn_s_det*wc_a(i,i)

end if

457

459

end if

end do

458

460

end do

deallocate(ipiv)

459

461

deallocate(ipiv)

deallocate(wc_a)

460

462

deallocate(wc_a)

461

463

end function

462

464

end function

463

465

function fvn_d_det(d,a,status)

464

466

function fvn_d_det(d,a,status)

!

465

467

!

! Evaluate the determinant of a square matrix using lapack LU factorization

466

468

! Evaluate the determinant of a square matrix using lapack LU factorization

!

467

469

!

! d (in) : matrix rank

468

470

! d (in) : matrix rank

! a (in) : The Matrix

469

471

! a (in) : The Matrix

! status (out) : =0 if LU factorization failed

470

472

! status (out) : =0 if LU factorization failed

!

471

473

!

implicit none

472

474

implicit none

integer, intent(in) :: d

473

475

integer, intent(in) :: d

double precision, intent(in) :: a(d,d)

474

476

double precision, intent(in) :: a(d,d)

integer, intent(out), optional :: status

475

477

integer, intent(out), optional :: status

double precision :: fvn_d_det

476

478

double precision :: fvn_d_det

477

479

double precision, allocatable :: wc_a(:,:)

478

480

double precision, allocatable :: wc_a(:,:)

integer, allocatable :: ipiv(:)

479

481

integer, allocatable :: ipiv(:)

integer :: info,i

480

482

integer :: info,i

481

483

if (present(status)) status=1

482

484

if (present(status)) status=1

allocate(wc_a(d,d))

483

485

allocate(wc_a(d,d))

allocate(ipiv(d))

484

486

allocate(ipiv(d))

wc_a(:,:)=a(:,:)

485

487

wc_a(:,:)=a(:,:)

call dgetrf(d,d,wc_a,d,ipiv,info)

486

488

call dgetrf(d,d,wc_a,d,ipiv,info)

if (info/= 0) then

487

489

if (info/= 0) then

if (present(status)) status=0

488

490

if (present(status)) status=0

fvn_d_det=0.d0

489

491

fvn_d_det=0.d0

deallocate(ipiv)

490

492

deallocate(ipiv)

deallocate(wc_a)

491

493

deallocate(wc_a)

return

492

494

return

end if

493

495

end if

fvn_d_det=1.d0

494

496

fvn_d_det=1.d0

do i=1,d

495

497

do i=1,d

if (ipiv(i)==i) then

496

498

if (ipiv(i)==i) then

fvn_d_det=fvn_d_det*wc_a(i,i)

497

499

fvn_d_det=fvn_d_det*wc_a(i,i)

else

498

500

else

fvn_d_det=-fvn_d_det*wc_a(i,i)

499

501

fvn_d_det=-fvn_d_det*wc_a(i,i)

end if

500

502

end if

end do

501

503

end do

deallocate(ipiv)

502

504

deallocate(ipiv)

deallocate(wc_a)

503

505

deallocate(wc_a)

504

506

end function

505

507

end function

506

508

function fvn_c_det(d,a,status) !

507

509

function fvn_c_det(d,a,status) !

! Evaluate the determinant of a square matrix using lapack LU factorization

508

510

! Evaluate the determinant of a square matrix using lapack LU factorization

!

509

511

!

! d (in) : matrix rank

510

512

! d (in) : matrix rank

! a (in) : The Matrix

511

513

! a (in) : The Matrix

! status (out) : =0 if LU factorization failed

512

514

! status (out) : =0 if LU factorization failed

!

513

515

!

implicit none

514

516

implicit none

integer, intent(in) :: d

515

517

integer, intent(in) :: d

complex, intent(in) :: a(d,d)

516

518

complex, intent(in) :: a(d,d)

integer, intent(out), optional :: status

517

519

integer, intent(out), optional :: status

complex :: fvn_c_det

518

520

complex :: fvn_c_det

519

521

complex, allocatable :: wc_a(:,:)

520

522

complex, allocatable :: wc_a(:,:)

integer, allocatable :: ipiv(:)

521

523

integer, allocatable :: ipiv(:)

integer :: info,i

522

524

integer :: info,i

523

525

if (present(status)) status=1

524

526

if (present(status)) status=1

allocate(wc_a(d,d))

525

527

allocate(wc_a(d,d))

allocate(ipiv(d))

526

528

allocate(ipiv(d))

wc_a(:,:)=a(:,:)

527

529

wc_a(:,:)=a(:,:)

call cgetrf(d,d,wc_a,d,ipiv,info)

528

530

call cgetrf(d,d,wc_a,d,ipiv,info)

if (info/= 0) then

529

531

if (info/= 0) then

if (present(status)) status=0

530

532

if (present(status)) status=0

fvn_c_det=(0.e0,0.e0)

531

533

fvn_c_det=(0.e0,0.e0)

deallocate(ipiv)

532

534

deallocate(ipiv)

deallocate(wc_a)

533

535

deallocate(wc_a)

return

534

536

return

end if

535

537

end if

fvn_c_det=(1.e0,0.e0)

536

538

fvn_c_det=(1.e0,0.e0)

do i=1,d

537

539

do i=1,d

if (ipiv(i)==i) then

538

540

if (ipiv(i)==i) then

fvn_c_det=fvn_c_det*wc_a(i,i)

539

541

fvn_c_det=fvn_c_det*wc_a(i,i)

else

540

542

else

fvn_c_det=-fvn_c_det*wc_a(i,i)

541

543

fvn_c_det=-fvn_c_det*wc_a(i,i)

end if

542

544

end if

end do

543

545

end do

deallocate(ipiv)

544

546

deallocate(ipiv)

deallocate(wc_a)

545

547

deallocate(wc_a)

546

548

end function

547

549

end function

548

550

function fvn_z_det(d,a,status)

549

551

function fvn_z_det(d,a,status)

!

550

552

!

! Evaluate the determinant of a square matrix using lapack LU factorization

551

553

! Evaluate the determinant of a square matrix using lapack LU factorization

!

552

554

!

! d (in) : matrix rank

553

555

! d (in) : matrix rank

! a (in) : The Matrix

554

556

! a (in) : The Matrix

! det (out) : determinant

555

557

! det (out) : determinant

! status (out) : =0 if LU factorization failed

556

558

! status (out) : =0 if LU factorization failed

!

557

559

!

implicit none

558

560

implicit none

integer, intent(in) :: d

559

561

integer, intent(in) :: d

double complex, intent(in) :: a(d,d)

560

562

double complex, intent(in) :: a(d,d)

integer, intent(out), optional :: status

561

563

integer, intent(out), optional :: status

double complex :: fvn_z_det

562

564

double complex :: fvn_z_det

563

565

double complex, allocatable :: wc_a(:,:)

564

566

double complex, allocatable :: wc_a(:,:)

integer, allocatable :: ipiv(:)

565

567

integer, allocatable :: ipiv(:)

integer :: info,i

566

568

integer :: info,i

567

569

if (present(status)) status=1

568

570

if (present(status)) status=1

allocate(wc_a(d,d))

569

571

allocate(wc_a(d,d))

allocate(ipiv(d))

570

572

allocate(ipiv(d))

wc_a(:,:)=a(:,:)

571

573

wc_a(:,:)=a(:,:)

call zgetrf(d,d,wc_a,d,ipiv,info)

572

574

call zgetrf(d,d,wc_a,d,ipiv,info)

if (info/= 0) then

573

575

if (info/= 0) then

if (present(status)) status=0

574

576

if (present(status)) status=0

fvn_z_det=(0.d0,0.d0)

575

577

fvn_z_det=(0.d0,0.d0)

deallocate(ipiv)

576

578

deallocate(ipiv)

deallocate(wc_a)

577

579

deallocate(wc_a)

return

578

580

return

end if

579

581

end if

fvn_z_det=(1.d0,0.d0)

580

582

fvn_z_det=(1.d0,0.d0)

do i=1,d

581

583

do i=1,d

if (ipiv(i)==i) then

582

584

if (ipiv(i)==i) then

fvn_z_det=fvn_z_det*wc_a(i,i)

583

585

fvn_z_det=fvn_z_det*wc_a(i,i)

else

584

586

else

fvn_z_det=-fvn_z_det*wc_a(i,i)

585

587

fvn_z_det=-fvn_z_det*wc_a(i,i)

end if

586

588

end if

end do

587

589

end do

deallocate(ipiv)

588

590

deallocate(ipiv)

deallocate(wc_a)

589

591

deallocate(wc_a)

590

592

end function

591

593

end function

592

594

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

593

595

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

594

596

!

! Condition test

595

597

! Condition test

!

596

598

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

597

599

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! 1-norm

598

600

! 1-norm

! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm

599

601

! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm

! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond

600

602

! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond

!

601

603

!

subroutine fvn_s_matcon(d,a,rcond,status)

602

604

subroutine fvn_s_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

603

605

! Matrix condition (reciprocal of condition number)

!

604

606

!

! d (in) : matrix rank

605

607

! d (in) : matrix rank

! a (in) : The Matrix

606

608

! a (in) : The Matrix

! rcond (out) : guess what

607

609

! rcond (out) : guess what

! status (out) : =0 if something went wrong

608

610

! status (out) : =0 if something went wrong

!

609

611

!

implicit none

610

612

implicit none

integer, intent(in) :: d

611

613

integer, intent(in) :: d

real, intent(in) :: a(d,d)

612

614

real, intent(in) :: a(d,d)

real, intent(out) :: rcond

613

615

real, intent(out) :: rcond

integer, intent(out), optional :: status

614

616

integer, intent(out), optional :: status

615

617

real, allocatable :: work(:)

616

618

real, allocatable :: work(:)

integer, allocatable :: iwork(:)

617

619

integer, allocatable :: iwork(:)

real :: anorm

618

620

real :: anorm

real, allocatable :: wc_a(:,:) ! working copy of a

619

621

real, allocatable :: wc_a(:,:) ! working copy of a

integer :: info

620

622

integer :: info

integer, allocatable :: ipiv(:)

621

623

integer, allocatable :: ipiv(:)

622

624

real, external :: slange

623

625

real, external :: slange

624

626

625

627

if (present(status)) status=1

626

628

if (present(status)) status=1

627

629

anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm

628

630

anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm

629

631

allocate(wc_a(d,d))

630

632

allocate(wc_a(d,d))

!call scopy(d*d,a,1,wc_a,1)

631

633

!call scopy(d*d,a,1,wc_a,1)

wc_a(:,:)=a(:,:)

632

634

wc_a(:,:)=a(:,:)

633

635

allocate(ipiv(d))

634

636

allocate(ipiv(d))

call sgetrf(d,d,wc_a,d,ipiv,info)

635

637

call sgetrf(d,d,wc_a,d,ipiv,info)

if (info /= 0) then

636

638

if (info /= 0) then

if (present(status)) status=0

637

639

if (present(status)) status=0

deallocate(ipiv)

638

640

deallocate(ipiv)

deallocate(wc_a)

639

641

deallocate(wc_a)

return

640

642

return

end if

641

643

end if

allocate(work(4*d))

642

644

allocate(work(4*d))

allocate(iwork(d))

643

645

allocate(iwork(d))

call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)

644

646

call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)

if (info /= 0) then

645

647

if (info /= 0) then

if (present(status)) status=0

646

648

if (present(status)) status=0

end if

647

649

end if

deallocate(iwork)

648

650

deallocate(iwork)

deallocate(work)

649

651

deallocate(work)

deallocate(ipiv)

650

652

deallocate(ipiv)

deallocate(wc_a)

651

653

deallocate(wc_a)

652

654

end subroutine

653

655

end subroutine

654

656

subroutine fvn_d_matcon(d,a,rcond,status)

655

657

subroutine fvn_d_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

656

658

! Matrix condition (reciprocal of condition number)

!

657

659

!

! d (in) : matrix rank

658

660

! d (in) : matrix rank

! a (in) : The Matrix

659

661

! a (in) : The Matrix

! rcond (out) : guess what

660

662

! rcond (out) : guess what

! status (out) : =0 if something went wrong

661

663

! status (out) : =0 if something went wrong

!

662

664

!

implicit none

663

665

implicit none

integer, intent(in) :: d

664

666

integer, intent(in) :: d

double precision, intent(in) :: a(d,d)

665

667

double precision, intent(in) :: a(d,d)

double precision, intent(out) :: rcond

666

668

double precision, intent(out) :: rcond

integer, intent(out), optional :: status

667

669

integer, intent(out), optional :: status

668

670

double precision, allocatable :: work(:)

669

671

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

670

672

integer, allocatable :: iwork(:)

double precision :: anorm

671

673

double precision :: anorm

double precision, allocatable :: wc_a(:,:) ! working copy of a

672

674

double precision, allocatable :: wc_a(:,:) ! working copy of a

integer :: info

673

675

integer :: info

integer, allocatable :: ipiv(:)

674

676

integer, allocatable :: ipiv(:)

675

677

double precision, external :: dlange

676

678

double precision, external :: dlange

677

679

678

680

if (present(status)) status=1

679

681

if (present(status)) status=1

680

682

anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm

681

683

anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm

682

684

allocate(wc_a(d,d))

683

685

allocate(wc_a(d,d))

!call dcopy(d*d,a,1,wc_a,1)

684

686

!call dcopy(d*d,a,1,wc_a,1)

wc_a(:,:)=a(:,:)

685

687

wc_a(:,:)=a(:,:)

686

688

allocate(ipiv(d))

687

689

allocate(ipiv(d))

call dgetrf(d,d,wc_a,d,ipiv,info)

688

690

call dgetrf(d,d,wc_a,d,ipiv,info)

if (info /= 0) then

689

691

if (info /= 0) then

if (present(status)) status=0

690

692

if (present(status)) status=0

deallocate(ipiv)

691

693

deallocate(ipiv)

deallocate(wc_a)

692

694

deallocate(wc_a)

return

693

695

return

end if

694

696

end if

695

697

allocate(work(4*d))

696

698

allocate(work(4*d))

allocate(iwork(d))

697

699

allocate(iwork(d))

call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)

698

700

call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)

if (info /= 0) then

699

701

if (info /= 0) then

if (present(status)) status=0

700

702

if (present(status)) status=0

end if

701

703

end if

deallocate(iwork)

702

704

deallocate(iwork)

deallocate(work)

703

705

deallocate(work)

deallocate(ipiv)

704

706

deallocate(ipiv)

deallocate(wc_a)

705

707

deallocate(wc_a)

706

708

end subroutine

707

709

end subroutine

708

710

subroutine fvn_c_matcon(d,a,rcond,status)

709

711

subroutine fvn_c_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

710

712

! Matrix condition (reciprocal of condition number)

!

711

713

!

! d (in) : matrix rank

712

714

! d (in) : matrix rank

! a (in) : The Matrix

713

715

! a (in) : The Matrix

! rcond (out) : guess what

714

716

! rcond (out) : guess what

! status (out) : =0 if something went wrong

715

717

! status (out) : =0 if something went wrong

!

716

718

!

implicit none

717

719

implicit none

integer, intent(in) :: d

718

720

integer, intent(in) :: d

complex, intent(in) :: a(d,d)

719

721

complex, intent(in) :: a(d,d)

real, intent(out) :: rcond

720

722

real, intent(out) :: rcond

integer, intent(out), optional :: status

721

723

integer, intent(out), optional :: status

722

724

real, allocatable :: rwork(:)

723

725

real, allocatable :: rwork(:)

complex, allocatable :: work(:)

724

726

complex, allocatable :: work(:)

integer, allocatable :: iwork(:)

725

727

integer, allocatable :: iwork(:)

real :: anorm

726

728

real :: anorm

complex, allocatable :: wc_a(:,:) ! working copy of a

727

729

complex, allocatable :: wc_a(:,:) ! working copy of a

integer :: info

728

730

integer :: info

integer, allocatable :: ipiv(:)

729

731

integer, allocatable :: ipiv(:)

730

732

real, external :: clange

731

733

real, external :: clange

732

734

733

735

if (present(status)) status=1

734

736

if (present(status)) status=1

735

737

anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm

736

738

anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm

737

739

allocate(wc_a(d,d))

738

740

allocate(wc_a(d,d))

!call ccopy(d*d,a,1,wc_a,1)

739

741

!call ccopy(d*d,a,1,wc_a,1)

wc_a(:,:)=a(:,:)

740

742

wc_a(:,:)=a(:,:)

741

743

allocate(ipiv(d))

742

744

allocate(ipiv(d))

call cgetrf(d,d,wc_a,d,ipiv,info)

743

745

call cgetrf(d,d,wc_a,d,ipiv,info)

if (info /= 0) then

744

746

if (info /= 0) then

if (present(status)) status=0

745

747

if (present(status)) status=0

deallocate(ipiv)

746

748

deallocate(ipiv)

deallocate(wc_a)

747

749

deallocate(wc_a)

return

748

750

return

end if

749

751

end if

allocate(work(2*d))

750

752

allocate(work(2*d))

allocate(rwork(2*d))

751

753

allocate(rwork(2*d))

call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)

752

754

call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)

if (info /= 0) then

753

755

if (info /= 0) then

if (present(status)) status=0

754

756

if (present(status)) status=0

end if

755

757

end if

deallocate(rwork)

756

758

deallocate(rwork)

deallocate(work)

757

759

deallocate(work)

deallocate(ipiv)

758

760

deallocate(ipiv)

deallocate(wc_a)

759

761

deallocate(wc_a)

end subroutine

760

762

end subroutine

761

763

subroutine fvn_z_matcon(d,a,rcond,status)

762

764

subroutine fvn_z_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

763

765

! Matrix condition (reciprocal of condition number)

!

764

766

!

! d (in) : matrix rank

765

767

! d (in) : matrix rank

! a (in) : The Matrix

766

768

! a (in) : The Matrix

! rcond (out) : guess what

767

769

! rcond (out) : guess what

! status (out) : =0 if something went wrong

768

770

! status (out) : =0 if something went wrong

!

769

771

!

implicit none

770

772

implicit none

integer, intent(in) :: d

771

773

integer, intent(in) :: d

double complex, intent(in) :: a(d,d)

772

774

double complex, intent(in) :: a(d,d)

double precision, intent(out) :: rcond

773

775

double precision, intent(out) :: rcond

integer, intent(out), optional :: status

774

776

integer, intent(out), optional :: status

775

777

double complex, allocatable :: work(:)

776

778

double complex, allocatable :: work(:)

double precision, allocatable :: rwork(:)

777

779

double precision, allocatable :: rwork(:)

double precision :: anorm

778

780

double precision :: anorm

double complex, allocatable :: wc_a(:,:) ! working copy of a

779

781

double complex, allocatable :: wc_a(:,:) ! working copy of a

integer :: info

780

782

integer :: info

integer, allocatable :: ipiv(:)

781

783

integer, allocatable :: ipiv(:)

782

784

double precision, external :: zlange

783

785

double precision, external :: zlange

784

786

785

787

if (present(status)) status=1

786

788

if (present(status)) status=1

787

789

anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm

788

790

anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm

789

791

allocate(wc_a(d,d))

790

792

allocate(wc_a(d,d))

!call zcopy(d*d,a,1,wc_a,1)

791

793

!call zcopy(d*d,a,1,wc_a,1)

wc_a(:,:)=a(:,:)

792

794

wc_a(:,:)=a(:,:)

793

795

allocate(ipiv(d))

794

796

allocate(ipiv(d))

call zgetrf(d,d,wc_a,d,ipiv,info)

795

797

call zgetrf(d,d,wc_a,d,ipiv,info)

if (info /= 0) then

796

798

if (info /= 0) then

if (present(status)) status=0

797

799

if (present(status)) status=0

deallocate(ipiv)

798

800

deallocate(ipiv)

deallocate(wc_a)

799

801

deallocate(wc_a)

return

800

802

return

end if

801

803

end if

802

804

allocate(work(2*d))

803

805

allocate(work(2*d))

allocate(rwork(2*d))

804

806

allocate(rwork(2*d))

call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)

805

807

call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)

if (info /= 0) then

806

808

if (info /= 0) then

if (present(status)) status=0

807

809

if (present(status)) status=0

end if

808

810

end if

deallocate(rwork)

809

811

deallocate(rwork)

deallocate(work)

810

812

deallocate(work)

deallocate(ipiv)

811

813

deallocate(ipiv)

deallocate(wc_a)

812

814

deallocate(wc_a)

end subroutine

813

815

end subroutine

814

816

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

815

817

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

816

818

!

! Valeurs propres/ Vecteurs propre

817

819

! Valeurs propres/ Vecteurs propre

!

818

820

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

819

821

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

820

822

subroutine fvn_s_matev(d,a,evala,eveca,status)

821

823

subroutine fvn_s_matev(d,a,evala,eveca,status)

!

822

824

!

! integer d (in) : matrice rank

823

825

! integer d (in) : matrice rank

! real a(d,d) (in) : The Matrix

824

826

! real a(d,d) (in) : The Matrix

! complex evala(d) (out) : eigenvalues

825

827

! complex evala(d) (out) : eigenvalues

! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector

826

828

! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector

! integer (out) : status =0 if something went wrong

827

829

! integer (out) : status =0 if something went wrong

!

828

830

!

! interfacing Lapack routine SGEEV

829

831

! interfacing Lapack routine SGEEV

implicit none

830

832

implicit none

integer, intent(in) :: d

831

833

integer, intent(in) :: d

real, intent(in) :: a(d,d)

832

834

real, intent(in) :: a(d,d)

complex, intent(out) :: evala(d)

833

835

complex, intent(out) :: evala(d)

complex, intent(out) :: eveca(d,d)

834

836

complex, intent(out) :: eveca(d,d)

integer, intent(out), optional :: status

835

837

integer, intent(out), optional :: status

836

838

real, allocatable :: wc_a(:,:) ! a working copy of a

837

839

real, allocatable :: wc_a(:,:) ! a working copy of a

integer :: info

838

840

integer :: info

integer :: lwork

839

841

integer :: lwork

real, allocatable :: wr(:),wi(:)

840

842

real, allocatable :: wr(:),wi(:)

real :: vl ! unused but necessary for the call

841

843

real :: vl ! unused but necessary for the call

real, allocatable :: vr(:,:)

842

844

real, allocatable :: vr(:,:)

real, allocatable :: work(:)

843

845

real, allocatable :: work(:)

real :: twork(1)

844

846

real :: twork(1)

integer i

845

847

integer i

integer j

846

848

integer j

847

849

if (present(status)) status=1

848

850

if (present(status)) status=1

849

851

! making a working copy of a

850

852

! making a working copy of a

allocate(wc_a(d,d))

851

853

allocate(wc_a(d,d))

!call scopy(d*d,a,1,wc_a,1)

852

854

!call scopy(d*d,a,1,wc_a,1)

wc_a(:,:)=a(:,:)

853

855

wc_a(:,:)=a(:,:)

854

856

allocate(wr(d))

855

857

allocate(wr(d))

allocate(wi(d))

856

858

allocate(wi(d))

allocate(vr(d,d))

857

859

allocate(vr(d,d))

! query optimal work size

858

860

! query optimal work size

call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)

859

861

call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)

lwork=int(twork(1))

860

862

lwork=int(twork(1))

allocate(work(lwork))

861

863

allocate(work(lwork))

call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)

862

864

call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)

863

865

if (info /= 0) then

864

866

if (info /= 0) then

if (present(status)) status=0

865

867

if (present(status)) status=0

deallocate(work)

866

868

deallocate(work)

deallocate(vr)

867

869

deallocate(vr)

deallocate(wi)

868

870

deallocate(wi)

deallocate(wr)

869

871

deallocate(wr)

deallocate(wc_a)

870

872

deallocate(wc_a)

return

871

873

return

end if

872

874

end if

873

875

! now fill in the results

874

876

! now fill in the results

i=1

875

877

i=1

do while(i<=d)

876

878

do while(i<=d)

evala(i)=cmplx(wr(i),wi(i))

877

879

evala(i)=cmplx(wr(i),wi(i))

if (wi(i) == 0.) then ! eigenvalue is real

878

880

if (wi(i) == 0.) then ! eigenvalue is real

eveca(:,i)=cmplx(vr(:,i),0.)

879

881

eveca(:,i)=cmplx(vr(:,i),0.)

else ! eigenvalue is complex

880

882

else ! eigenvalue is complex

evala(i+1)=cmplx(wr(i+1),wi(i+1))

881

883

evala(i+1)=cmplx(wr(i+1),wi(i+1))

eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))

882

884

eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))

eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))

883

885

eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))

i=i+1

884

886

i=i+1

end if

885

887

end if

i=i+1

886

888

i=i+1

enddo

887

889

enddo

deallocate(work)

888

890

deallocate(work)

deallocate(vr)

889

891

deallocate(vr)

deallocate(wi)

890

892

deallocate(wi)

deallocate(wr)

891

893

deallocate(wr)

deallocate(wc_a)

892

894

deallocate(wc_a)

893

895

end subroutine

894

896

end subroutine

895

897

subroutine fvn_d_matev(d,a,evala,eveca,status)

896

898

subroutine fvn_d_matev(d,a,evala,eveca,status)

!

897

899

!

! integer d (in) : matrice rank

898

900

! integer d (in) : matrice rank

! double precision a(d,d) (in) : The Matrix

899

901

! double precision a(d,d) (in) : The Matrix

! double complex evala(d) (out) : eigenvalues

900

902

! double complex evala(d) (out) : eigenvalues

! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector

901

903

! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector

! integer (out) : status =0 if something went wrong

902

904

! integer (out) : status =0 if something went wrong

!

903

905

!

! interfacing Lapack routine DGEEV

904

906

! interfacing Lapack routine DGEEV

implicit none

905

907

implicit none

integer, intent(in) :: d

906

908

integer, intent(in) :: d

double precision, intent(in) :: a(d,d)

907

909

double precision, intent(in) :: a(d,d)

double complex, intent(out) :: evala(d)

908

910

double complex, intent(out) :: evala(d)

double complex, intent(out) :: eveca(d,d)

909

911

double complex, intent(out) :: eveca(d,d)

integer, intent(out), optional :: status

910

912

integer, intent(out), optional :: status

911

913

double precision, allocatable :: wc_a(:,:) ! a working copy of a

912

914

double precision, allocatable :: wc_a(:,:) ! a working copy of a

integer :: info

913

915

integer :: info

integer :: lwork

914

916

integer :: lwork

double precision, allocatable :: wr(:),wi(:)

915

917

double precision, allocatable :: wr(:),wi(:)

double precision :: vl ! unused but necessary for the call

916

918

double precision :: vl ! unused but necessary for the call

double precision, allocatable :: vr(:,:)

917

919

double precision, allocatable :: vr(:,:)

double precision, allocatable :: work(:)

918

920

double precision, allocatable :: work(:)

double precision :: twork(1)

919

921

double precision :: twork(1)

integer i

920

922

integer i

integer j

921

923

integer j

922

924

if (present(status)) status=1

923

925

if (present(status)) status=1

924

926

! making a working copy of a

925

927

! making a working copy of a

allocate(wc_a(d,d))

926

928

allocate(wc_a(d,d))

!call dcopy(d*d,a,1,wc_a,1)

927

929

!call dcopy(d*d,a,1,wc_a,1)

wc_a(:,:)=a(:,:)

928

930

wc_a(:,:)=a(:,:)

929

931

allocate(wr(d))

930

932

allocate(wr(d))

allocate(wi(d))

931

933

allocate(wi(d))

allocate(vr(d,d))

932

934

allocate(vr(d,d))

! query optimal work size

933

935

! query optimal work size

call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)

934

936

call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)

lwork=int(twork(1))

935

937

lwork=int(twork(1))

allocate(work(lwork))

936

938

allocate(work(lwork))

call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)

937

939

call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)

938

940

if (info /= 0) then

939

941

if (info /= 0) then

if (present(status)) status=0

940

942

if (present(status)) status=0

deallocate(work)

941

943

deallocate(work)

deallocate(vr)

942

944

deallocate(vr)

deallocate(wi)

943

945

deallocate(wi)

deallocate(wr)

944

946

deallocate(wr)

deallocate(wc_a)

945

947

deallocate(wc_a)

return

946

948

return

end if

947

949

end if

948

950

! now fill in the results

949

951

! now fill in the results

i=1

950

952

i=1

do while(i<=d)

951

953

do while(i<=d)

evala(i)=dcmplx(wr(i),wi(i))

952

954

evala(i)=dcmplx(wr(i),wi(i))

if (wi(i) == 0.) then ! eigenvalue is real

953

955

if (wi(i) == 0.) then ! eigenvalue is real

eveca(:,i)=dcmplx(vr(:,i),0.)

954

956

eveca(:,i)=dcmplx(vr(:,i),0.)

else ! eigenvalue is complex

955

957

else ! eigenvalue is complex

evala(i+1)=dcmplx(wr(i+1),wi(i+1))

956

958

evala(i+1)=dcmplx(wr(i+1),wi(i+1))

eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))

957

959

eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))

eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))

958

960

eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))

i=i+1

959

961

i=i+1

end if

960

962

end if

i=i+1

961

963

i=i+1

enddo

962

964

enddo

963

965

deallocate(work)

964

966

deallocate(work)

deallocate(vr)

965

967

deallocate(vr)

deallocate(wi)

966

968

deallocate(wi)

deallocate(wr)

967

969

deallocate(wr)

deallocate(wc_a)

968

970

deallocate(wc_a)

969

971

end subroutine

970

972

end subroutine

971

973

subroutine fvn_c_matev(d,a,evala,eveca,status)

972

974

subroutine fvn_c_matev(d,a,evala,eveca,status)

!

973

975

!

! integer d (in) : matrice rank

974

976

! integer d (in) : matrice rank

! complex a(d,d) (in) : The Matrix

975

977

! complex a(d,d) (in) : The Matrix

! complex evala(d) (out) : eigenvalues

976

978

! complex evala(d) (out) : eigenvalues

! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector

977

979

! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector

! integer (out) : status =0 if something went wrong

978

980

! integer (out) : status =0 if something went wrong

!

979

981

!

! interfacing Lapack routine CGEEV

980

982

! interfacing Lapack routine CGEEV

implicit none

981

983

implicit none

integer, intent(in) :: d

982

984

integer, intent(in) :: d

complex, intent(in) :: a(d,d)

983

985

complex, intent(in) :: a(d,d)

complex, intent(out) :: evala(d)

984

986

complex, intent(out) :: evala(d)

complex, intent(out) :: eveca(d,d)

985

987

complex, intent(out) :: eveca(d,d)

integer, intent(out), optional :: status

986

988

integer, intent(out), optional :: status

987

989

complex, allocatable :: wc_a(:,:) ! a working copy of a

988

990

complex, allocatable :: wc_a(:,:) ! a working copy of a

integer :: info

989

991

integer :: info

integer :: lwork

990

992

integer :: lwork

complex, allocatable :: work(:)

991

993

complex, allocatable :: work(:)

complex :: twork(1)

992

994

complex :: twork(1)

real, allocatable :: rwork(:)

993

995

real, allocatable :: rwork(:)

complex :: vl ! unused but necessary for the call

994

996

complex :: vl ! unused but necessary for the call

995

997

if (present(status)) status=1

996

998

if (present(status)) status=1

997

999

! making a working copy of a

998

1000

! making a working copy of a

allocate(wc_a(d,d))

999

1001

allocate(wc_a(d,d))

!call ccopy(d*d,a,1,wc_a,1)

1000

1002

!call ccopy(d*d,a,1,wc_a,1)

wc_a(:,:)=a(:,:)

1001

1003

wc_a(:,:)=a(:,:)

1002

1004

1003

1005

! query optimal work size

1004

1006

! query optimal work size

call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)

1005

1007

call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)

lwork=int(twork(1))

1006

1008

lwork=int(twork(1))

allocate(work(lwork))

1007

1009

allocate(work(lwork))

allocate(rwork(2*d))

1008

1010

allocate(rwork(2*d))

call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)

1009

1011

call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)

1010

1012

if (info /= 0) then

1011

1013

if (info /= 0) then

if (present(status)) status=0

1012

1014

if (present(status)) status=0

end if

1013

1015

end if

deallocate(rwork)

1014

1016

deallocate(rwork)

deallocate(work)

1015

1017

deallocate(work)

deallocate(wc_a)

1016

1018

deallocate(wc_a)

1017

1019

end subroutine

1018

1020

end subroutine

1019

1021

subroutine fvn_z_matev(d,a,evala,eveca,status)

1020

1022

subroutine fvn_z_matev(d,a,evala,eveca,status)

!

1021

1023

!

! integer d (in) : matrice rank

1022

1024

! integer d (in) : matrice rank

! double complex a(d,d) (in) : The Matrix

1023

1025

! double complex a(d,d) (in) : The Matrix

! double complex evala(d) (out) : eigenvalues

1024

1026

! double complex evala(d) (out) : eigenvalues

! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector

1025

1027

! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector

! integer (out) : status =0 if something went wrong

1026

1028

! integer (out) : status =0 if something went wrong

!

1027

1029

!

! interfacing Lapack routine ZGEEV

1028

1030

! interfacing Lapack routine ZGEEV

implicit none

1029

1031

implicit none

integer, intent(in) :: d

1030

1032

integer, intent(in) :: d

double complex, intent(in) :: a(d,d)

1031

1033

double complex, intent(in) :: a(d,d)

double complex, intent(out) :: evala(d)

1032

1034

double complex, intent(out) :: evala(d)

double complex, intent(out) :: eveca(d,d)

1033

1035

double complex, intent(out) :: eveca(d,d)

integer, intent(out), optional :: status

1034

1036

integer, intent(out), optional :: status

1035

1037

double complex, allocatable :: wc_a(:,:) ! a working copy of a

1036

1038

double complex, allocatable :: wc_a(:,:) ! a working copy of a

integer :: info

1037

1039

integer :: info

integer :: lwork

1038

1040

integer :: lwork

double complex, allocatable :: work(:)

1039

1041

double complex, allocatable :: work(:)

double complex :: twork(1)

1040

1042

double complex :: twork(1)

double precision, allocatable :: rwork(:)

1041

1043

double precision, allocatable :: rwork(:)

double complex :: vl ! unused but necessary for the call

1042

1044

double complex :: vl ! unused but necessary for the call

1043

1045

if (present(status)) status=1

1044

1046

if (present(status)) status=1

1045

1047

! making a working copy of a

1046

1048

! making a working copy of a

allocate(wc_a(d,d))

1047

1049

allocate(wc_a(d,d))

!call zcopy(d*d,a,1,wc_a,1)

1048

1050

!call zcopy(d*d,a,1,wc_a,1)

wc_a(:,:)=a(:,:)

1049

1051

wc_a(:,:)=a(:,:)

1050

1052

! query optimal work size

1051

1053

! query optimal work size

call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)

1052

1054

call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)

lwork=int(twork(1))

1053

1055

lwork=int(twork(1))

allocate(work(lwork))

1054

1056

allocate(work(lwork))

allocate(rwork(2*d))

1055

1057

allocate(rwork(2*d))

call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)

1056

1058

call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)

1057

1059

if (info /= 0) then

1058

1060

if (info /= 0) then

if (present(status)) status=0

1059

1061

if (present(status)) status=0

end if

1060

1062

end if

deallocate(rwork)

1061

1063

deallocate(rwork)

deallocate(work)

1062

1064

deallocate(work)

deallocate(wc_a)

1063

1065

deallocate(wc_a)

1064

1066

end subroutine

1065

1067

end subroutine

1066

1068

1067

1069

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

1068

1070

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

1069

1071

!

! Quadratic interpolation of tabulated function of 1,2 or 3 variables

1070

1072

! Quadratic interpolation of tabulated function of 1,2 or 3 variables

!

1071

1073

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

1072

1074

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

1073

1075

subroutine fvn_s_find_interval(x,i,xdata,n)

1074

1076

subroutine fvn_s_find_interval(x,i,xdata,n)

implicit none

1075

1077

implicit none

! This routine find the indice i where xdata(i) <= x < xdata(i+1)

1076

1078

! This routine find the indice i where xdata(i) <= x < xdata(i+1)

! xdata(n) must contains a set of increasingly ordered values

1077

1079

! xdata(n) must contains a set of increasingly ordered values

! if x < xdata(1) i=0 is returned

1078

1080

! if x < xdata(1) i=0 is returned

! if x > xdata(n) i=n is returned

1079

1081

! if x > xdata(n) i=n is returned

! special case is where x=xdata(n) then n-1 is returned so

1080

1082

! special case is where x=xdata(n) then n-1 is returned so

! we will not exclude the upper limit

1081

1083

! we will not exclude the upper limit

! a simple dichotomy method is used

1082

1084

! a simple dichotomy method is used

1083

1085

real(kind=4), intent(in) :: x

1084

1086

real(kind=4), intent(in) :: x

real(kind=4), intent(in), dimension(n) :: xdata

1085

1087

real(kind=4), intent(in), dimension(n) :: xdata

integer(kind=4), intent(in) :: n

1086

1088

integer(kind=4), intent(in) :: n

integer(kind=4), intent(out) :: i

1087

1089

integer(kind=4), intent(out) :: i

1088

1090

integer(kind=4) :: imin,imax,imoyen

1089

1091

integer(kind=4) :: imin,imax,imoyen

1090

1092

! special case is where x=xdata(n) then n-1 is returned so

1091

1093

! special case is where x=xdata(n) then n-1 is returned so

! we will not exclude the upper limit

1092

1094

! we will not exclude the upper limit

if (x == xdata(n)) then

1093

1095

if (x == xdata(n)) then

i=n-1

1094

1096

i=n-1

return

1095

1097

return

end if

1096

1098

end if

1097

1099

! if x < xdata(1) i=0 is returned

1098

1100

! if x < xdata(1) i=0 is returned

if (x < xdata(1)) then

1099

1101

if (x < xdata(1)) then

i=0

1100

1102

i=0

return

1101

1103

return

end if

1102

1104

end if

1103

1105

! if x > xdata(n) i=n is returned

1104

1106

! if x > xdata(n) i=n is returned

if (x > xdata(n)) then

1105

1107

if (x > xdata(n)) then

i=n

1106

1108

i=n

return

1107

1109

return

end if

1108

1110

end if

1109

1111

! here xdata(1) <= x <= xdata(n)

1110

1112

! here xdata(1) <= x <= xdata(n)

imin=0

1111

1113

imin=0

imax=n+1

1112

1114

imax=n+1

1113

1115

do while((imax-imin) > 1)

1114

1116

do while((imax-imin) > 1)

imoyen=(imax+imin)/2

1115

1117

imoyen=(imax+imin)/2

if (x >= xdata(imoyen)) then

1116

1118

if (x >= xdata(imoyen)) then

imin=imoyen

1117

1119

imin=imoyen

else

1118

1120

else

imax=imoyen

1119

1121

imax=imoyen

end if

1120

1122

end if

end do

1121

1123

end do

1122

1124

i=imin

1123

1125

i=imin

1124

1126

end subroutine

1125

1127

end subroutine

1126

1128

1127

1129

subroutine fvn_d_find_interval(x,i,xdata,n)

1128

1130

subroutine fvn_d_find_interval(x,i,xdata,n)

implicit none

1129

1131

implicit none

! This routine find the indice i where xdata(i) <= x < xdata(i+1)

1130

1132

! This routine find the indice i where xdata(i) <= x < xdata(i+1)

! xdata(n) must contains a set of increasingly ordered values

1131

1133

! xdata(n) must contains a set of increasingly ordered values

! if x < xdata(1) i=0 is returned

1132

1134

! if x < xdata(1) i=0 is returned

! if x > xdata(n) i=n is returned

1133

1135

! if x > xdata(n) i=n is returned

! special case is where x=xdata(n) then n-1 is returned so

1134

1136

! special case is where x=xdata(n) then n-1 is returned so

! we will not exclude the upper limit

1135

1137

! we will not exclude the upper limit

! a simple dichotomy method is used

1136

1138

! a simple dichotomy method is used

1137

1139

real(kind=8), intent(in) :: x

1138

1140

real(kind=8), intent(in) :: x

real(kind=8), intent(in), dimension(n) :: xdata

1139

1141

real(kind=8), intent(in), dimension(n) :: xdata

integer(kind=4), intent(in) :: n

1140

1142

integer(kind=4), intent(in) :: n

integer(kind=4), intent(out) :: i

1141

1143

integer(kind=4), intent(out) :: i

1142

1144

integer(kind=4) :: imin,imax,imoyen

1143

1145

integer(kind=4) :: imin,imax,imoyen

1144

1146

! special case is where x=xdata(n) then n-1 is returned so

1145

1147

! special case is where x=xdata(n) then n-1 is returned so

! we will not exclude the upper limit

1146

1148

! we will not exclude the upper limit

if (x == xdata(n)) then

1147

1149

if (x == xdata(n)) then

i=n-1

1148

1150

i=n-1

return

1149

1151

return

end if

1150

1152

end if

1151

1153

! if x < xdata(1) i=0 is returned

1152

1154

! if x < xdata(1) i=0 is returned

if (x < xdata(1)) then

1153

1155

if (x < xdata(1)) then

i=0

1154

1156

i=0

return

1155

1157

return

end if

1156

1158

end if

1157

1159

! if x > xdata(n) i=n is returned

1158

1160

! if x > xdata(n) i=n is returned

if (x > xdata(n)) then

1159

1161

if (x > xdata(n)) then

i=n

1160

1162

i=n

return

1161

1163

return

end if

1162

1164

end if

1163

1165

! here xdata(1) <= x <= xdata(n)

1164

1166

! here xdata(1) <= x <= xdata(n)

imin=0

1165

1167

imin=0

imax=n+1

1166

1168

imax=n+1

1167

1169

do while((imax-imin) > 1)

1168

1170

do while((imax-imin) > 1)

imoyen=(imax+imin)/2

1169

1171

imoyen=(imax+imin)/2

if (x >= xdata(imoyen)) then

1170

1172

if (x >= xdata(imoyen)) then

imin=imoyen

1171

1173

imin=imoyen

else

1172

1174

else

imax=imoyen

1173

1175

imax=imoyen

end if

1174

1176

end if

end do

1175

1177

end do

1176

1178

i=imin

1177

1179

i=imin

1178

1180

end subroutine

1179

1181

end subroutine

1180

1182

1181

1183

function fvn_s_quad_interpol(x,n,xdata,ydata)

1182

1184

function fvn_s_quad_interpol(x,n,xdata,ydata)

implicit none

1183

1185

implicit none

! This function evaluate the value of a function defined by a set of points

1184

1186

! This function evaluate the value of a function defined by a set of points

! and values, using a quadratic interpolation

1185

1187

! and values, using a quadratic interpolation

! xdata must be increasingly ordered

1186

1188

! xdata must be increasingly ordered

! x must be within xdata(1) and xdata(n) to actually do interpolation

1187

1189

! x must be within xdata(1) and xdata(n) to actually do interpolation

! otherwise extrapolation is done

1188

1190

! otherwise extrapolation is done

integer(kind=4), intent(in) :: n

1189

1191

integer(kind=4), intent(in) :: n

real(kind=4), intent(in), dimension(n) :: xdata,ydata

1190

1192

real(kind=4), intent(in), dimension(n) :: xdata,ydata

real(kind=4), intent(in) :: x

1191

1193

real(kind=4), intent(in) :: x

real(kind=4) :: fvn_s_quad_interpol

1192

1194

real(kind=4) :: fvn_s_quad_interpol

1193

1195

integer(kind=4) :: iinf,base,i,j

1194

1196

integer(kind=4) :: iinf,base,i,j

real(kind=4) :: p

1195

1197

real(kind=4) :: p

1196

1198

call fvn_s_find_interval(x,iinf,xdata,n)

1197

1199

call fvn_s_find_interval(x,iinf,xdata,n)

1198

1200

! Settings for extrapolation

1199

1201

! Settings for extrapolation

if (iinf==0) then

1200

1202

if (iinf==0) then

! TODO -> Lower bound extrapolation warning

1201

1203

! TODO -> Lower bound extrapolation warning

iinf=1

1202

1204

iinf=1

end if

1203

1205

end if

1204

1206

if (iinf==n) then

1205

1207

if (iinf==n) then

! TODO -> Higher bound extrapolation warning

1206

1208

! TODO -> Higher bound extrapolation warning

iinf=n-1

1207

1209

iinf=n-1

end if

1208

1210

end if

1209

1211

! The three points we will use are iinf-1,iinf and iinf+1 with the

1210

1212

! The three points we will use are iinf-1,iinf and iinf+1 with the

! exception of the first interval, where iinf=1 we will use 1,2 and 3

1211

1213

! exception of the first interval, where iinf=1 we will use 1,2 and 3

if (iinf==1) then

1212

1214

if (iinf==1) then

base=0

1213

1215

base=0

else

1214

1216

else

base=iinf-2

1215

1217

base=iinf-2

end if

1216

1218

end if

1217

1219

! The three points we will use are :

1218

1220

! The three points we will use are :

! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)

1219

1221

! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)

1220

1222

! Straight forward Lagrange polynomial

1221

1223

! Straight forward Lagrange polynomial

fvn_s_quad_interpol=0.

1222

1224

fvn_s_quad_interpol=0.

do i=1,3

1223

1225

do i=1,3

! polynome i

1224

1226

! polynome i

p=ydata(base+i)

1225

1227

p=ydata(base+i)

do j=1,3

1226

1228

do j=1,3

if (j /= i) then

1227

1229

if (j /= i) then

p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))

1228

1230

p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))

end if

1229

1231

end if

end do

1230

1232

end do

fvn_s_quad_interpol=fvn_s_quad_interpol+p

1231

1233

fvn_s_quad_interpol=fvn_s_quad_interpol+p

end do

1232

1234

end do

1233

1235

end function

1234

1236

end function

1235

1237

1236

1238

function fvn_d_quad_interpol(x,n,xdata,ydata)

1237

1239

function fvn_d_quad_interpol(x,n,xdata,ydata)

implicit none

1238

1240

implicit none

! This function evaluate the value of a function defined by a set of points

1239

1241

! This function evaluate the value of a function defined by a set of points

! and values, using a quadratic interpolation

1240

1242

! and values, using a quadratic interpolation

! xdata must be increasingly ordered

1241

1243

! xdata must be increasingly ordered

! x must be within xdata(1) and xdata(n) to actually do interpolation

1242

1244

! x must be within xdata(1) and xdata(n) to actually do interpolation

! otherwise extrapolation is done

1243

1245

! otherwise extrapolation is done

integer(kind=4), intent(in) :: n

1244

1246

integer(kind=4), intent(in) :: n

real(kind=8), intent(in), dimension(n) :: xdata,ydata

1245

1247

real(kind=8), intent(in), dimension(n) :: xdata,ydata

real(kind=8), intent(in) :: x

1246

1248

real(kind=8), intent(in) :: x

real(kind=8) :: fvn_d_quad_interpol

1247

1249

real(kind=8) :: fvn_d_quad_interpol

1248

1250

integer(kind=4) :: iinf,base,i,j

1249

1251

integer(kind=4) :: iinf,base,i,j

real(kind=8) :: p

1250

1252

real(kind=8) :: p

1251

1253

call fvn_d_find_interval(x,iinf,xdata,n)

1252

1254

call fvn_d_find_interval(x,iinf,xdata,n)

1253

1255

! Settings for extrapolation

1254

1256

! Settings for extrapolation

if (iinf==0) then

1255

1257

if (iinf==0) then

! TODO -> Lower bound extrapolation warning

1256

1258

! TODO -> Lower bound extrapolation warning

iinf=1

1257

1259

iinf=1

end if

1258

1260

end if

1259

1261

if (iinf==n) then

1260

1262

if (iinf==n) then

! TODO Higher bound extrapolation warning

1261

1263

! TODO Higher bound extrapolation warning

iinf=n-1

1262

1264

iinf=n-1

end if

1263

1265

end if

1264

1266

! The three points we will use are iinf-1,iinf and iinf+1 with the

1265

1267

! The three points we will use are iinf-1,iinf and iinf+1 with the

! exception of the first interval, where iinf=1 we will use 1,2 and 3

1266

1268

! exception of the first interval, where iinf=1 we will use 1,2 and 3

if (iinf==1) then

1267

1269

if (iinf==1) then

base=0

1268

1270

base=0

else

1269

1271

else

base=iinf-2

1270

1272

base=iinf-2

end if

1271

1273

end if

1272

1274

! The three points we will use are :

1273

1275

! The three points we will use are :

! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)

1274

1276

! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)

1275

1277

! Straight forward Lagrange polynomial

1276

1278

! Straight forward Lagrange polynomial

fvn_d_quad_interpol=0.

1277

1279

fvn_d_quad_interpol=0.

do i=1,3

1278

1280

do i=1,3

! polynome i

1279

1281

! polynome i

p=ydata(base+i)

1280

1282

p=ydata(base+i)

do j=1,3

1281

1283

do j=1,3

if (j /= i) then

1282

1284

if (j /= i) then

p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))

1283

1285

p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))

end if

1284

1286

end if

end do

1285

1287

end do

fvn_d_quad_interpol=fvn_d_quad_interpol+p

1286

1288

fvn_d_quad_interpol=fvn_d_quad_interpol+p

end do

1287

1289

end do

1288

1290

end function

1289

1291

end function

1290

1292

1291

1293

function fvn_s_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)

1292

1294

function fvn_s_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)

implicit none

1293

1295

implicit none

! This function evaluate the value of a two variable function defined by a

1294

1296

! This function evaluate the value of a two variable function defined by a

! set of points and values, using a quadratic interpolation

1295

1297

! set of points and values, using a quadratic interpolation

! xdata and ydata must be increasingly ordered

1296

1298

! xdata and ydata must be increasingly ordered

! the couple (x,y) must be as x within xdata(1) and xdata(nx) and

1297

1299

! the couple (x,y) must be as x within xdata(1) and xdata(nx) and

! y within ydata(1) and ydata(ny) to actually do interpolation

1298

1300

! y within ydata(1) and ydata(ny) to actually do interpolation

! otherwise extrapolation is done

1299

1301

! otherwise extrapolation is done

integer(kind=4), intent(in) :: nx,ny

1300

1302

integer(kind=4), intent(in) :: nx,ny

real(kind=4), intent(in) :: x,y

1301

1303

real(kind=4), intent(in) :: x,y

real(kind=4), intent(in), dimension(nx) :: xdata

1302

1304

real(kind=4), intent(in), dimension(nx) :: xdata

real(kind=4), intent(in), dimension(ny) :: ydata

1303

1305

real(kind=4), intent(in), dimension(ny) :: ydata

real(kind=4), intent(in), dimension(nx,ny) :: zdata

1304

1306

real(kind=4), intent(in), dimension(nx,ny) :: zdata

real(kind=4) :: fvn_s_quad_2d_interpol

1305

1307

real(kind=4) :: fvn_s_quad_2d_interpol

1306

1308

integer(kind=4) :: ixinf,iyinf,basex,basey,i

1307

1309

integer(kind=4) :: ixinf,iyinf,basex,basey,i

real(kind=4),dimension(3) :: ztmp

1308

1310

real(kind=4),dimension(3) :: ztmp

!real(kind=4), external :: fvn_s_quad_interpol

1309

1311

!real(kind=4), external :: fvn_s_quad_interpol

1310

1312

call fvn_s_find_interval(x,ixinf,xdata,nx)

1311

1313

call fvn_s_find_interval(x,ixinf,xdata,nx)

call fvn_s_find_interval(y,iyinf,ydata,ny)

1312

1314

call fvn_s_find_interval(y,iyinf,ydata,ny)

1313

1315

! Settings for extrapolation

1314

1316

! Settings for extrapolation

if (ixinf==0) then

1315

1317

if (ixinf==0) then

! TODO -> Lower x bound extrapolation warning

1316

1318

! TODO -> Lower x bound extrapolation warning

ixinf=1

1317

1319

ixinf=1

end if

1318

1320

end if

1319

1321

if (ixinf==nx) then

1320

1322

if (ixinf==nx) then

! TODO -> Higher x bound extrapolation warning

1321

1323

! TODO -> Higher x bound extrapolation warning

ixinf=nx-1

1322

1324

ixinf=nx-1

end if

1323

1325

end if

1324

1326

if (iyinf==0) then

1325

1327

if (iyinf==0) then

! TODO -> Lower y bound extrapolation warning

1326

1328

! TODO -> Lower y bound extrapolation warning

iyinf=1

1327

1329

iyinf=1

end if

1328

1330

end if

1329

1331

if (iyinf==ny) then

1330

1332

if (iyinf==ny) then

! TODO -> Higher y bound extrapolation warning

1331

1333

! TODO -> Higher y bound extrapolation warning

iyinf=ny-1

1332

1334

iyinf=ny-1

end if

1333

1335

end if

1334

1336

! The three points we will use are iinf-1,iinf and iinf+1 with the

1335

1337

! The three points we will use are iinf-1,iinf and iinf+1 with the

! exception of the first interval, where iinf=1 we will use 1,2 and 3

1336

1338

! exception of the first interval, where iinf=1 we will use 1,2 and 3

if (ixinf==1) then

1337

1339

if (ixinf==1) then

basex=0

1338

1340

basex=0

else

1339

1341

else

basex=ixinf-2

1340

1342

basex=ixinf-2

end if

1341

1343

end if

1342

1344

if (iyinf==1) then

1343

1345

if (iyinf==1) then

basey=0

1344

1346

basey=0

else

1345

1347

else

basey=iyinf-2

1346

1348

basey=iyinf-2

end if

1347

1349

end if

1348

1350

! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)

1349

1351

! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)

! stored in ztmp(1:3)

1350

1352

! stored in ztmp(1:3)

do i=1,3

1351

1353

do i=1,3

ztmp(i)=fvn_s_quad_interpol(x,nx,xdata,zdata(:,basey+i))

1352

1354

ztmp(i)=fvn_s_quad_interpol(x,nx,xdata,zdata(:,basey+i))

end do

1353

1355

end do

1354

1356

! Then we make an interpolation for y using previous interpolations

1355

1357

! Then we make an interpolation for y using previous interpolations

fvn_s_quad_2d_interpol=fvn_s_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)

1356

1358

fvn_s_quad_2d_interpol=fvn_s_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)

end function

1357

1359

end function

1358

1360

1359

1361

function fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)

1360

1362

function fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)

implicit none

1361

1363

implicit none

! This function evaluate the value of a two variable function defined by a

1362

1364

! This function evaluate the value of a two variable function defined by a

! set of points and values, using a quadratic interpolation

1363

1365

! set of points and values, using a quadratic interpolation

! xdata and ydata must be increasingly ordered

1364

1366

! xdata and ydata must be increasingly ordered

! the couple (x,y) must be as x within xdata(1) and xdata(nx) and

1365

1367

! the couple (x,y) must be as x within xdata(1) and xdata(nx) and

! y within ydata(1) and ydata(ny) to actually do interpolation

1366

1368

! y within ydata(1) and ydata(ny) to actually do interpolation

! otherwise extrapolation is done

1367

1369

! otherwise extrapolation is done

integer(kind=4), intent(in) :: nx,ny

1368

1370

integer(kind=4), intent(in) :: nx,ny

real(kind=8), intent(in) :: x,y

1369

1371

real(kind=8), intent(in) :: x,y

real(kind=8), intent(in), dimension(nx) :: xdata

1370

1372

real(kind=8), intent(in), dimension(nx) :: xdata

real(kind=8), intent(in), dimension(ny) :: ydata

1371

1373

real(kind=8), intent(in), dimension(ny) :: ydata

real(kind=8), intent(in), dimension(nx,ny) :: zdata

1372

1374

real(kind=8), intent(in), dimension(nx,ny) :: zdata

real(kind=8) :: fvn_d_quad_2d_interpol

1373

1375

real(kind=8) :: fvn_d_quad_2d_interpol

1374

1376

integer(kind=4) :: ixinf,iyinf,basex,basey,i

1375

1377

integer(kind=4) :: ixinf,iyinf,basex,basey,i

real(kind=8),dimension(3) :: ztmp

1376

1378

real(kind=8),dimension(3) :: ztmp

!real(kind=8), external :: fvn_d_quad_interpol

1377

1379

!real(kind=8), external :: fvn_d_quad_interpol

1378

1380

call fvn_d_find_interval(x,ixinf,xdata,nx)

1379

1381

call fvn_d_find_interval(x,ixinf,xdata,nx)

call fvn_d_find_interval(y,iyinf,ydata,ny)

1380

1382

call fvn_d_find_interval(y,iyinf,ydata,ny)

1381

1383

! Settings for extrapolation

1382

1384

! Settings for extrapolation

if (ixinf==0) then

1383

1385

if (ixinf==0) then

! TODO -> Lower x bound extrapolation warning

1384

1386

! TODO -> Lower x bound extrapolation warning

ixinf=1

1385

1387

ixinf=1

end if

1386

1388

end if

1387

1389

if (ixinf==nx) then

1388

1390

if (ixinf==nx) then

! TODO -> Higher x bound extrapolation warning

1389

1391

! TODO -> Higher x bound extrapolation warning

ixinf=nx-1

1390

1392

ixinf=nx-1

end if

1391

1393

end if

1392

1394

if (iyinf==0) then

1393

1395

if (iyinf==0) then

! TODO -> Lower y bound extrapolation warning

1394

1396

! TODO -> Lower y bound extrapolation warning

iyinf=1

1395

1397

iyinf=1

end if

1396

1398

end if

1397

1399

if (iyinf==ny) then

1398

1400

if (iyinf==ny) then

! TODO -> Higher y bound extrapolation warning

1399

1401

! TODO -> Higher y bound extrapolation warning

iyinf=ny-1

1400

1402

iyinf=ny-1

end if

1401

1403

end if

1402

1404

! The three points we will use are iinf-1,iinf and iinf+1 with the

1403

1405

! The three points we will use are iinf-1,iinf and iinf+1 with the

! exception of the first interval, where iinf=1 we will use 1,2 and 3

1404

1406

! exception of the first interval, where iinf=1 we will use 1,2 and 3

if (ixinf==1) then

1405

1407

if (ixinf==1) then

basex=0

1406

1408

basex=0

else

1407

1409

else

basex=ixinf-2

1408

1410

basex=ixinf-2

end if

1409

1411

end if

1410

1412

if (iyinf==1) then

1411

1413

if (iyinf==1) then

basey=0

1412

1414

basey=0

else

1413

1415

else

basey=iyinf-2

1414

1416

basey=iyinf-2

end if

1415

1417

end if

1416

1418

! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)

1417

1419

! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)

! stored in ztmp(1:3)

1418

1420

! stored in ztmp(1:3)

do i=1,3

1419

1421

do i=1,3

ztmp(i)=fvn_d_quad_interpol(x,nx,xdata,zdata(:,basey+i))

1420

1422

ztmp(i)=fvn_d_quad_interpol(x,nx,xdata,zdata(:,basey+i))

end do

1421

1423

end do

1422

1424

! Then we make an interpolation for y using previous interpolations

1423

1425

! Then we make an interpolation for y using previous interpolations

fvn_d_quad_2d_interpol=fvn_d_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)

1424

1426

fvn_d_quad_2d_interpol=fvn_d_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)

end function

1425

1427

end function

1426

1428

1427

1429

function fvn_s_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)

1428

1430

function fvn_s_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)

implicit none

1429

1431

implicit none

! This function evaluate the value of a 3 variables function defined by a

1430

1432

! This function evaluate the value of a 3 variables function defined by a

! set of points and values, using a quadratic interpolation

1431

1433

! set of points and values, using a quadratic interpolation

! xdata, ydata and zdata must be increasingly ordered

1432

1434

! xdata, ydata and zdata must be increasingly ordered

! The triplet (x,y,z) must be within xdata,ydata and zdata to actually

1433

1435

! The triplet (x,y,z) must be within xdata,ydata and zdata to actually

! perform an interpolation, otherwise extrapolation is done

1434

1436

! perform an interpolation, otherwise extrapolation is done

integer(kind=4), intent(in) :: nx,ny,nz

1435

1437

integer(kind=4), intent(in) :: nx,ny,nz

real(kind=4), intent(in) :: x,y,z

1436

1438

real(kind=4), intent(in) :: x,y,z

real(kind=4), intent(in), dimension(nx) :: xdata

1437

1439

real(kind=4), intent(in), dimension(nx) :: xdata

real(kind=4), intent(in), dimension(ny) :: ydata

1438

1440

real(kind=4), intent(in), dimension(ny) :: ydata

real(kind=4), intent(in), dimension(nz) :: zdata

1439

1441

real(kind=4), intent(in), dimension(nz) :: zdata

real(kind=4), intent(in), dimension(nx,ny,nz) :: tdata

1440

1442

real(kind=4), intent(in), dimension(nx,ny,nz) :: tdata

real(kind=4) :: fvn_s_quad_3d_interpol

1441

1443

real(kind=4) :: fvn_s_quad_3d_interpol

1442

1444

integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j

1443

1445

integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j

!real(kind=4), external :: fvn_s_quad_interpol,fvn_s_quad_2d_interpol

1444

1446

!real(kind=4), external :: fvn_s_quad_interpol,fvn_s_quad_2d_interpol

real(kind=4),dimension(3,3) :: ttmp

1445

1447

real(kind=4),dimension(3,3) :: ttmp

1446

1448

call fvn_s_find_interval(x,ixinf,xdata,nx)

1447

1449

call fvn_s_find_interval(x,ixinf,xdata,nx)

call fvn_s_find_interval(y,iyinf,ydata,ny)

1448

1450

call fvn_s_find_interval(y,iyinf,ydata,ny)

call fvn_s_find_interval(z,izinf,zdata,nz)

1449

1451

call fvn_s_find_interval(z,izinf,zdata,nz)

1450

1452

! Settings for extrapolation

1451

1453

! Settings for extrapolation

if (ixinf==0) then

1452

1454

if (ixinf==0) then

! TODO -> Lower x bound extrapolation warning

1453

1455

! TODO -> Lower x bound extrapolation warning

ixinf=1

1454

1456

ixinf=1

end if

1455

1457

end if

1456

1458

if (ixinf==nx) then

1457

1459

if (ixinf==nx) then

! TODO -> Higher x bound extrapolation warning

1458

1460

! TODO -> Higher x bound extrapolation warning

ixinf=nx-1

1459

1461

ixinf=nx-1

end if

1460

1462

end if

1461

1463

if (iyinf==0) then

1462

1464

if (iyinf==0) then

! TODO -> Lower y bound extrapolation warning

1463

1465

! TODO -> Lower y bound extrapolation warning

iyinf=1

1464

1466

iyinf=1

end if

1465

1467

end if

1466

1468

if (iyinf==ny) then

1467

1469

if (iyinf==ny) then

! TODO -> Higher y bound extrapolation warning

1468

1470

! TODO -> Higher y bound extrapolation warning

iyinf=ny-1

1469

1471

iyinf=ny-1

end if

1470

1472

end if

1471

1473

if (izinf==0) then

1472

1474

if (izinf==0) then

! TODO -> Lower z bound extrapolation warning

1473

1475

! TODO -> Lower z bound extrapolation warning

izinf=1

1474

1476

izinf=1

end if

1475

1477

end if

1476

1478

if (izinf==nz) then

1477

1479

if (izinf==nz) then

! TODO -> Higher z bound extrapolation warning

1478

1480

! TODO -> Higher z bound extrapolation warning

izinf=nz-1

1479

1481

izinf=nz-1

end if

1480

1482

end if

1481

1483

! The three points we will use are iinf-1,iinf and iinf+1 with the

1482

1484

! The three points we will use are iinf-1,iinf and iinf+1 with the

! exception of the first interval, where iinf=1 we will use 1,2 and 3

1483

1485

! exception of the first interval, where iinf=1 we will use 1,2 and 3

if (ixinf==1) then

1484

1486

if (ixinf==1) then

basex=0

1485

1487

basex=0

else

1486

1488

else

basex=ixinf-2

1487

1489

basex=ixinf-2

end if

1488

1490

end if

1489

1491

if (iyinf==1) then

1490

1492

if (iyinf==1) then

basey=0

1491

1493

basey=0

else

1492

1494

else

basey=iyinf-2

1493

1495

basey=iyinf-2

end if

1494

1496

end if

1495

1497

if (izinf==1) then

1496

1498

if (izinf==1) then

basez=0

1497

1499

basez=0

else

1498

1500

else

basez=izinf-2

1499

1501

basez=izinf-2

end if

1500

1502

end if

1501

1503

! We first make 9 one dimensional interpolation on variable x.

1502

1504

! We first make 9 one dimensional interpolation on variable x.

! results are stored in ttmp

1503

1505

! results are stored in ttmp

do i=1,3

1504

1506

do i=1,3

do j=1,3

1505

1507

do j=1,3

ttmp(i,j)=fvn_s_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))

1506

1508

ttmp(i,j)=fvn_s_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))

end do

1507

1509

end do

1508

1510

end do

1509

1511

! We then make a 2 dimensionnal interpolation on variables y and z

1510

1512

! We then make a 2 dimensionnal interpolation on variables y and z

fvn_s_quad_3d_interpol=fvn_s_quad_2d_interpol(y,z, &

1511

1513

fvn_s_quad_3d_interpol=fvn_s_quad_2d_interpol(y,z, &

3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)

1512

1514

3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)

end function

1513

1515

end function

1514

1516

1515

1517

function fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)

1516

1518

function fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)

implicit none

1517

1519

implicit none

! This function evaluate the value of a 3 variables function defined by a

1518

1520

! This function evaluate the value of a 3 variables function defined by a

! set of points and values, using a quadratic interpolation

1519

1521

! set of points and values, using a quadratic interpolation

! xdata, ydata and zdata must be increasingly ordered

1520

1522

! xdata, ydata and zdata must be increasingly ordered

! The triplet (x,y,z) must be within xdata,ydata and zdata to actually

1521

1523

! The triplet (x,y,z) must be within xdata,ydata and zdata to actually

! perform an interpolation, otherwise extrapolation is done

1522

1524

! perform an interpolation, otherwise extrapolation is done

integer(kind=4), intent(in) :: nx,ny,nz

1523

1525

integer(kind=4), intent(in) :: nx,ny,nz

real(kind=8), intent(in) :: x,y,z

1524

1526

real(kind=8), intent(in) :: x,y,z

real(kind=8), intent(in), dimension(nx) :: xdata

1525

1527

real(kind=8), intent(in), dimension(nx) :: xdata

real(kind=8), intent(in), dimension(ny) :: ydata

1526

1528

real(kind=8), intent(in), dimension(ny) :: ydata

real(kind=8), intent(in), dimension(nz) :: zdata

1527

1529

real(kind=8), intent(in), dimension(nz) :: zdata

real(kind=8), intent(in), dimension(nx,ny,nz) :: tdata

1528

1530

real(kind=8), intent(in), dimension(nx,ny,nz) :: tdata

real(kind=8) :: fvn_d_quad_3d_interpol

1529

1531

real(kind=8) :: fvn_d_quad_3d_interpol

1530

1532

integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j

1531

1533

integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j

!real(kind=8), external :: fvn_d_quad_interpol,fvn_d_quad_2d_interpol

1532

1534

!real(kind=8), external :: fvn_d_quad_interpol,fvn_d_quad_2d_interpol

real(kind=8),dimension(3,3) :: ttmp

1533

1535

real(kind=8),dimension(3,3) :: ttmp

1534

1536

call fvn_d_find_interval(x,ixinf,xdata,nx)

1535

1537

call fvn_d_find_interval(x,ixinf,xdata,nx)

call fvn_d_find_interval(y,iyinf,ydata,ny)

1536

1538

call fvn_d_find_interval(y,iyinf,ydata,ny)

call fvn_d_find_interval(z,izinf,zdata,nz)

1537

1539

call fvn_d_find_interval(z,izinf,zdata,nz)

1538

1540

! Settings for extrapolation

1539

1541

! Settings for extrapolation

if (ixinf==0) then

1540

1542

if (ixinf==0) then

! TODO -> Lower x bound extrapolation warning

1541

1543

! TODO -> Lower x bound extrapolation warning

ixinf=1

1542

1544

ixinf=1

end if

1543

1545

end if

1544

1546

if (ixinf==nx) then

1545

1547

if (ixinf==nx) then

! TODO -> Higher x bound extrapolation warning

1546

1548

! TODO -> Higher x bound extrapolation warning

ixinf=nx-1

1547

1549

ixinf=nx-1

end if

1548

1550

end if

1549

1551

if (iyinf==0) then

1550

1552

if (iyinf==0) then

! TODO -> Lower y bound extrapolation warning

1551

1553

! TODO -> Lower y bound extrapolation warning

iyinf=1

1552

1554

iyinf=1

end if

1553

1555

end if

1554

1556

if (iyinf==ny) then

1555

1557

if (iyinf==ny) then

! TODO -> Higher y bound extrapolation warning

1556

1558

! TODO -> Higher y bound extrapolation warning

iyinf=ny-1

1557

1559

iyinf=ny-1

end if

1558

1560

end if

1559

1561

if (izinf==0) then

1560

1562

if (izinf==0) then

! TODO -> Lower z bound extrapolation warning

1561

1563

! TODO -> Lower z bound extrapolation warning

izinf=1

1562

1564

izinf=1

end if

1563

1565

end if

1564

1566

if (izinf==nz) then

1565

1567

if (izinf==nz) then

! TODO -> Higher z bound extrapolation warning

1566

1568

! TODO -> Higher z bound extrapolation warning

izinf=nz-1

1567

1569

izinf=nz-1

end if

1568

1570

end if

1569

1571

! The three points we will use are iinf-1,iinf and iinf+1 with the

1570

1572

! The three points we will use are iinf-1,iinf and iinf+1 with the

! exception of the first interval, where iinf=1 we will use 1,2 and 3

1571

1573

! exception of the first interval, where iinf=1 we will use 1,2 and 3

if (ixinf==1) then

1572

1574

if (ixinf==1) then

basex=0

1573

1575

basex=0

else

1574

1576

else

basex=ixinf-2

1575

1577

basex=ixinf-2

end if

1576

1578

end if

1577

1579

if (iyinf==1) then

1578

1580

if (iyinf==1) then

basey=0

1579

1581

basey=0

else

1580

1582

else

basey=iyinf-2

1581

1583

basey=iyinf-2

end if

1582

1584

end if

1583

1585

if (izinf==1) then

1584

1586

if (izinf==1) then

basez=0

1585

1587

basez=0

else

1586

1588

else

basez=izinf-2

1587

1589

basez=izinf-2

end if

1588

1590

end if

1589

1591

! We first make 9 one dimensional interpolation on variable x.

1590

1592

! We first make 9 one dimensional interpolation on variable x.

! results are stored in ttmp

1591

1593

! results are stored in ttmp

do i=1,3

1592

1594

do i=1,3

do j=1,3

1593

1595

do j=1,3

ttmp(i,j)=fvn_d_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))

1594

1596

ttmp(i,j)=fvn_d_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))

end do

1595

1597

end do

1596

1598

end do

1597

1599

! We then make a 2 dimensionnal interpolation on variables y and z

1598

1600

! We then make a 2 dimensionnal interpolation on variables y and z

fvn_d_quad_3d_interpol=fvn_d_quad_2d_interpol(y,z, &

1599

1601

fvn_d_quad_3d_interpol=fvn_d_quad_2d_interpol(y,z, &

3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)

1600

1602

3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)

end function

1601

1603

end function

1602

1604

1603

1605

1604

1606

1605

1607

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

1606

1608

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

1607

1609

!

! Akima spline interpolation and spline evaluation

1608

1610

! Akima spline interpolation and spline evaluation

!

1609

1611

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

1610

1612

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

1611

1613

! Single precision

1612

1614

! Single precision

subroutine fvn_s_akima(n,x,y,br,co)

1613

1615

subroutine fvn_s_akima(n,x,y,br,co)

implicit none

1614

1616

implicit none

integer, intent(in) :: n

1615

1617

integer, intent(in) :: n

real, intent(in) :: x(n)

1616

1618

real, intent(in) :: x(n)

real, intent(in) :: y(n)

1617

1619

real, intent(in) :: y(n)

real, intent(out) :: br(n)

1618

1620

real, intent(out) :: br(n)

real, intent(out) :: co(4,n)

1619

1621

real, intent(out) :: co(4,n)

1620

1622

real, allocatable :: var(:),z(:)

1621

1623

real, allocatable :: var(:),z(:)

real :: wi_1,wi

1622

1624

real :: wi_1,wi

integer :: i

1623

1625

integer :: i

real :: dx,a,b

1624

1626

real :: dx,a,b

1625

1627

! br is just a copy of x

1626

1628

! br is just a copy of x

br(:)=x(:)

1627

1629

br(:)=x(:)

1628

1630

allocate(var(n+3))

1629

1631

allocate(var(n+3))

allocate(z(n))

1630

1632

allocate(z(n))

! evaluate the variations

1631

1633

! evaluate the variations

do i=1, n-1

1632

1634

do i=1, n-1

var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))

1633

1635

var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))

end do

1634

1636

end do

var(n+2)=2.e0*var(n+1)-var(n)

1635

1637

var(n+2)=2.e0*var(n+1)-var(n)

var(n+3)=2.e0*var(n+2)-var(n+1)

1636

1638

var(n+3)=2.e0*var(n+2)-var(n+1)

var(2)=2.e0*var(3)-var(4)

1637

1639

var(2)=2.e0*var(3)-var(4)

var(1)=2.e0*var(2)-var(3)

1638

1640

var(1)=2.e0*var(2)-var(3)

1639

1641

do i = 1, n

1640

1642

do i = 1, n

wi_1=abs(var(i+3)-var(i+2))

1641

1643

wi_1=abs(var(i+3)-var(i+2))

wi=abs(var(i+1)-var(i))

1642

1644

wi=abs(var(i+1)-var(i))

if ((wi_1+wi).eq.0.e0) then

1643

1645

if ((wi_1+wi).eq.0.e0) then

z(i)=(var(i+2)+var(i+1))/2.e0

1644

1646

z(i)=(var(i+2)+var(i+1))/2.e0

else

1645

1647

else

z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)

1646

1648

z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)

end if

1647

1649

end if

end do

1648

1650

end do

1649

1651

do i=1, n-1

1650

1652

do i=1, n-1

dx=x(i+1)-x(i)

1651

1653

dx=x(i+1)-x(i)

a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd

1652

1654

a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd

b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd

1653

1655

b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd

co(1,i)=y(i)

1654

1656

co(1,i)=y(i)

co(2,i)=z(i)

1655

1657

co(2,i)=z(i)

!co(3,i)=-(a-3.*b)/dx**2 ! méthode wd

1656

1658

!co(3,i)=-(a-3.*b)/dx**2 ! méthode wd

!co(4,i)=(a-2.*b)/dx**3 ! méthode wd

1657

1659

!co(4,i)=(a-2.*b)/dx**3 ! méthode wd

co(3,i)=(3.e0*var(i+2)-2.e0*z(i)-z(i+1))/dx ! méthode JP Moreau

1658

1660

co(3,i)=(3.e0*var(i+2)-2.e0*z(i)-z(i+1))/dx ! méthode JP Moreau

co(4,i)=(z(i)+z(i+1)-2.e0*var(i+2))/dx**2 !

1659

1661

co(4,i)=(z(i)+z(i+1)-2.e0*var(i+2))/dx**2 !

! les coefficients donnés par imsl sont co(3,i)*2 et co(4,i)*6

1660

1662

! les coefficients donnés par imsl sont co(3,i)*2 et co(4,i)*6

! etrangement la fonction csval corrige et donne la bonne valeur ...

1661

1663

! etrangement la fonction csval corrige et donne la bonne valeur ...

end do

1662

1664

end do

co(1,n)=y(n)

1663

1665

co(1,n)=y(n)

co(2,n)=z(n)

1664

1666

co(2,n)=z(n)

co(3,n)=0.e0

1665

1667

co(3,n)=0.e0

co(4,n)=0.e0

1666

1668

co(4,n)=0.e0

1667

1669

deallocate(z)

1668

1670

deallocate(z)

deallocate(var)

1669

1671

deallocate(var)

1670

1672

end subroutine

1671

1673

end subroutine

1672

1674

! Double precision

1673

1675

! Double precision

subroutine fvn_d_akima(n,x,y,br,co)

1674

1676

subroutine fvn_d_akima(n,x,y,br,co)

1675

1677

implicit none

1676

1678

implicit none

integer, intent(in) :: n

1677

1679

integer, intent(in) :: n

double precision, intent(in) :: x(n)

1678

1680

double precision, intent(in) :: x(n)

double precision, intent(in) :: y(n)

1679

1681

double precision, intent(in) :: y(n)

double precision, intent(out) :: br(n)

1680

1682

double precision, intent(out) :: br(n)

double precision, intent(out) :: co(4,n)

1681

1683

double precision, intent(out) :: co(4,n)

1682

1684

double precision, allocatable :: var(:),z(:)

1683

1685

double precision, allocatable :: var(:),z(:)

double precision :: wi_1,wi

1684

1686

double precision :: wi_1,wi

integer :: i

1685

1687

integer :: i

double precision :: dx,a,b

1686

1688

double precision :: dx,a,b

1687

1689

! br is just a copy of x

1688

1690

! br is just a copy of x

br(:)=x(:)

1689

1691

br(:)=x(:)

1690

1692

allocate(var(n+3))

1691

1693

allocate(var(n+3))

allocate(z(n))

1692

1694

allocate(z(n))

! evaluate the variations

1693

1695

! evaluate the variations

do i=1, n-1

1694

1696

do i=1, n-1

var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))

1695

1697

var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))

end do

1696

1698

end do

var(n+2)=2.d0*var(n+1)-var(n)

1697

1699

var(n+2)=2.d0*var(n+1)-var(n)

var(n+3)=2.d0*var(n+2)-var(n+1)

1698

1700

var(n+3)=2.d0*var(n+2)-var(n+1)

var(2)=2.d0*var(3)-var(4)

1699

1701

var(2)=2.d0*var(3)-var(4)

var(1)=2.d0*var(2)-var(3)

1700

1702

var(1)=2.d0*var(2)-var(3)

1701

1703

do i = 1, n

1702

1704

do i = 1, n

wi_1=dabs(var(i+3)-var(i+2))

1703

1705

wi_1=dabs(var(i+3)-var(i+2))

wi=dabs(var(i+1)-var(i))

1704

1706

wi=dabs(var(i+1)-var(i))

if ((wi_1+wi).eq.0.d0) then

1705

1707

if ((wi_1+wi).eq.0.d0) then

z(i)=(var(i+2)+var(i+1))/2.d0

1706

1708

z(i)=(var(i+2)+var(i+1))/2.d0

else

1707

1709

else

z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)

1708

1710

z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)

end if

1709

1711

end if

end do

1710

1712

end do

1711

1713

do i=1, n-1

1712

1714

do i=1, n-1

dx=x(i+1)-x(i)

1713

1715

dx=x(i+1)-x(i)

a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd

1714

1716

a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd

b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd

1715

1717

b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd

co(1,i)=y(i)

1716

1718

co(1,i)=y(i)

co(2,i)=z(i)

1717

1719

co(2,i)=z(i)

!co(3,i)=-(a-3.*b)/dx**2 ! méthode wd

1718

1720

!co(3,i)=-(a-3.*b)/dx**2 ! méthode wd

!co(4,i)=(a-2.*b)/dx**3 ! méthode wd

1719

1721

!co(4,i)=(a-2.*b)/dx**3 ! méthode wd

co(3,i)=(3.d0*var(i+2)-2.d0*z(i)-z(i+1))/dx ! méthode JP Moreau

1720

1722

co(3,i)=(3.d0*var(i+2)-2.d0*z(i)-z(i+1))/dx ! méthode JP Moreau

co(4,i)=(z(i)+z(i+1)-2.d0*var(i+2))/dx**2 !

1721

1723

co(4,i)=(z(i)+z(i+1)-2.d0*var(i+2))/dx**2 !

! les coefficients donnés par imsl sont co(3,i)*2 et co(4,i)*6

1722

1724

! les coefficients donnés par imsl sont co(3,i)*2 et co(4,i)*6

! etrangement la fonction csval corrige et donne la bonne valeur ...

1723

1725

! etrangement la fonction csval corrige et donne la bonne valeur ...

end do

1724

1726

end do

co(1,n)=y(n)

1725

1727

co(1,n)=y(n)

co(2,n)=z(n)

1726

1728

co(2,n)=z(n)

co(3,n)=0.d0

1727

1729

co(3,n)=0.d0

co(4,n)=0.d0

1728

1730

co(4,n)=0.d0

1729

1731

deallocate(z)

1730

1732

deallocate(z)

deallocate(var)

1731

1733

deallocate(var)

1732

1734

end subroutine

1733

1735

end subroutine

1734

1736

!

1735

1737

!

! Single precision spline evaluation

1736

1738

! Single precision spline evaluation

!

1737

1739

!

function fvn_s_spline_eval(x,n,br,co)

1738

1740

function fvn_s_spline_eval(x,n,br,co)

implicit none

1739

1741

implicit none

real, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated

1740

1742

real, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated

integer, intent(in) :: n ! number of intervals

1741

1743

integer, intent(in) :: n ! number of intervals

real, intent(in) :: br(n+1) ! breakpoints

1742

1744

real, intent(in) :: br(n+1) ! breakpoints

real, intent(in) :: co(4,n+1) ! spline coeeficients

1743

1745

real, intent(in) :: co(4,n+1) ! spline coeeficients

real :: fvn_s_spline_eval

1744

1746

real :: fvn_s_spline_eval

1745

1747

integer :: i

1746

1748

integer :: i

real :: dx

1747

1749

real :: dx

1748

1750

if (x<=br(1)) then

1749

1751

if (x<=br(1)) then

i=1

1750

1752

i=1

else if (x>=br(n+1)) then

1751

1753

else if (x>=br(n+1)) then

i=n

1752

1754

i=n

else

1753

1755

else

i=1

1754

1756

i=1

do while(x>=br(i))

1755

1757

do while(x>=br(i))

i=i+1

1756

1758

i=i+1

end do

1757

1759

end do

i=i-1

1758

1760

i=i-1

end if

1759

1761

end if

dx=x-br(i)

1760

1762

dx=x-br(i)

fvn_s_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3

1761

1763

fvn_s_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3

1762

1764

end function

1763

1765

end function

1764

1766

! Double precision spline evaluation

1765

1767

! Double precision spline evaluation

function fvn_d_spline_eval(x,n,br,co)

1766

1768

function fvn_d_spline_eval(x,n,br,co)

implicit none

1767

1769

implicit none

double precision, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated

1768

1770

double precision, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated

integer, intent(in) :: n ! number of intervals

1769

1771

integer, intent(in) :: n ! number of intervals

double precision, intent(in) :: br(n+1) ! breakpoints

1770

1772

double precision, intent(in) :: br(n+1) ! breakpoints

double precision, intent(in) :: co(4,n+1) ! spline coeeficients

1771

1773

double precision, intent(in) :: co(4,n+1) ! spline coeeficients

double precision :: fvn_d_spline_eval

1772

1774

double precision :: fvn_d_spline_eval

1773

1775

integer :: i

1774

1776

integer :: i

double precision :: dx

1775

1777

double precision :: dx

1776

1778

1777

1779

if (x<=br(1)) then

1778

1780

if (x<=br(1)) then

i=1

1779

1781

i=1

else if (x>=br(n+1)) then

1780

1782

else if (x>=br(n+1)) then

i=n

1781

1783

i=n

else

1782

1784

else

i=1

1783

1785

i=1

do while(x>=br(i))

1784

1786

do while(x>=br(i))

i=i+1

1785

1787

i=i+1

end do

1786

1788

end do

i=i-1

1787

1789

i=i-1

end if

1788

1790

end if

1789

1791

dx=x-br(i)

1790

1792

dx=x-br(i)

fvn_d_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3

1791

1793

fvn_d_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3

1792

1794

end function

1793

1795

end function

1794

1796

1795

1797

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

1796

1798

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

1797

1799

!

! Least square problem

1798

1800

! Least square problem

!

1799

1801

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

1800

1802

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

1801

1803

!

1802

1804

!

1803

1805

1804

1806

1805

1807

1806

1808

subroutine fvn_d_lspoly(np,x,y,deg,coeff,status)

1807

1809

subroutine fvn_d_lspoly(np,x,y,deg,coeff,status)

!

1808

1810

!

! Least square polynomial fitting

1809

1811

! Least square polynomial fitting

!

1810

1812

!

! Find the coefficients of the least square polynomial of a given degree

1811

1813

! Find the coefficients of the least square polynomial of a given degree

! for a set of coordinates.

1812

1814

! for a set of coordinates.

!

1813

1815

!

! The degree must be lower than the number of points

1814

1816

! The degree must be lower than the number of points

!

1815

1817

!

! np (in) : number of points

1816

1818

! np (in) : number of points

! x(np) (in) : x data

1817

1819

! x(np) (in) : x data

! y(np) (in) : y data

1818

1820

! y(np) (in) : y data

! deg (in) : polynomial's degree

1819

1821

! deg (in) : polynomial's degree

! coeff(deg+1) (out) : polynomial coefficients

1820

1822

! coeff(deg+1) (out) : polynomial coefficients

! status (out) : =0 if a problem occurs

1821

1823

! status (out) : =0 if a problem occurs

!

1822

1824

!

implicit none

1823

1825

implicit none

1824

1826

integer, intent(in) :: np,deg

1825

1827

integer, intent(in) :: np,deg

real(kind=8), intent(in), dimension(np) :: x,y

1826

1828

real(kind=8), intent(in), dimension(np) :: x,y

real(kind=8), intent(out), dimension(deg+1) :: coeff

1827

1829

real(kind=8), intent(out), dimension(deg+1) :: coeff

integer, intent(out), optional :: status

1828

1830

integer, intent(out), optional :: status

1829

1831

real(kind=8), allocatable, dimension(:,:) :: mat,bmat

1830

1832

real(kind=8), allocatable, dimension(:,:) :: mat,bmat

real(kind=8),dimension(:),allocatable :: work

1831

1833

real(kind=8),dimension(:),allocatable :: work

real(kind=8),dimension(1) :: twork

1832

1834

real(kind=8),dimension(1) :: twork

integer :: lwork,info

1833

1835

integer :: lwork,info

1834

1836

integer :: i,j

1835

1837

integer :: i,j

1836

1838

if (present(status)) status=1

1837

1839

if (present(status)) status=1

allocate(mat(np,deg+1),bmat(np,1))

1838

1840

allocate(mat(np,deg+1),bmat(np,1))

1839

1841

! Design matrix valorisation

1840

1842

! Design matrix valorisation

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1841

1843

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1842

1844

! second member valorisation

1843

1845

! second member valorisation

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1844

1846

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1845

1847

! query workspace size

1846

1848

! query workspace size

call dgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)

1847

1849

call dgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)

lwork=twork(1)

1848

1850

lwork=twork(1)

allocate(work(int(lwork)))

1849

1851

allocate(work(int(lwork)))

! real call

1850

1852

! real call

call dgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)

1851

1853

call dgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)

1852

1854

if (info /= 0) then

1853

1855

if (info /= 0) then

if (present(status)) status=0

1854

1856

if (present(status)) status=0

end if

1855

1857

end if

1856

1858

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1857

1859

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1858

1860

deallocate(work)

1859

1861

deallocate(work)

deallocate(mat,bmat)

1860

1862

deallocate(mat,bmat)

end subroutine

1861

1863

end subroutine

1862

1864

subroutine fvn_s_lspoly(np,x,y,deg,coeff,status)

1863

1865

subroutine fvn_s_lspoly(np,x,y,deg,coeff,status)

!

1864

1866

!

! Least square polynomial fitting

1865

1867

! Least square polynomial fitting

!

1866

1868

!

! Find the coefficients of the least square polynomial of a given degree

1867

1869

! Find the coefficients of the least square polynomial of a given degree

! for a set of coordinates.

1868

1870

! for a set of coordinates.

!

1869

1871

!

! The degree must be lower than the number of points

1870

1872

! The degree must be lower than the number of points

!

1871

1873

!

! np (in) : number of points

1872

1874

! np (in) : number of points

! x(np) (in) : x data

1873

1875

! x(np) (in) : x data

! y(np) (in) : y data

1874

1876

! y(np) (in) : y data

! deg (in) : polynomial's degree

1875

1877

! deg (in) : polynomial's degree

! coeff(deg+1) (out) : polynomial coefficients

1876

1878

! coeff(deg+1) (out) : polynomial coefficients

! status (out) : =0 if a problem occurs

1877

1879

! status (out) : =0 if a problem occurs

!

1878

1880

!

implicit none

1879

1881

implicit none

1880

1882

integer, intent(in) :: np,deg

1881

1883

integer, intent(in) :: np,deg

real(kind=4), intent(in), dimension(np) :: x,y

1882

1884

real(kind=4), intent(in), dimension(np) :: x,y

real(kind=4), intent(out), dimension(deg+1) :: coeff

1883

1885

real(kind=4), intent(out), dimension(deg+1) :: coeff

integer, intent(out), optional :: status

1884

1886

integer, intent(out), optional :: status

1885

1887

real(kind=4), allocatable, dimension(:,:) :: mat,bmat

1886

1888

real(kind=4), allocatable, dimension(:,:) :: mat,bmat

real(kind=4),dimension(:),allocatable :: work

1887

1889

real(kind=4),dimension(:),allocatable :: work

real(kind=4),dimension(1) :: twork

1888

1890

real(kind=4),dimension(1) :: twork

integer :: lwork,info

1889

1891

integer :: lwork,info

1890

1892

integer :: i,j

1891

1893

integer :: i,j

1892

1894

if (present(status)) status=1

1893

1895

if (present(status)) status=1

allocate(mat(np,deg+1),bmat(np,1))

1894

1896

allocate(mat(np,deg+1),bmat(np,1))

1895

1897

! Design matrix valorisation

1896

1898

! Design matrix valorisation

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1897

1899

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1898

1900

! second member valorisation

1899

1901

! second member valorisation

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1900

1902

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1901

1903

! query workspace size

1902

1904

! query workspace size

call sgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)

1903

1905

call sgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)

lwork=twork(1)

1904

1906

lwork=twork(1)

allocate(work(int(lwork)))

1905

1907

allocate(work(int(lwork)))

! real call

1906

1908

! real call

call sgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)

1907

1909

call sgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)

1908

1910

if (info /= 0) then

1909

1911

if (info /= 0) then

if (present(status)) status=0

1910

1912

if (present(status)) status=0

end if

1911

1913

end if

1912

1914

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1913

1915

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1914

1916

deallocate(work)

1915

1917

deallocate(work)

deallocate(mat,bmat)

1916

1918

deallocate(mat,bmat)

end subroutine

1917

1919

end subroutine

1918

1920

1919

1921

1920

1922

1921

1923

1922

1924

1923

1925

1924

1926

1925

1927

subroutine fvn_d_lspoly_svd(np,x,y,deg,coeff,status)

1926

1928

subroutine fvn_d_lspoly_svd(np,x,y,deg,coeff,status)

!

1927

1929

!

! Least square polynomial fitting using singular value decomposition

1928

1930

! Least square polynomial fitting using singular value decomposition

!

1929

1931

!

! Find the coefficients of the least square polynomial of a given degree

1930

1932

! Find the coefficients of the least square polynomial of a given degree

! for a set of coordinates.

1931

1933

! for a set of coordinates.

!

1932

1934

!

! The degree must be lower than the number of points

1933

1935

! The degree must be lower than the number of points

!

1934

1936

!

! np (in) : number of points

1935

1937

! np (in) : number of points

! x(np) (in) : x data

1936

1938

! x(np) (in) : x data

! y(np) (in) : y data

1937

1939

! y(np) (in) : y data

! deg (in) : polynomial's degree

1938

1940

! deg (in) : polynomial's degree

! coeff(deg+1) (out) : polynomial coefficients

1939

1941

! coeff(deg+1) (out) : polynomial coefficients

! status (out) : =0 if a problem occurs

1940

1942

! status (out) : =0 if a problem occurs

!

1941

1943

!

implicit none

1942

1944

implicit none

1943

1945

integer, intent(in) :: np,deg

1944

1946

integer, intent(in) :: np,deg

real(kind=8), intent(in), dimension(np) :: x,y

1945

1947

real(kind=8), intent(in), dimension(np) :: x,y

real(kind=8), intent(out), dimension(deg+1) :: coeff

1946

1948

real(kind=8), intent(out), dimension(deg+1) :: coeff

integer, intent(out), optional :: status

1947

1949

integer, intent(out), optional :: status

1948

1950

real(kind=8), allocatable, dimension(:,:) :: mat,bmat

1949

1951

real(kind=8), allocatable, dimension(:,:) :: mat,bmat

real(kind=8),dimension(:),allocatable :: work,singval

1950

1952

real(kind=8),dimension(:),allocatable :: work,singval

real(kind=8),dimension(1) :: twork

1951

1953

real(kind=8),dimension(1) :: twork

integer :: lwork,info,rank

1952

1954

integer :: lwork,info,rank

1953

1955

integer :: i,j

1954

1956

integer :: i,j

1955

1957

if (present(status)) status=1

1956

1958

if (present(status)) status=1

allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))

1957

1959

allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))

1958

1960

! Design matrix valorisation

1959

1961

! Design matrix valorisation

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1960

1962

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1961

1963

! second member valorisation

1962

1964

! second member valorisation

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1963

1965

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1964

1966

! query workspace size

1965

1967

! query workspace size

call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)

1966

1968

call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)

lwork=twork(1)

1967

1969

lwork=twork(1)

allocate(work(int(lwork)))

1968

1970

allocate(work(int(lwork)))

! real call

1969

1971

! real call

call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)

1970

1972

call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)

1971

1973

if (info /= 0) then

1972

1974

if (info /= 0) then

if (present(status)) status=0

1973

1975

if (present(status)) status=0

end if

1974

1976

end if

1975

1977

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1976

1978

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1977

1979

deallocate(work)

1978

1980

deallocate(work)

deallocate(mat,bmat,singval)

1979

1981

deallocate(mat,bmat,singval)

end subroutine

1980

1982

end subroutine

1981

1983

subroutine fvn_s_lspoly_svd(np,x,y,deg,coeff,status)

1982

1984

subroutine fvn_s_lspoly_svd(np,x,y,deg,coeff,status)

!

1983

1985

!

! Least square polynomial fitting using singular value decomposition

1984

1986

! Least square polynomial fitting using singular value decomposition

!

1985

1987

!

! Find the coefficients of the least square polynomial of a given degree

1986

1988

! Find the coefficients of the least square polynomial of a given degree

! for a set of coordinates.

1987

1989

! for a set of coordinates.

!

1988

1990

!

! The degree must be lower than the number of points

1989

1991

! The degree must be lower than the number of points

!

1990

1992

!

! np (in) : number of points

1991

1993

! np (in) : number of points

! x(np) (in) : x data

1992

1994

! x(np) (in) : x data

! y(np) (in) : y data

1993

1995

! y(np) (in) : y data

! deg (in) : polynomial's degree

1994

1996

! deg (in) : polynomial's degree

! coeff(deg+1) (out) : polynomial coefficients

1995

1997

! coeff(deg+1) (out) : polynomial coefficients

! status (out) : =0 if a problem occurs

1996

1998

! status (out) : =0 if a problem occurs

!

1997

1999

!

implicit none

1998

2000

implicit none

1999

2001

integer, intent(in) :: np,deg

2000

2002

integer, intent(in) :: np,deg

real(kind=4), intent(in), dimension(np) :: x,y

2001

2003

real(kind=4), intent(in), dimension(np) :: x,y

real(kind=4), intent(out), dimension(deg+1) :: coeff

2002

2004

real(kind=4), intent(out), dimension(deg+1) :: coeff

integer, intent(out), optional :: status

2003

2005

integer, intent(out), optional :: status

2004

2006

real(kind=4), allocatable, dimension(:,:) :: mat,bmat

2005

2007

real(kind=4), allocatable, dimension(:,:) :: mat,bmat

real(kind=4),dimension(:),allocatable :: work,singval

2006

2008

real(kind=4),dimension(:),allocatable :: work,singval

real(kind=4),dimension(1) :: twork

2007

2009

real(kind=4),dimension(1) :: twork

integer :: lwork,info,rank

2008

2010

integer :: lwork,info,rank

2009

2011

integer :: i,j

2010

2012

integer :: i,j

2011

2013

if (present(status)) status=1

2012

2014

if (present(status)) status=1

allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))

2013

2015

allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))

2014

2016

! Design matrix valorisation

2015

2017

! Design matrix valorisation

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

2016

2018

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

2017

2019

! second member valorisation

2018

2020

! second member valorisation

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

2019

2021

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

2020

2022

! query workspace size

2021

2023

! query workspace size

call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)

2022

2024

call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)

lwork=twork(1)

2023

2025

lwork=twork(1)

allocate(work(int(lwork)))

2024

2026

allocate(work(int(lwork)))

! real call

2025

2027

! real call

call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)

2026

2028

call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)

2027

2029

if (info /= 0) then

2028

2030

if (info /= 0) then

if (present(status)) status=0

2029

2031

if (present(status)) status=0

end if

2030

2032

end if

2031

2033

coeff = (/ (bmat(i,1),i=1,deg+1) /)

2032

2034

coeff = (/ (bmat(i,1),i=1,deg+1) /)

2033

2035

deallocate(work)

2034

2036

deallocate(work)

deallocate(mat,bmat,singval)

2035

2037

deallocate(mat,bmat,singval)

end subroutine

2036

2038

end subroutine

2037

2039

2038

2040

!

2039

2041

!

! Muller

2040

2042

! Muller

!

2041

2043

!

2042

2044

!

2043

2045

!

! William Daniau 2007

2044

2046

! William Daniau 2007

!

2045

2047

!

! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f

2046

2048

! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f

! http://plato.asu.edu/ftp/other_software/muller.f

2047

2049

! http://plato.asu.edu/ftp/other_software/muller.f

!

2048

2050

!

! it can be used as a replacement for imsl routine dzanly with minor changes

2049

2051

! it can be used as a replacement for imsl routine dzanly with minor changes

!

2050

2052

!

!-----------------------------------------------------------------------

2051

2053

!-----------------------------------------------------------------------

!

2052

2054

!

! purpose - zeros of an analytic complex function

2053

2055

! purpose - zeros of an analytic complex function

! using the muller method with deflation

2054

2056

! using the muller method with deflation

!

2055

2057

!

! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,

2056

2058

! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,

! infer,ier)

2057

2059

! infer,ier)

!

2058

2060

!

! arguments f - a complex function subprogram, f(z), written

2059

2061

! arguments f - a complex function subprogram, f(z), written

! by the user specifying the equation whose

2060

2062

! by the user specifying the equation whose

! roots are to be found. f must appear in

2061

2063

! roots are to be found. f must appear in

! an external statement in the calling pro-

2062

2064

! an external statement in the calling pro-

! gram.

2063

2065

! gram.

! eps - 1st stopping criterion. let fp(z)=f(z)/p

2064

2066

! eps - 1st stopping criterion. let fp(z)=f(z)/p

! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))

2065

2067

! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))

! and z(1),...,z(k-1) are previously found

2066

2068

! and z(1),...,z(k-1) are previously found

! roots. if ((cdabs(f(z)).le.eps) .and.

2067

2069

! roots. if ((cdabs(f(z)).le.eps) .and.

! (cdabs(fp(z)).le.eps)), then z is accepted

2068

2070

! (cdabs(fp(z)).le.eps)), then z is accepted

! as a root. (input)

2069

2071

! as a root. (input)

! eps1 - 2nd stopping criterion. a root is accepted

2070

2072

! eps1 - 2nd stopping criterion. a root is accepted

! if two successive approximations to a given

2071

2073

! if two successive approximations to a given

! root agree within eps1. (input)

2072

2074

! root agree within eps1. (input)

! note. if either or both of the stopping

2073

2075

! note. if either or both of the stopping

! criteria are fulfilled, the root is

2074

2076

! criteria are fulfilled, the root is

! accepted.

2075

2077

! accepted.

! kn - the number of known roots which must be stored

2076

2078

! kn - the number of known roots which must be stored

! in x(1),...,x(kn), prior to entry to muller

2077

2079

! in x(1),...,x(kn), prior to entry to muller

! nguess - the number of initial guesses provided. these

2078

2080

! nguess - the number of initial guesses provided. these

! guesses must be stored in x(kn+1),...,

2079

2081

! guesses must be stored in x(kn+1),...,

! x(kn+nguess). nguess must be set equal

2080

2082

! x(kn+nguess). nguess must be set equal

! to zero if no guesses are provided. (input)

2081

2083

! to zero if no guesses are provided. (input)

! n - the number of new roots to be found by

2082

2084

! n - the number of new roots to be found by

! muller (input)

2083

2085

! muller (input)

! x - a complex vector of length kn+n. x(1),...,

2084

2086

! x - a complex vector of length kn+n. x(1),...,

! x(kn) on input must contain any known

2085

2087

! x(kn) on input must contain any known

! roots. x(kn+1),..., x(kn+n) on input may,

2086

2088

! roots. x(kn+1),..., x(kn+n) on input may,

! on user option, contain initial guesses for

2087

2089

! on user option, contain initial guesses for

! the n new roots which are to be computed.

2088

2090

! the n new roots which are to be computed.

! if the user does not provide an initial

2089

2091

! if the user does not provide an initial

! guess, zero is used.

2090

2092

! guess, zero is used.

! on output, x(kn+1),...,x(kn+n) contain the

2091

2093

! on output, x(kn+1),...,x(kn+n) contain the

! approximate roots found by muller.

2092

2094

! approximate roots found by muller.

! itmax - the maximum allowable number of iterations

2093

2095

! itmax - the maximum allowable number of iterations

! per root (input)

2094

2096

! per root (input)

! infer - an integer vector of length kn+n. on

2095

2097

! infer - an integer vector of length kn+n. on

! output infer(j) contains the number of

2096

2098

! output infer(j) contains the number of

! iterations used in finding the j-th root

2097

2099

! iterations used in finding the j-th root

! when convergence was achieved. if

2098

2100

! when convergence was achieved. if

! convergence was not obtained in itmax

2099

2101

! convergence was not obtained in itmax

! iterations, infer(j) will be greater than

2100

2102

! iterations, infer(j) will be greater than

! itmax (output).

2101

2103

! itmax (output).

! ier - error parameter (output)

2102

2104

! ier - error parameter (output)

! warning error

2103

2105

! warning error

! ier = 33 indicates failure to converge with-

2104

2106

! ier = 33 indicates failure to converge with-

! in itmax iterations for at least one of

2105

2107

! in itmax iterations for at least one of

! the (n) new roots.

2106

2108

! the (n) new roots.

!

2107

2109

!

2108

2110

!

! remarks muller always returns the last approximation for root j

2109

2111

! remarks muller always returns the last approximation for root j

! in x(j). if the convergence criterion is satisfied,

2110

2112

! in x(j). if the convergence criterion is satisfied,

! then infer(j) is less than or equal to itmax. if the

2111

2113

! then infer(j) is less than or equal to itmax. if the

! convergence criterion is not satisified, then infer(j)

2112

2114

! convergence criterion is not satisified, then infer(j)

! is set to either itmax+1 or itmax+k, with k greater

2113

2115

! is set to either itmax+1 or itmax+k, with k greater

! than 1. infer(j) = itmax+1 indicates that muller did

2114

2116

! than 1. infer(j) = itmax+1 indicates that muller did

! not obtain convergence in the allowed number of iter-

2115

2117

! not obtain convergence in the allowed number of iter-

! ations. in this case, the user may wish to set itmax

2116

2118

! ations. in this case, the user may wish to set itmax

! to a larger value. infer(j) = itmax+k means that con-

2117

2119

! to a larger value. infer(j) = itmax+k means that con-

! vergence was obtained (on iteration k) for the defla-

2118

2120

! vergence was obtained (on iteration k) for the defla-

! ted function

2119

2121

! ted function

! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))

2120

2122

! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))

!

2121

2123

!

! but failed for f(z). in this case, better initial

2122

2124

! but failed for f(z). in this case, better initial

! guesses might help or, it might be necessary to relax

2123

2125

! guesses might help or, it might be necessary to relax

! the convergence criterion.

2124

2126

! the convergence criterion.

!

2125

2127

!

!-----------------------------------------------------------------------

2126

2128

!-----------------------------------------------------------------------

!

2127

2129

!

subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)

2128

2130

subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)

implicit none

2129

2131

implicit none

double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq,eps1w

2130

2132

double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq,eps1w

double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &

2131

2133

double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &

tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &

2132

2134

tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &

zero,p1,one,four,p5

2133

2135

zero,p1,one,four,p5

2134

2136

double complex, external :: f

2135

2137

double complex, external :: f

integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &

2136

2138

integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &

knpng,jk,ick,nn,lm1,errcode

2137

2139

knpng,jk,ick,nn,lm1,errcode

double complex :: x(kn+n)

2138

2140

double complex :: x(kn+n)

integer :: infer(kn+n)

2139

2141

integer :: infer(kn+n)

2140

2142

2141

2143

data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &

2142

2144

data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &

one/(1.0,0.0)/,four/(4.0,0.0)/, &

2143

2145

one/(1.0,0.0)/,four/(4.0,0.0)/, &

p5/(0.5,0.0)/, &

2144

2146

p5/(0.5,0.0)/, &

rzero/0.0/,rten/10.0/,rhun/100.0/, &

2145

2147

rzero/0.0/,rten/10.0/,rhun/100.0/, &

ax/0.1/,ickmax/3/,rp01/0.01/

2146

2148

ax/0.1/,ickmax/3/,rp01/0.01/

2147

2149

ier = 0

2148

2150

ier = 0

if (n .lt. 1) then ! What the hell are doing here then ...

2149

2151

if (n .lt. 1) then ! What the hell are doing here then ...

return

2150

2152

return

end if

2151

2153

end if

!eps1 = rten **(-nsig)

2152

2154

!eps1 = rten **(-nsig)

eps1w = min(eps1,rp01)

2153

2155

eps1w = min(eps1,rp01)

2154

2156

knp1 = kn+1

2155

2157

knp1 = kn+1

knpn = kn+n

2156

2158

knpn = kn+n

knpng = kn+nguess

2157

2159

knpng = kn+nguess

do i=1,knpn

2158

2160

do i=1,knpn

infer(i) = 0

2159

2161

infer(i) = 0

if (i .gt. knpng) x(i) = zero

2160

2162

if (i .gt. knpng) x(i) = zero

end do

2161

2163

end do

l= knp1

2162

2164

l= knp1

2163

2165

ic=0

2164

2166

ic=0

rloop: do while (l<=knpn) ! Main loop over new roots

2165

2167

rloop: do while (l<=knpn) ! Main loop over new roots

jk = 0

2166

2168

jk = 0

ick = 0

2167

2169

ick = 0

xl = x(l)

2168

2170

xl = x(l)

icloop: do

2169

2171

icloop: do

ic = 0

2170

2172

ic = 0

h = ax

2171

2173

h = ax

h = p1*h

2172

2174

h = p1*h

if (cdabs(xl) .gt. ax) h = p1*xl

2173

2175

if (cdabs(xl) .gt. ax) h = p1*xl

! first three points are

2174

2176

! first three points are

! xl+h, xl-h, xl

2175

2177

! xl+h, xl-h, xl

rt = xl+h

2176

2178

rt = xl+h

call deflated_work(errcode)

2177

2179

call deflated_work(errcode)

if (errcode == 1) then

2178

2180

if (errcode == 1) then

exit icloop

2179

2181

exit icloop

end if

2180

2182

end if

2181

2183

z0 = fprt

2182

2184

z0 = fprt

y0 = frt

2183

2185

y0 = frt

x0 = rt

2184

2186

x0 = rt

rt = xl-h

2185

2187

rt = xl-h

call deflated_work(errcode)

2186

2188

call deflated_work(errcode)

if (errcode == 1) then

2187

2189

if (errcode == 1) then

exit icloop

2188

2190

exit icloop

end if

2189

2191

end if

2190

2192

z1 = fprt

2191

2193

z1 = fprt

y1 = frt

2192

2194

y1 = frt

h = xl-rt

2193

2195

h = xl-rt

d = h/(rt-x0)

2194

2196

d = h/(rt-x0)

rt = xl

2195

2197

rt = xl

2196

2198

call deflated_work(errcode)

2197

2199

call deflated_work(errcode)

if (errcode == 1) then

2198

2200

if (errcode == 1) then

exit icloop

2199

2201

exit icloop

end if

2200

2202

end if

2201

2203

2202

2204

z2 = fprt

2203

2205

z2 = fprt

y2 = frt

2204

2206

y2 = frt

! begin main algorithm

2205

2207

! begin main algorithm

iloop: do

2206

2208

iloop: do

dd = one + d

2207

2209

dd = one + d

t1 = z0*d*d

2208

2210

t1 = z0*d*d

t2 = z1*dd*dd

2209

2211

t2 = z1*dd*dd

xx = z2*dd

2210

2212

xx = z2*dd

t3 = z2*d

2211

2213

t3 = z2*d

bi = t1-t2+xx+t3

2212

2214

bi = t1-t2+xx+t3

den = bi*bi-four*(xx*t1-t3*(t2-xx))

2213

2215

den = bi*bi-four*(xx*t1-t3*(t2-xx))

! use denominator of maximum amplitude

2214

2216

! use denominator of maximum amplitude

t1 = cdsqrt(den)

2215

2217

t1 = cdsqrt(den)

qz = rhun*max(cdabs(bi),cdabs(t1))

2216

2218

qz = rhun*max(cdabs(bi),cdabs(t1))

t2 = bi + t1

2217

2219

t2 = bi + t1

tpq = cdabs(t2)+qz

2218

2220

tpq = cdabs(t2)+qz

if (tpq .eq. qz) t2 = zero

2219

2221

if (tpq .eq. qz) t2 = zero

t3 = bi - t1

2220

2222

t3 = bi - t1

tpq = cdabs(t3) + qz

2221

2223

tpq = cdabs(t3) + qz

if (tpq .eq. qz) t3 = zero

2222

2224

if (tpq .eq. qz) t3 = zero

den = t2

2223

2225

den = t2

qz = cdabs(t3)-cdabs(t2)

2224

2226

qz = cdabs(t3)-cdabs(t2)

if (qz .gt. rzero) den = t3

2225

2227

if (qz .gt. rzero) den = t3

! test for zero denominator

2226

2228

! test for zero denominator

if (cdabs(den) .eq. rzero) then

2227

2229

if (cdabs(den) .eq. rzero) then

call trans_rt()

2228

2230

call trans_rt()

call deflated_work(errcode)

2229

2231

call deflated_work(errcode)

if (errcode == 1) then

2230

2232

if (errcode == 1) then

exit icloop

2231

2233

exit icloop

end if

2232

2234

end if

z2 = fprt

2233

2235

z2 = fprt

y2 = frt

2234

2236

y2 = frt

cycle iloop

2235

2237

cycle iloop

end if

2236

2238

end if

2237

2239

2238

2240

d = -xx/den

2239

2241

d = -xx/den

d = d+d

2240

2242

d = d+d

h = d*h

2241

2243

h = d*h

rt = rt + h

2242

2244

rt = rt + h

! check convergence of the first kind

2243

2245

! check convergence of the first kind

if (cdabs(h) .le. eps1w*max(cdabs(rt),ax)) then

2244

2246

if (cdabs(h) .le. eps1w*max(cdabs(rt),ax)) then

if (ic .ne. 0) then

2245

2247

if (ic .ne. 0) then

exit icloop

2246

2248

exit icloop

end if

2247

2249

end if

ic = 1

2248

2250

ic = 1

z0 = y1

2249

2251

z0 = y1

z1 = y2

2250

2252

z1 = y2

z2 = f(rt)

2251

2253

z2 = f(rt)

xl = rt

2252

2254

xl = rt

ick = ick+1

2253

2255

ick = ick+1

if (ick .le. ickmax) then

2254

2256

if (ick .le. ickmax) then

cycle iloop

2255

2257

cycle iloop

end if

2256

2258

end if

! warning error, itmax = maximum

2257

2259

! warning error, itmax = maximum

jk = itmax + jk

2258

2260

jk = itmax + jk

ier = 33

2259

2261

ier = 33

end if

2260

2262

end if

if (ic .ne. 0) then

2261

2263

if (ic .ne. 0) then

cycle icloop

2262

2264

cycle icloop

end if

2263

2265

end if

call deflated_work(errcode)

2264

2266

call deflated_work(errcode)

if (errcode == 1) then

2265

2267

if (errcode == 1) then

exit icloop

2266

2268

exit icloop

end if

2267

2269

end if

2268

2270

do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)

2269

2271

do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)

! take remedial action to induce

2270

2272

! take remedial action to induce

! convergence

2271

2273

! convergence

d = d*p5

2272

2274

d = d*p5

h = h*p5

2273

2275

h = h*p5

rt = rt-h

2274

2276

rt = rt-h

call deflated_work(errcode)

2275

2277

call deflated_work(errcode)

if (errcode == 1) then

2276

2278

if (errcode == 1) then

exit icloop

2277

2279

exit icloop

end if

2278

2280

end if

end do

2279

2281

end do

z0 = z1

2280

2282

z0 = z1

z1 = z2

2281

2283

z1 = z2

z2 = fprt

2282

2284

z2 = fprt

y0 = y1

2283

2285

y0 = y1

y1 = y2

2284

2286

y1 = y2

y2 = frt

2285

2287

y2 = frt

end do iloop

2286

2288

end do iloop

end do icloop

2287

2289

end do icloop

x(l) = rt

2288

2290

x(l) = rt

infer(l) = jk

2289

2291

infer(l) = jk

l = l+1

2290

2292

l = l+1

end do rloop

2291

2293

end do rloop

2292

2294

contains

2293

2295

contains

subroutine trans_rt()

2294

2296

subroutine trans_rt()

tem = rten*eps1w

2295

2297

tem = rten*eps1w

if (cdabs(rt) .gt. ax) tem = tem*rt

2296

2298

if (cdabs(rt) .gt. ax) tem = tem*rt

rt = rt+tem

2297

2299

rt = rt+tem

d = (h+tem)*d/h

2298

2300

d = (h+tem)*d/h

h = h+tem

2299

2301

h = h+tem

end subroutine trans_rt

2300

2302

end subroutine trans_rt

2301

2303

subroutine deflated_work(errcode)

2302

2304

subroutine deflated_work(errcode)

! errcode=0 => no errors

2303

2305

! errcode=0 => no errors

! errcode=1 => jk>itmax or convergence of second kind achieved

2304

2306

! errcode=1 => jk>itmax or convergence of second kind achieved

integer :: errcode,flag

2305

2307

integer :: errcode,flag

2306

2308

flag=1

2307

2309

flag=1

loop1: do while(flag==1)

2308

2310

loop1: do while(flag==1)

errcode=0

2309

2311

errcode=0

jk = jk+1

2310

2312

jk = jk+1

if (jk .gt. itmax) then

2311

2313

if (jk .gt. itmax) then

ier=33

2312

2314

ier=33

errcode=1

2313

2315

errcode=1

return

2314

2316

return

end if

2315

2317

end if

frt = f(rt)

2316

2318

frt = f(rt)

fprt = frt

2317

2319

fprt = frt

if (l /= 1) then

2318

2320

if (l /= 1) then

lm1 = l-1

2319

2321

lm1 = l-1

do i=1,lm1

2320

2322

do i=1,lm1

tem = rt - x(i)

2321

2323

tem = rt - x(i)

if (cdabs(tem) .eq. rzero) then

2322

2324

if (cdabs(tem) .eq. rzero) then

!if (ic .ne. 0) go to 15 !! ?? possible?

2323

2325

!if (ic .ne. 0) go to 15 !! ?? possible?

call trans_rt()

2324

2326

call trans_rt()

cycle loop1

2325

2327

cycle loop1

end if

2326

2328

end if

fprt = fprt/tem

2327

2329

fprt = fprt/tem

end do

2328

2330

end do

end if

2329

2331

end if

flag=0

2330

2332

flag=0

end do loop1

2331

2333

end do loop1

2332

2334

if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then

2333

2335

if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then

errcode=1

2334

2336

errcode=1

return

2335

2337

return

end if

2336

2338

end if

2337

2339

end subroutine deflated_work

2338

2340

end subroutine deflated_work

2339

2341

end subroutine

2340

2342

end subroutine

2341

2343

2342

2344

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

2343

2345

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

2344

2346

!

! Integration

2345

2347

! Integration

!

2346

2348

!

! Only double precision coded atm

2347

2349

! Only double precision coded atm

!

2348

2350

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

2349

2351

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

2350

2352

2351

2353

subroutine fvn_d_gauss_legendre(n,qx,qw)

2352

2354

subroutine fvn_d_gauss_legendre(n,qx,qw)

!

2353

2355

!

! This routine compute the n Gauss Legendre abscissas and weights

2354

2356

! This routine compute the n Gauss Legendre abscissas and weights

! Adapted from Numerical Recipes routine gauleg

2355

2357

! Adapted from Numerical Recipes routine gauleg

!

2356

2358

!

! n (in) : number of points

2357

2359

! n (in) : number of points

! qx(out) : abscissas

2358

2360

! qx(out) : abscissas

! qw(out) : weights

2359

2361

! qw(out) : weights

!

2360

2362

!

implicit none

2361

2363

implicit none

double precision,parameter :: pi=3.141592653589793d0

2362

2364

double precision,parameter :: pi=3.141592653589793d0

integer, intent(in) :: n

2363

2365

integer, intent(in) :: n

double precision, intent(out) :: qx(n),qw(n)

2364

2366

double precision, intent(out) :: qx(n),qw(n)

2365

2367

integer :: m,i,j

2366

2368

integer :: m,i,j

double precision :: z,z1,p1,p2,p3,pp

2367

2369

double precision :: z,z1,p1,p2,p3,pp

m=(n+1)/2

2368

2370

m=(n+1)/2

do i=1,m

2369

2371

do i=1,m

z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))

2370

2372

z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))

iloop: do

2371

2373

iloop: do

p1=1.d0

2372

2374

p1=1.d0

p2=0.d0

2373

2375

p2=0.d0

do j=1,n

2374

2376

do j=1,n

p3=p2

2375

2377

p3=p2

p2=p1

2376

2378

p2=p1

p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)

2377

2379

p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)

end do

2378

2380

end do

pp=dble(n)*(z*p1-p2)/(z*z-1.d0)

2379

2381

pp=dble(n)*(z*p1-p2)/(z*z-1.d0)

z1=z

2380

2382

z1=z

z=z1-p1/pp

2381

2383

z=z1-p1/pp

if (dabs(z-z1)<=epsilon(z)) then

2382

2384

if (dabs(z-z1)<=epsilon(z)) then

exit iloop

2383

2385

exit iloop

end if

2384

2386

end if

end do iloop

2385

2387

end do iloop

qx(i)=-z

2386

2388

qx(i)=-z

qx(n+1-i)=z

2387

2389

qx(n+1-i)=z

qw(i)=2.d0/((1.d0-z*z)*pp*pp)

2388

2390

qw(i)=2.d0/((1.d0-z*z)*pp*pp)

qw(n+1-i)=qw(i)

2389

2391

qw(n+1-i)=qw(i)

end do

2390

2392

end do

end subroutine

2391

2393

end subroutine

2392

2394

2393

2395

2394

2396

subroutine fvn_d_gl_integ(f,a,b,n,res)

2395

2397

subroutine fvn_d_gl_integ(f,a,b,n,res)

!

2396

2398

!

! This is a simple non adaptative integration routine

2397

2399

! This is a simple non adaptative integration routine

! using n gauss legendre abscissas and weights

2398

2400

! using n gauss legendre abscissas and weights

!

2399

2401

!

! f(in) : the function to integrate

2400

2402

! f(in) : the function to integrate

! a(in) : lower bound

2401

2403

! a(in) : lower bound

! b(in) : higher bound

2402

2404

! b(in) : higher bound

! n(in) : number of gauss legendre pairs

2403

2405

! n(in) : number of gauss legendre pairs

! res(out): the evaluation of the integral

2404

2406

! res(out): the evaluation of the integral

!

2405

2407

!

double precision,external :: f

2406

2408

double precision,external :: f

double precision, intent(in) :: a,b

2407

2409

double precision, intent(in) :: a,b

integer, intent(in):: n

2408

2410

integer, intent(in):: n

double precision, intent(out) :: res

2409

2411

double precision, intent(out) :: res

2410

2412

double precision, allocatable :: qx(:),qw(:)

2411

2413

double precision, allocatable :: qx(:),qw(:)

double precision :: xm,xr

2412

2414

double precision :: xm,xr

integer :: i

2413

2415

integer :: i

2414

2416

! First compute n gauss legendre abs and weight

2415

2417

! First compute n gauss legendre abs and weight

allocate(qx(n))

2416

2418

allocate(qx(n))

allocate(qw(n))

2417

2419

allocate(qw(n))

call fvn_d_gauss_legendre(n,qx,qw)

2418

2420

call fvn_d_gauss_legendre(n,qx,qw)

2419

2421

xm=0.5d0*(b+a)

2420

2422

xm=0.5d0*(b+a)

xr=0.5d0*(b-a)

2421

2423

xr=0.5d0*(b-a)

2422

2424

res=0.d0

2423

2425

res=0.d0

2424

2426

do i=1,n

2425

2427

do i=1,n

res=res+qw(i)*f(xm+xr*qx(i))

2426

2428

res=res+qw(i)*f(xm+xr*qx(i))

end do

2427

2429

end do

2428

2430

res=xr*res

2429

2431

res=xr*res

2430

2432

deallocate(qw)

2431

2433

deallocate(qw)

deallocate(qx)

2432

2434

deallocate(qx)

2433

2435

end subroutine

2434

2436

end subroutine

2435

2437

!!!!!!!!!!!!!!!!!!!!!!!!

2436

2438

!!!!!!!!!!!!!!!!!!!!!!!!

!

2437

2439

!

! Simple and double adaptative Gauss Kronrod integration based on

2438

2440

! Simple and double adaptative Gauss Kronrod integration based on

! a modified version of quadpack ( http://www.netlib.org/quadpack

2439

2441

! a modified version of quadpack ( http://www.netlib.org/quadpack

!

2440

2442

!

! Common parameters :

2441

2443

! Common parameters :

!

2442

2444

!

! key (in)

2443

2445

! key (in)

! epsabs

2444

2446

! epsabs

! epsrel

2445

2447

! epsrel

!

2446

2448

!

2447

2449

!

!!!!!!!!!!!!!!!!!!!!!!!!

2448

2450

!!!!!!!!!!!!!!!!!!!!!!!!

2449

2451

subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

2450

2452

subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

!

2451

2453

!

! Evaluate the integral of function f(x) between a and b

2452

2454

! Evaluate the integral of function f(x) between a and b

!

2453

2455

!

! f(in) : the function

2454

2456

! f(in) : the function

! a(in) : lower bound

2455

2457

! a(in) : lower bound

! b(in) : higher bound

2456

2458

! b(in) : higher bound

! epsabs(in) : desired absolute error

2457

2459

! epsabs(in) : desired absolute error

! epsrel(in) : desired relative error

2458

2460

! epsrel(in) : desired relative error

! key(in) : gauss kronrod rule

2459

2461

! key(in) : gauss kronrod rule

! 1: 7 - 15 points

2460

2462

! 1: 7 - 15 points

! 2: 10 - 21 points

2461

2463

! 2: 10 - 21 points

! 3: 15 - 31 points

2462

2464

! 3: 15 - 31 points

! 4: 20 - 41 points

2463

2465

! 4: 20 - 41 points

! 5: 25 - 51 points

2464

2466

! 5: 25 - 51 points

! 6: 30 - 61 points

2465

2467

! 6: 30 - 61 points

!

2466

2468

!

! limit(in) : maximum number of subintervals in the partition of the

2467

2469

! limit(in) : maximum number of subintervals in the partition of the

! given integration interval (a,b). A value of 500 will give the same

2468

2470

! given integration interval (a,b). A value of 500 will give the same

! behaviour as the imsl routine dqdag

2469

2471

! behaviour as the imsl routine dqdag

!

2470

2472

!

! res(out) : estimated integral value

2471

2473

! res(out) : estimated integral value

! abserr(out) : estimated absolute error

2472

2474

! abserr(out) : estimated absolute error

! ier(out) : error flag from quadpack routines

2473

2475

! ier(out) : error flag from quadpack routines

! 0 : no error

2474

2476

! 0 : no error

! 1 : maximum number of subdivisions allowed

2475

2477

! 1 : maximum number of subdivisions allowed

! has been achieved. one can allow more

2476

2478

! has been achieved. one can allow more

! subdivisions by increasing the value of

2477

2479

! subdivisions by increasing the value of

! limit (and taking the according dimension

2478

2480

! limit (and taking the according dimension

! adjustments into account). however, if

2479

2481

! adjustments into account). however, if

! this yield no improvement it is advised

2480

2482

! this yield no improvement it is advised

! to analyze the integrand in order to

2481

2483

! to analyze the integrand in order to

! determine the integration difficulaties.

2482

2484

! determine the integration difficulaties.

! if the position of a local difficulty can

2483

2485

! if the position of a local difficulty can

! be determined (i.e.singularity,

2484

2486

! be determined (i.e.singularity,

! discontinuity within the interval) one

2485

2487

! discontinuity within the interval) one

! will probably gain from splitting up the

2486

2488

! will probably gain from splitting up the

! interval at this point and calling the

2487

2489

! interval at this point and calling the

! integrator on the subranges. if possible,

2488

2490

! integrator on the subranges. if possible,

! an appropriate special-purpose integrator

2489

2491

! an appropriate special-purpose integrator

! should be used which is designed for

2490

2492

! should be used which is designed for

! handling the type of difficulty involved.

2491

2493

! handling the type of difficulty involved.

! 2 : the occurrence of roundoff error is

2492

2494

! 2 : the occurrence of roundoff error is

! detected, which prevents the requested

2493

2495

! detected, which prevents the requested

! tolerance from being achieved.

2494

2496

! tolerance from being achieved.

! 3 : extremely bad integrand behaviour occurs

2495

2497

! 3 : extremely bad integrand behaviour occurs

! at some points of the integration

2496

2498

! at some points of the integration

! interval.

2497

2499

! interval.

! 6 : the input is invalid, because

2498

2500

! 6 : the input is invalid, because

! (epsabs.le.0 and

2499

2501

! (epsabs.le.0 and

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

2500

2502

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

! or limit.lt.1 or lenw.lt.limit*4.

2501

2503

! or limit.lt.1 or lenw.lt.limit*4.

! result, abserr, neval, last are set

2502

2504

! result, abserr, neval, last are set

! to zero.

2503

2505

! to zero.

! except when lenw is invalid, iwork(1),

2504

2506

! except when lenw is invalid, iwork(1),

! work(limit*2+1) and work(limit*3+1) are

2505

2507

! work(limit*2+1) and work(limit*3+1) are

! set to zero, work(1) is set to a and

2506

2508

! set to zero, work(1) is set to a and

! work(limit+1) to b.

2507

2509

! work(limit+1) to b.

2508

2510

implicit none

2509

2511

implicit none

double precision, external :: f

2510

2512

double precision, external :: f

double precision, intent(in) :: a,b,epsabs,epsrel

2511

2513

double precision, intent(in) :: a,b,epsabs,epsrel

integer, intent(in) :: key

2512

2514

integer, intent(in) :: key

integer, intent(in),optional :: limit

2513

2515

integer, intent(in),optional :: limit

double precision, intent(out) :: res,abserr

2514

2516

double precision, intent(out) :: res,abserr

integer, intent(out) :: ier

2515

2517

integer, intent(out) :: ier

2516

2518

double precision, allocatable :: work(:)

2517

2519

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

2518

2520

integer, allocatable :: iwork(:)

integer :: lenw,neval,last

2519

2521

integer :: lenw,neval,last

integer :: limitw

2520

2522

integer :: limitw

2521

2523

! imsl value for limit is 500

2522

2524

! imsl value for limit is 500

limitw=500

2523

2525

limitw=500

if (present(limit)) limitw=limit

2524

2526

if (present(limit)) limitw=limit

2525

2527

lenw=limitw*4

2526

2528

lenw=limitw*4

2527

2529

allocate(iwork(limitw))

2528

2530

allocate(iwork(limitw))

allocate(work(lenw))

2529

2531

allocate(work(lenw))

2530

2532

call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limitw,lenw,last,iwork,work)

2531

2533

call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limitw,lenw,last,iwork,work)

2532

2534

deallocate(work)

2533

2535

deallocate(work)

deallocate(iwork)

2534

2536

deallocate(iwork)

2535

2537

end subroutine

2536

2538

end subroutine

2537

2539

2538

2540

2539

2541

subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)

2540

2542

subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)

!

2541

2543

!

! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)

2542

2544

! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)

!

2543

2545

!

! f(in) : the function

2544

2546

! f(in) : the function

! a(in) : lower bound

2545

2547

! a(in) : lower bound

! b(in) : higher bound

2546

2548

! b(in) : higher bound

! g(in) : external function describing lower bound for y

2547

2549

! g(in) : external function describing lower bound for y

! h(in) : external function describing higher bound for y

2548

2550

! h(in) : external function describing higher bound for y

! epsabs(in) : desired absolute error

2549

2551

! epsabs(in) : desired absolute error

! epsrel(in) : desired relative error

2550

2552

! epsrel(in) : desired relative error

! key(in) : gauss kronrod rule

2551

2553

! key(in) : gauss kronrod rule

! 1: 7 - 15 points

2552

2554

! 1: 7 - 15 points

! 2: 10 - 21 points

2553

2555

! 2: 10 - 21 points

! 3: 15 - 31 points

2554

2556

! 3: 15 - 31 points

! 4: 20 - 41 points

2555

2557

! 4: 20 - 41 points

! 5: 25 - 51 points

2556

2558

! 5: 25 - 51 points

! 6: 30 - 61 points

2557

2559

! 6: 30 - 61 points

!

2558

2560

!

! limit(in) : maximum number of subintervals in the partition of the

2559

2561

! limit(in) : maximum number of subintervals in the partition of the

! given integration interval (a,b). A value of 500 will give the same

2560

2562

! given integration interval (a,b). A value of 500 will give the same

! behaviour as the imsl routine dqdag

2561

2563

! behaviour as the imsl routine dqdag

!

2562

2564

!

! res(out) : estimated integral value

2563

2565

! res(out) : estimated integral value

! abserr(out) : estimated absolute error

2564

2566

! abserr(out) : estimated absolute error

! ier(out) : error flag from quadpack routines

2565

2567

! ier(out) : error flag from quadpack routines

! 0 : no error

2566

2568

! 0 : no error

! 1 : maximum number of subdivisions allowed

2567

2569

! 1 : maximum number of subdivisions allowed

! has been achieved. one can allow more

2568

2570

! has been achieved. one can allow more

! subdivisions by increasing the value of

2569

2571

! subdivisions by increasing the value of

! limit (and taking the according dimension

2570

2572

! limit (and taking the according dimension

! adjustments into account). however, if

2571

2573

! adjustments into account). however, if

! this yield no improvement it is advised

2572

2574

! this yield no improvement it is advised

! to analyze the integrand in order to

2573

2575

! to analyze the integrand in order to

! determine the integration difficulaties.

2574

2576

! determine the integration difficulaties.

! if the position of a local difficulty can

2575

2577

! if the position of a local difficulty can

! be determined (i.e.singularity,

2576

2578

! be determined (i.e.singularity,

! discontinuity within the interval) one

2577

2579

! discontinuity within the interval) one

! will probably gain from splitting up the

2578

2580

! will probably gain from splitting up the

! interval at this point and calling the

2579

2581

! interval at this point and calling the

! integrator on the subranges. if possible,

2580

2582

! integrator on the subranges. if possible,

! an appropriate special-purpose integrator

2581

2583

! an appropriate special-purpose integrator

! should be used which is designed for

2582

2584

! should be used which is designed for

! handling the type of difficulty involved.

2583

2585

! handling the type of difficulty involved.

! 2 : the occurrence of roundoff error is

2584

2586

! 2 : the occurrence of roundoff error is

! detected, which prevents the requested

2585

2587

! detected, which prevents the requested

! tolerance from being achieved.

2586

2588

! tolerance from being achieved.

! 3 : extremely bad integrand behaviour occurs

2587

2589

! 3 : extremely bad integrand behaviour occurs

! at some points of the integration

2588

2590

! at some points of the integration

! interval.

2589

2591

! interval.

! 6 : the input is invalid, because

2590

2592

! 6 : the input is invalid, because

! (epsabs.le.0 and

2591

2593

! (epsabs.le.0 and

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

2592

2594

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

! or limit.lt.1 or lenw.lt.limit*4.

2593

2595

! or limit.lt.1 or lenw.lt.limit*4.

! result, abserr, neval, last are set

2594

2596

! result, abserr, neval, last are set

! to zero.

2595

2597

! to zero.

! except when lenw is invalid, iwork(1),

2596

2598

! except when lenw is invalid, iwork(1),

! work(limit*2+1) and work(limit*3+1) are

2597

2599

! work(limit*2+1) and work(limit*3+1) are

! set to zero, work(1) is set to a and

2598

2600

! set to zero, work(1) is set to a and

! work(limit+1) to b.

2599

2601

! work(limit+1) to b.

2600

2602

implicit none

2601

2603

implicit none

double precision, external:: f,g,h

2602

2604

double precision, external:: f,g,h

double precision, intent(in) :: a,b,epsabs,epsrel

2603

2605

double precision, intent(in) :: a,b,epsabs,epsrel

integer, intent(in) :: key

2604

2606

integer, intent(in) :: key

integer, intent(in), optional :: limit

2605

2607

integer, intent(in), optional :: limit

integer, intent(out) :: ier

2606

2608

integer, intent(out) :: ier

double precision, intent(out) :: res,abserr

2607

2609

double precision, intent(out) :: res,abserr

2608

2610

2609

2611

double precision, allocatable :: work(:)

2610

2612

double precision, allocatable :: work(:)

integer :: limitw

2611

2613

integer :: limitw

integer, allocatable :: iwork(:)

2612

2614

integer, allocatable :: iwork(:)

integer :: lenw,neval,last

2613

2615

integer :: lenw,neval,last

2614

2616

! imsl value for limit is 500

2615

2617

! imsl value for limit is 500

limitw=500

2616

2618

limitw=500

if (present(limit)) limitw=limit

2617

2619

if (present(limit)) limitw=limit

2618

2620

lenw=limitw*4

2619

2621

lenw=limitw*4

allocate(work(lenw))

2620

2622

allocate(work(lenw))

allocate(iwork(limitw))

2621

2623

allocate(iwork(limitw))

2622

2624

call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limitw,lenw,last,iwork,work)

2623

2625

call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limitw,lenw,last,iwork,work)

2624

2626

deallocate(iwork)

2625

2627

deallocate(iwork)

deallocate(work)

2626

2628

deallocate(work)

end subroutine

2627

2629

end subroutine

2628

2630

2629

2631

2630

2632

subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

2631

2633

subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

!

2632

2634

!

! Evaluate the single integral of function f(x,y) for y between a and b with a

2633

2635

! Evaluate the single integral of function f(x,y) for y between a and b with a

! given x value

2634

2636

! given x value

!

2635

2637

!

! This function is used for the evaluation of the double integral fvn_d_integ_2_gk

2636

2638

! This function is used for the evaluation of the double integral fvn_d_integ_2_gk

!

2637

2639

!

! f(in) : the function

2638

2640

! f(in) : the function

! x(in) : x

2639

2641

! x(in) : x

! a(in) : lower bound

2640

2642

! a(in) : lower bound

! b(in) : higher bound

2641

2643

! b(in) : higher bound

! epsabs(in) : desired absolute error

2642

2644

! epsabs(in) : desired absolute error

! epsrel(in) : desired relative error

2643

2645

! epsrel(in) : desired relative error

! key(in) : gauss kronrod rule

2644

2646

! key(in) : gauss kronrod rule

! 1: 7 - 15 points

2645

2647

! 1: 7 - 15 points

! 2: 10 - 21 points

2646

2648

! 2: 10 - 21 points

! 3: 15 - 31 points

2647

2649

! 3: 15 - 31 points

! 4: 20 - 41 points

2648

2650

! 4: 20 - 41 points

! 5: 25 - 51 points

2649

2651

! 5: 25 - 51 points

! 6: 30 - 61 points

2650

2652

! 6: 30 - 61 points

!

2651

2653

!

! limit(in) : maximum number of subintervals in the partition of the

2652

2654

! limit(in) : maximum number of subintervals in the partition of the

! given integration interval (a,b). A value of 500 will give the same

2653

2655

! given integration interval (a,b). A value of 500 will give the same

! behaviour as the imsl routine dqdag

2654

2656

! behaviour as the imsl routine dqdag

!

2655

2657

!

! res(out) : estimated integral value

2656

2658

! res(out) : estimated integral value

! abserr(out) : estimated absolute error

2657

2659

! abserr(out) : estimated absolute error

! ier(out) : error flag from quadpack routines

2658

2660

! ier(out) : error flag from quadpack routines

! 0 : no error

2659

2661

! 0 : no error

! 1 : maximum number of subdivisions allowed

2660

2662

! 1 : maximum number of subdivisions allowed

! has been achieved. one can allow more

2661

2663

! has been achieved. one can allow more

! subdivisions by increasing the value of

2662

2664

! subdivisions by increasing the value of

! limit (and taking the according dimension

2663

2665

! limit (and taking the according dimension

! adjustments into account). however, if

2664

2666

! adjustments into account). however, if

! this yield no improvement it is advised

2665

2667

! this yield no improvement it is advised

! to analyze the integrand in order to

2666

2668

! to analyze the integrand in order to

! determine the integration difficulaties.

2667

2669

! determine the integration difficulaties.

! if the position of a local difficulty can

2668

2670

! if the position of a local difficulty can

! be determined (i.e.singularity,

2669

2671

! be determined (i.e.singularity,

! discontinuity within the interval) one

2670

2672

! discontinuity within the interval) one

! will probably gain from splitting up the

2671

2673

! will probably gain from splitting up the

! interval at this point and calling the

2672

2674

! interval at this point and calling the

! integrator on the subranges. if possible,

2673

2675

! integrator on the subranges. if possible,

! an appropriate special-purpose integrator

2674

2676

! an appropriate special-purpose integrator

! should be used which is designed for

2675

2677

! should be used which is designed for

! handling the type of difficulty involved.

2676

2678

! handling the type of difficulty involved.

! 2 : the occurrence of roundoff error is

2677

2679

! 2 : the occurrence of roundoff error is

! detected, which prevents the requested

2678

2680

! detected, which prevents the requested

! tolerance from being achieved.

2679

2681

! tolerance from being achieved.

! 3 : extremely bad integrand behaviour occurs

2680

2682

! 3 : extremely bad integrand behaviour occurs

! at some points of the integration

2681

2683

! at some points of the integration

! interval.

2682

2684

! interval.

! 6 : the input is invalid, because

2683

2685

! 6 : the input is invalid, because

! (epsabs.le.0 and

2684

2686

! (epsabs.le.0 and

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

2685

2687

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

! or limit.lt.1 or lenw.lt.limit*4.

2686

2688

! or limit.lt.1 or lenw.lt.limit*4.

! result, abserr, neval, last are set

2687

2689

! result, abserr, neval, last are set

! to zero.

2688

2690

! to zero.

! except when lenw is invalid, iwork(1),

2689

2691

! except when lenw is invalid, iwork(1),

! work(limit*2+1) and work(limit*3+1) are

2690

2692

! work(limit*2+1) and work(limit*3+1) are

! set to zero, work(1) is set to a and

2691

2693

! set to zero, work(1) is set to a and

! work(limit+1) to b.

2692

2694

! work(limit+1) to b.

2693

2695

implicit none

2694

2696

implicit none

double precision, external:: f

2695

2697

double precision, external:: f

double precision, intent(in) :: x,a,b,epsabs,epsrel

2696

2698

double precision, intent(in) :: x,a,b,epsabs,epsrel

integer, intent(in) :: key

2697

2699

integer, intent(in) :: key

integer, intent(in),optional :: limit

2698

2700

integer, intent(in),optional :: limit

integer, intent(out) :: ier

2699

2701

integer, intent(out) :: ier

double precision, intent(out) :: res,abserr

2700

2702

double precision, intent(out) :: res,abserr

2701

2703

2702

2704

double precision, allocatable :: work(:)

2703

2705

double precision, allocatable :: work(:)

integer :: limitw

2704

2706

integer :: limitw

integer, allocatable :: iwork(:)

2705

2707

integer, allocatable :: iwork(:)

integer :: lenw,neval,last

2706

2708

integer :: lenw,neval,last

2707

2709

! imsl value for limit is 500

2708

2710

! imsl value for limit is 500

limitw=500

2709

2711

limitw=500

if (present(limit)) limitw=limit

2710

2712

if (present(limit)) limitw=limit

2711

2713

lenw=limitw*4

2712

2714

lenw=limitw*4

allocate(work(lenw))

2713

2715

allocate(work(lenw))

allocate(iwork(limitw))

2714

2716

allocate(iwork(limitw))

2715

2717

call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limitw,lenw,last,iwork,work)

2716

2718

call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limitw,lenw,last,iwork,work)

2717

2719

deallocate(iwork)

2718

2720

deallocate(iwork)

deallocate(work)

2719

2721

deallocate(work)

end subroutine

2720

2722

end subroutine

2721

2723

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

2722

2724

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! Include the modified quadpack files

2723

2725

! Include the modified quadpack files

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

2724

2726

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

include "fvn_quadpack/dqag_2d_inner.f"

2725

2727

include "fvn_quadpack/dqag_2d_inner.f"

include "fvn_quadpack/dqk15_2d_inner.f"

2726

2728

include "fvn_quadpack/dqk15_2d_inner.f"

include "fvn_quadpack/dqk31_2d_outer.f"

2727

2729

include "fvn_quadpack/dqk31_2d_outer.f"

include "fvn_quadpack/d1mach.f"

2728

2730

include "fvn_quadpack/d1mach.f"

include "fvn_quadpack/dqk31_2d_inner.f"

2729

2731

include "fvn_quadpack/dqk31_2d_inner.f"

include "fvn_quadpack/dqage.f"

2730

2732

include "fvn_quadpack/dqage.f"

include "fvn_quadpack/dqk15.f"

2731

2733

include "fvn_quadpack/dqk15.f"

include "fvn_quadpack/dqk21.f"

2732

2734

include "fvn_quadpack/dqk21.f"

include "fvn_quadpack/dqk31.f"

2733

2735

include "fvn_quadpack/dqk31.f"

include "fvn_quadpack/dqk41.f"

2734

2736

include "fvn_quadpack/dqk41.f"

include "fvn_quadpack/dqk51.f"

2735

2737

include "fvn_quadpack/dqk51.f"

include "fvn_quadpack/dqk61.f"

2736

2738

include "fvn_quadpack/dqk61.f"

include "fvn_quadpack/dqk41_2d_outer.f"

2737

2739

include "fvn_quadpack/dqk41_2d_outer.f"

include "fvn_quadpack/dqk41_2d_inner.f"

2738

2740

include "fvn_quadpack/dqk41_2d_inner.f"

include "fvn_quadpack/dqag_2d_outer.f"

2739

2741

include "fvn_quadpack/dqag_2d_outer.f"

include "fvn_quadpack/dqpsrt.f"

2740

2742

include "fvn_quadpack/dqpsrt.f"

include "fvn_quadpack/dqag.f"

2741

2743

include "fvn_quadpack/dqag.f"

include "fvn_quadpack/dqage_2d_outer.f"

2742

2744

include "fvn_quadpack/dqage_2d_outer.f"

include "fvn_quadpack/dqage_2d_inner.f"

2743

2745

include "fvn_quadpack/dqage_2d_inner.f"

include "fvn_quadpack/dqk51_2d_outer.f"

2744

2746

include "fvn_quadpack/dqk51_2d_outer.f"

include "fvn_quadpack/dqk51_2d_inner.f"

2745

2747

include "fvn_quadpack/dqk51_2d_inner.f"

include "fvn_quadpack/dqk61_2d_outer.f"

2746

2748

include "fvn_quadpack/dqk61_2d_outer.f"

include "fvn_quadpack/dqk21_2d_outer.f"

2747

2749

include "fvn_quadpack/dqk21_2d_outer.f"

include "fvn_quadpack/dqk61_2d_inner.f"

2748

2750

include "fvn_quadpack/dqk61_2d_inner.f"

include "fvn_quadpack/dqk21_2d_inner.f"

2749

2751

include "fvn_quadpack/dqk21_2d_inner.f"

include "fvn_quadpack/dqk15_2d_outer.f"

2750

2752

include "fvn_quadpack/dqk15_2d_outer.f"

2751

2753

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

2752

2754

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

2753

2755

!

! Trigonometric functions

2754

2756

! Trigonometric functions

!

2755

2757

!

! fvn_z_acos, fvn_z_asin : complex arc cosine and sine

2756

2758

! fvn_z_acos, fvn_z_asin : complex arc cosine and sine

! fvn_d_acosh : arc cosinus hyperbolic

2757

2759

! fvn_d_acosh : arc cosinus hyperbolic

! fvn_d_asinh : arc sinus hyperbolic

2758

2760

! fvn_d_asinh : arc sinus hyperbolic

!

2759

2761

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

2760

2762

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! February 2008

2761

2763

! February 2008

! All Trigonometric functions removed due to implementation of fnlib

2762

2764

! All Trigonometric functions removed due to implementation of fnlib

2763

2765

function fvn_z_acos(z)

2764

2766

function fvn_z_acos(z)

! double complex arccos function adapted from

2765

2767

! double complex arccos function adapted from

! the c gsl library

2766

2768

! the c gsl library

! http://www.gnu.org/software/gsl/

2767

2769

! http://www.gnu.org/software/gsl/

implicit none

2768

2770

implicit none

complex(kind=8) :: fvn_z_acos

2769

2771

complex(kind=8) :: fvn_z_acos

complex(kind=8) :: z

2770

2772

complex(kind=8) :: z

real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1

2771

2773

real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1

real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8

2772

2774

real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8

complex(kind=8),parameter :: i=(0._8,1._8)

2773

2775

complex(kind=8),parameter :: i=(0._8,1._8)

real(kind=8) :: r_res,i_res

2774

2776

real(kind=8) :: r_res,i_res

2775

2777

rz=dreal(z)

2776

2778

rz=dreal(z)

iz=dimag(z)

2777

2779

iz=dimag(z)

if ( iz == 0._8 ) then

2778

2780

if ( iz == 0._8 ) then

fvn_z_acos=fvn_z_acos_real(rz)

2779

2781

fvn_z_acos=fvn_z_acos_real(rz)

return

2780

2782

return

end if

2781

2783

end if

2782

2784

x=dabs(rz)

2783

2785

x=dabs(rz)

y=dabs(iz)

2784

2786

y=dabs(iz)

r=fvn_d_hypot(x+1.,y)

2785

2787

r=fvn_d_hypot(x+1.,y)

s=fvn_d_hypot(x-1.,y)

2786

2788

s=fvn_d_hypot(x-1.,y)

a=0.5*(r + s)

2787

2789

a=0.5*(r + s)

b=x/a

2788

2790

b=x/a

y2=y*y

2789

2791

y2=y*y

2790

2792

if (b <= b_crossover) then

2791

2793

if (b <= b_crossover) then

r_res=dacos(b)

2792

2794

r_res=dacos(b)

else

2793

2795

else

if (x <= 1.) then

2794

2796

if (x <= 1.) then

d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))

2795

2797

d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))

r_res=datan(dsqrt(d)/x)

2796

2798

r_res=datan(dsqrt(d)/x)

else

2797

2799

else

apx=a+x

2798

2800

apx=a+x

d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))

2799

2801

d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))

r_res=datan((y*dsqrt(d))/x);

2800

2802

r_res=datan((y*dsqrt(d))/x);

end if

2801

2803

end if

2802

2804

end if

2803

2805

if (a <= a_crossover) then

2804

2806

if (a <= a_crossover) then

if (x < 1.) then

2805

2807

if (x < 1.) then

am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))

2806

2808

am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))

else

2807

2809

else

am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))

2808

2810

am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))

end if

2809

2811

end if

i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));

2810

2812

i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));

else

2811

2813

else

i_res = dlog (a + dsqrt (a*a - 1.));

2812

2814

i_res = dlog (a + dsqrt (a*a - 1.));

end if

2813

2815

end if

if (rz <0.) then

2814

2816

if (rz <0.) then

r_res=fvn_pi-r_res

2815

2817

r_res=fvn_pi-r_res

end if

2816

2818

end if

i_res=-sign(1._8,iz)*i_res

2817

2819

i_res=-sign(1._8,iz)*i_res

fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res)

2818

2820

fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res)

2819

2821

end function fvn_z_acos

2820

2822

end function fvn_z_acos

2821

2823

function fvn_z_acos_real(r)

2822

2824

function fvn_z_acos_real(r)

! return the double complex arc cosinus for a

2823

2825

! return the double complex arc cosinus for a

! double precision argument

2824

2826

! double precision argument

implicit none

2825

2827

implicit none

real(kind=8) :: r

2826

2828

real(kind=8) :: r

complex(kind=8) :: fvn_z_acos_real

2827

2829

complex(kind=8) :: fvn_z_acos_real

2828

2830

if (dabs(r)<=1._8) then

2829

2831

if (dabs(r)<=1._8) then

fvn_z_acos_real=dcmplx(dacos(r))

2830

2832

fvn_z_acos_real=dcmplx(dacos(r))

return

2831

2833

return

end if

2832

2834

end if

if (r < 0._8) then

2833

2835

if (r < 0._8) then

fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r))

2834

2836

fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r))

else

2835

2837

else

fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r))

2836

2838

fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r))

end if

2837

2839

end if

end function

2838

2840

end function

2839

2841

2840

2842

function fvn_z_asin(z)

2841

2843

function fvn_z_asin(z)

! double complex arcsin function derived from

2842

2844

! double complex arcsin function derived from

! the c gsl library

2843

2845

! the c gsl library

! http://www.gnu.org/software/gsl/

2844

2846

! http://www.gnu.org/software/gsl/

implicit none

2845

2847

implicit none

complex(kind=8) :: fvn_z_asin

2846

2848

complex(kind=8) :: fvn_z_asin

complex(kind=8) :: z

2847

2849

complex(kind=8) :: z

real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1

2848

2850

real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1

real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8

2849

2851

real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8

real(kind=8) :: r_res,i_res

2850

2852

real(kind=8) :: r_res,i_res

2851

2853

rz=dreal(z)

2852

2854

rz=dreal(z)

iz=dimag(z)

2853

2855

iz=dimag(z)

if ( iz == 0._8 ) then

2854

2856

if ( iz == 0._8 ) then

! z is real

2855

2857

! z is real

fvn_z_asin=fvn_z_asin_real(rz)

2856

2858

fvn_z_asin=fvn_z_asin_real(rz)

return

2857

2859

return

end if

2858

2860

end if

2859

2861

x=dabs(rz)

2860

2862

x=dabs(rz)

y=dabs(iz)

2861

2863

y=dabs(iz)

r=fvn_d_hypot(x+1.,y)

2862

2864

r=fvn_d_hypot(x+1.,y)

s=fvn_d_hypot(x-1.,y)

2863

2865

s=fvn_d_hypot(x-1.,y)

a=0.5*(r + s)

2864

2866

a=0.5*(r + s)

b=x/a

2865

2867

b=x/a

y2=y*y

2866

2868

y2=y*y

2867

2869

if (b <= b_crossover) then

2868

2870

if (b <= b_crossover) then

r_res=dasin(b)

2869

2871

r_res=dasin(b)

else

2870

2872

else

if (x <= 1.) then

2871

2873

if (x <= 1.) then

d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))

2872

2874

d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))

r_res=datan(x/dsqrt(d))

2873

2875

r_res=datan(x/dsqrt(d))

else

2874

2876

else

apx=a+x

2875

2877

apx=a+x

d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))

2876

2878

d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))

r_res=datan(x/(y*dsqrt(d)));

2877

2879

r_res=datan(x/(y*dsqrt(d)));

end if

2878

2880

end if

2879

2881

end if

2880

2882

if (a <= a_crossover) then

2881

2883

if (a <= a_crossover) then

if (x < 1.) then

2882

2884

if (x < 1.) then

am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))

2883

2885

am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))

else

2884

2886

else

am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))

2885

2887

am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))

end if

2886

2888

end if

i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));

2887

2889

i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));

else

2888

2890

else

i_res = dlog (a + dsqrt (a*a - 1.));

2889

2891

i_res = dlog (a + dsqrt (a*a - 1.));

end if

2890

2892

end if

r_res=sign(1._8,rz)*r_res

2891

2893

r_res=sign(1._8,rz)*r_res

i_res=sign(1._8,iz)*i_res

2892

2894

i_res=sign(1._8,iz)*i_res

fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res)

2893

2895

fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res)

2894

2896

end function fvn_z_asin

2895

2897

end function fvn_z_asin

2896

2898

function fvn_z_asin_real(r)

2897

2899

function fvn_z_asin_real(r)

! return the double complex arc sinus for a

2898

2900

! return the double complex arc sinus for a

! double precision argument

2899

2901

! double precision argument

implicit none

2900

2902

implicit none

real(kind=8) :: r

2901

2903

real(kind=8) :: r

complex(kind=8) :: fvn_z_asin_real

2902

2904

complex(kind=8) :: fvn_z_asin_real

2903

2905

if (dabs(r)<=1._8) then

2904

2906

if (dabs(r)<=1._8) then

fvn_z_asin_real=dcmplx(dasin(r))

2905

2907

fvn_z_asin_real=dcmplx(dasin(r))

return

2906

2908

return

end if

2907

2909

end if

if (r < 0._8) then

2908

2910

if (r < 0._8) then

fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r))

2909

2911

fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r))

else

2910

2912

else

fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r))

2911

2913

fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r))

end if

2912

2914

end if

end function fvn_z_asin_real

2913

2915

end function fvn_z_asin_real

2914

2916

function fvn_d_acosh(r)

2915

2917

function fvn_d_acosh(r)

! return the arc hyperbolic cosine

2916

2918

! return the arc hyperbolic cosine

implicit none

2917

2919

implicit none

real(kind=8) :: r

2918

2920

real(kind=8) :: r

GITLAB

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