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

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

25

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

26

! a function which return identity matrix,

27

! evaluation of nterm chebyshev series

! September 2007 : added sparse system solving by interfacing umfpack

25

28

! September 2007 : added sparse system solving by interfacing umfpack

! June 2007 : added some complex trigonometric functions

26

29

! June 2007 : added some complex trigonometric functions

!

27

30

!

! TO DO LIST :

28

31

! TO DO LIST :

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

29

32

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

! eigenvalues are given with no particular order.

30

33

! eigenvalues are given with no particular order.

! + Generic interface for fvn_x_name family -> fvn_name

31

34

! + Generic interface for fvn_x_name family -> fvn_name

! + Make some parameters optional, status for example

32

35

! + Make some parameters optional, status for example

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

33

36

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

! + unify quadpack routines

34

37

! + unify quadpack routines

! + ...

35

38

! + ...

!

36

39

!

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

37

40

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

38

41

implicit none

39

42

implicit none

! We define pi and i for the module

40

43

! We define pi and i for the module

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

41

44

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

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

42

45

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

43

46

44

47

! All quadpack routines are private to the module

45

48

! All quadpack routines are private to the module

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

46

49

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, &

47

50

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, &

48

51

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

49

52

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

50

53

51

54

contains

52

55

contains

53

56

57

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

58

!

59

! Identity Matrix

60

!

61

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

62

function fvn_d_ident(n)

63

implicit none

64

integer(kind=4) :: n

65

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

66

67

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

68

integer(kind=4) :: i

69

70

vect=0._8

71

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

72

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

73

end function

74

75

function fvn_s_ident(n)

76

implicit none

77

integer(kind=4) :: n

78

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

79

80

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

81

integer(kind=4) :: i

82

83

vect=0._4

84

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

85

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

86

end function

87

88

function fvn_c_ident(n)

89

implicit none

90

integer(kind=4) :: n

91

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

92

93

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

94

integer(kind=4) :: i

95

96

vect=(0._4,0._4)

97

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

98

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

99

end function

100

101

function fvn_z_ident(n)

102

implicit none

103

integer(kind=4) :: n

104

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

105

106

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

107

integer(kind=4) :: i

108

109

vect=(0._8,0._8)

110

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

111

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

112

end function

113

114

115

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

54

116

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

!

55

117

!

! Matrix inversion subroutines

56

118

! Matrix inversion subroutines

!

57

119

!

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

58

120

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

subroutine fvn_s_matinv(d,a,inva,status)

59

121

subroutine fvn_s_matinv(d,a,inva,status)

!

60

122

!

! Matrix inversion of a real matrix using BLAS and LAPACK

61

123

! Matrix inversion of a real matrix using BLAS and LAPACK

!

62

124

!

! d (in) : matrix rank

63

125

! d (in) : matrix rank

! a (in) : input matrix

64

126

! a (in) : input matrix

! inva (out) : inversed matrix

65

127

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

66

128

! status (ou) : =0 if something failed

!

67

129

!

implicit none

68

130

implicit none

integer, intent(in) :: d

69

131

integer, intent(in) :: d

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

70

132

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

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

71

133

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

integer, intent(out) :: status

72

134

integer, intent(out) :: status

73

135

integer, allocatable :: ipiv(:)

74

136

integer, allocatable :: ipiv(:)

real, allocatable :: work(:)

75

137

real, allocatable :: work(:)

real twork(1)

76

138

real twork(1)

integer :: info

77

139

integer :: info

integer :: lwork

78

140

integer :: lwork

79

141

status=1

80

142

status=1

81

143

allocate(ipiv(d))

82

144

allocate(ipiv(d))

! copy a into inva using BLAS

83

145

! copy a into inva using BLAS

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

84

146

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

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

85

147

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

! LU factorization using LAPACK

86

148

! LU factorization using LAPACK

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

87

149

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

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

88

150

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

if (info /= 0) then

89

151

if (info /= 0) then

status=0

90

152

status=0

deallocate(ipiv)

91

153

deallocate(ipiv)

return

92

154

return

end if

93

155

end if

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

94

156

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

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

95

157

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

lwork=int(twork(1))

96

158

lwork=int(twork(1))

allocate(work(lwork))

97

159

allocate(work(lwork))

! Matrix inversion using LAPACK

98

160

! Matrix inversion using LAPACK

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

99

161

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

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

100

162

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

if (info /= 0) then

101

163

if (info /= 0) then

status=0

102

164

status=0

end if

103

165

end if

deallocate(work)

104

166

deallocate(work)

deallocate(ipiv)

105

167

deallocate(ipiv)

end subroutine

106

168

end subroutine

107

169

subroutine fvn_d_matinv(d,a,inva,status)

108

170

subroutine fvn_d_matinv(d,a,inva,status)

!

109

171

!

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

110

172

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

!

111

173

!

! d (in) : matrix rank

112

174

! d (in) : matrix rank

! a (in) : input matrix

113

175

! a (in) : input matrix

! inva (out) : inversed matrix

114

176

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

115

177

! status (ou) : =0 if something failed

!

116

178

!

implicit none

117

179

implicit none

integer, intent(in) :: d

118

180

integer, intent(in) :: d

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

119

181

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

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

120

182

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

integer, intent(out) :: status

121

183

integer, intent(out) :: status

122

184

integer, allocatable :: ipiv(:)

123

185

integer, allocatable :: ipiv(:)

double precision, allocatable :: work(:)

124

186

double precision, allocatable :: work(:)

double precision :: twork(1)

125

187

double precision :: twork(1)

integer :: info

126

188

integer :: info

integer :: lwork

127

189

integer :: lwork

128

190

status=1

129

191

status=1

130

192

allocate(ipiv(d))

131

193

allocate(ipiv(d))

! copy a into inva using BLAS

132

194

! copy a into inva using BLAS

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

133

195

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

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

134

196

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

! LU factorization using LAPACK

135

197

! LU factorization using LAPACK

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

136

198

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

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

137

199

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

if (info /= 0) then

138

200

if (info /= 0) then

status=0

139

201

status=0

deallocate(ipiv)

140

202

deallocate(ipiv)

return

141

203

return

end if

142

204

end if

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

143

205

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

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

144

206

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

lwork=int(twork(1))

145

207

lwork=int(twork(1))

allocate(work(lwork))

146

208

allocate(work(lwork))

! Matrix inversion using LAPACK

147

209

! Matrix inversion using LAPACK

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

148

210

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

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

149

211

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

if (info /= 0) then

150

212

if (info /= 0) then

status=0

151

213

status=0

end if

152

214

end if

deallocate(work)

153

215

deallocate(work)

deallocate(ipiv)

154

216

deallocate(ipiv)

end subroutine

155

217

end subroutine

156

218

subroutine fvn_c_matinv(d,a,inva,status)

157

219

subroutine fvn_c_matinv(d,a,inva,status)

!

158

220

!

! Matrix inversion of a complex matrix using BLAS and LAPACK

159

221

! Matrix inversion of a complex matrix using BLAS and LAPACK

!

160

222

!

! d (in) : matrix rank

161

223

! d (in) : matrix rank

! a (in) : input matrix

162

224

! a (in) : input matrix

! inva (out) : inversed matrix

163

225

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

164

226

! status (ou) : =0 if something failed

!

165

227

!

implicit none

166

228

implicit none

integer, intent(in) :: d

167

229

integer, intent(in) :: d

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

168

230

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

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

169

231

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

integer, intent(out) :: status

170

232

integer, intent(out) :: status

171

233

integer, allocatable :: ipiv(:)

172

234

integer, allocatable :: ipiv(:)

complex, allocatable :: work(:)

173

235

complex, allocatable :: work(:)

complex :: twork(1)

174

236

complex :: twork(1)

integer :: info

175

237

integer :: info

integer :: lwork

176

238

integer :: lwork

177

239

status=1

178

240

status=1

179

241

allocate(ipiv(d))

180

242

allocate(ipiv(d))

! copy a into inva using BLAS

181

243

! copy a into inva using BLAS

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

182

244

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

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

183

245

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

184

246

! LU factorization using LAPACK

185

247

! LU factorization using LAPACK

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

186

248

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

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

187

249

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

if (info /= 0) then

188

250

if (info /= 0) then

status=0

189

251

status=0

deallocate(ipiv)

190

252

deallocate(ipiv)

return

191

253

return

end if

192

254

end if

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

193

255

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

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

194

256

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

lwork=int(twork(1))

195

257

lwork=int(twork(1))

allocate(work(lwork))

196

258

allocate(work(lwork))

! Matrix inversion using LAPACK

197

259

! Matrix inversion using LAPACK

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

198

260

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

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

199

261

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

if (info /= 0) then

200

262

if (info /= 0) then

status=0

201

263

status=0

end if

202

264

end if

deallocate(work)

203

265

deallocate(work)

deallocate(ipiv)

204

266

deallocate(ipiv)

end subroutine

205

267

end subroutine

206

268

subroutine fvn_z_matinv(d,a,inva,status)

207

269

subroutine fvn_z_matinv(d,a,inva,status)

!

208

270

!

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

209

271

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

!

210

272

!

! d (in) : matrix rank

211

273

! d (in) : matrix rank

! a (in) : input matrix

212

274

! a (in) : input matrix

! inva (out) : inversed matrix

213

275

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

214

276

! status (ou) : =0 if something failed

!

215

277

!

implicit none

216

278

implicit none

integer, intent(in) :: d

217

279

integer, intent(in) :: d

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

218

280

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

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

219

281

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

integer, intent(out) :: status

220

282

integer, intent(out) :: status

221

283

integer, allocatable :: ipiv(:)

222

284

integer, allocatable :: ipiv(:)

double complex, allocatable :: work(:)

223

285

double complex, allocatable :: work(:)

double complex :: twork(1)

224

286

double complex :: twork(1)

integer :: info

225

287

integer :: info

integer :: lwork

226

288

integer :: lwork

227

289

status=1

228

290

status=1

229

291

allocate(ipiv(d))

230

292

allocate(ipiv(d))

! copy a into inva using BLAS

231

293

! copy a into inva using BLAS

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

232

294

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

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

233

295

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

234

296

! LU factorization using LAPACK

235

297

! LU factorization using LAPACK

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

236

298

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

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

237

299

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

if (info /= 0) then

238

300

if (info /= 0) then

status=0

239

301

status=0

deallocate(ipiv)

240

302

deallocate(ipiv)

return

241

303

return

end if

242

304

end if

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

243

305

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

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

244

306

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

lwork=int(twork(1))

245

307

lwork=int(twork(1))

allocate(work(lwork))

246

308

allocate(work(lwork))

! Matrix inversion using LAPACK

247

309

! Matrix inversion using LAPACK

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

248

310

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

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

249

311

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

if (info /= 0) then

250

312

if (info /= 0) then

status=0

251

313

status=0

end if

252

314

end if

deallocate(work)

253

315

deallocate(work)

deallocate(ipiv)

254

316

deallocate(ipiv)

end subroutine

255

317

end subroutine

256

318

257

319

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

258

320

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

!

259

321

!

! Determinants

260

322

! Determinants

!

261

323

!

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

262

324

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

function fvn_s_det(d,a,status)

263

325

function fvn_s_det(d,a,status)

!

264

326

!

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

265

327

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

!

266

328

!

! d (in) : matrix rank

267

329

! d (in) : matrix rank

! a (in) : The Matrix

268

330

! a (in) : The Matrix

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

269

331

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

!

270

332

!

implicit none

271

333

implicit none

integer, intent(in) :: d

272

334

integer, intent(in) :: d

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

273

335

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

integer, intent(out) :: status

274

336

integer, intent(out) :: status

real :: fvn_s_det

275

337

real :: fvn_s_det

276

338

real, allocatable :: wc_a(:,:)

277

339

real, allocatable :: wc_a(:,:)

integer, allocatable :: ipiv(:)

278

340

integer, allocatable :: ipiv(:)

integer :: info,i

279

341

integer :: info,i

280

342

status=1

281

343

status=1

allocate(wc_a(d,d))

282

344

allocate(wc_a(d,d))

allocate(ipiv(d))

283

345

allocate(ipiv(d))

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

284

346

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

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

285

347

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

if (info/= 0) then

286

348

if (info/= 0) then

status=0

287

349

status=0

fvn_s_det=0.e0

288

350

fvn_s_det=0.e0

deallocate(ipiv)

289

351

deallocate(ipiv)

deallocate(wc_a)

290

352

deallocate(wc_a)

return

291

353

return

end if

292

354

end if

fvn_s_det=1.e0

293

355

fvn_s_det=1.e0

do i=1,d

294

356

do i=1,d

if (ipiv(i)==i) then

295

357

if (ipiv(i)==i) then

fvn_s_det=fvn_s_det*wc_a(i,i)

296

358

fvn_s_det=fvn_s_det*wc_a(i,i)

else

297

359

else

fvn_s_det=-fvn_s_det*wc_a(i,i)

298

360

fvn_s_det=-fvn_s_det*wc_a(i,i)

end if

299

361

end if

end do

300

362

end do

deallocate(ipiv)

301

363

deallocate(ipiv)

deallocate(wc_a)

302

364

deallocate(wc_a)

303

365

end function

304

366

end function

305

367

function fvn_d_det(d,a,status)

306

368

function fvn_d_det(d,a,status)

!

307

369

!

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

308

370

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

!

309

371

!

! d (in) : matrix rank

310

372

! d (in) : matrix rank

! a (in) : The Matrix

311

373

! a (in) : The Matrix

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

312

374

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

!

313

375

!

implicit none

314

376

implicit none

integer, intent(in) :: d

315

377

integer, intent(in) :: d

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

316

378

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

integer, intent(out) :: status

317

379

integer, intent(out) :: status

double precision :: fvn_d_det

318

380

double precision :: fvn_d_det

319

381

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

320

382

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

integer, allocatable :: ipiv(:)

321

383

integer, allocatable :: ipiv(:)

integer :: info,i

322

384

integer :: info,i

323

385

status=1

324

386

status=1

allocate(wc_a(d,d))

325

387

allocate(wc_a(d,d))

allocate(ipiv(d))

326

388

allocate(ipiv(d))

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

327

389

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

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

328

390

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

if (info/= 0) then

329

391

if (info/= 0) then

status=0

330

392

status=0

fvn_d_det=0.d0

331

393

fvn_d_det=0.d0

deallocate(ipiv)

332

394

deallocate(ipiv)

deallocate(wc_a)

333

395

deallocate(wc_a)

return

334

396

return

end if

335

397

end if

fvn_d_det=1.d0

336

398

fvn_d_det=1.d0

do i=1,d

337

399

do i=1,d

if (ipiv(i)==i) then

338

400

if (ipiv(i)==i) then

fvn_d_det=fvn_d_det*wc_a(i,i)

339

401

fvn_d_det=fvn_d_det*wc_a(i,i)

else

340

402

else

fvn_d_det=-fvn_d_det*wc_a(i,i)

341

403

fvn_d_det=-fvn_d_det*wc_a(i,i)

end if

342

404

end if

end do

343

405

end do

deallocate(ipiv)

344

406

deallocate(ipiv)

deallocate(wc_a)

345

407

deallocate(wc_a)

346

408

end function

347

409

end function

348

410

function fvn_c_det(d,a,status) !

349

411

function fvn_c_det(d,a,status) !

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

350

412

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

!

351

413

!

! d (in) : matrix rank

352

414

! d (in) : matrix rank

! a (in) : The Matrix

353

415

! a (in) : The Matrix

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

354

416

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

!

355

417

!

implicit none

356

418

implicit none

integer, intent(in) :: d

357

419

integer, intent(in) :: d

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

358

420

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

integer, intent(out) :: status

359

421

integer, intent(out) :: status

complex :: fvn_c_det

360

422

complex :: fvn_c_det

361

423

complex, allocatable :: wc_a(:,:)

362

424

complex, allocatable :: wc_a(:,:)

integer, allocatable :: ipiv(:)

363

425

integer, allocatable :: ipiv(:)

integer :: info,i

364

426

integer :: info,i

365

427

status=1

366

428

status=1

allocate(wc_a(d,d))

367

429

allocate(wc_a(d,d))

allocate(ipiv(d))

368

430

allocate(ipiv(d))

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

369

431

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

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

370

432

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

if (info/= 0) then

371

433

if (info/= 0) then

status=0

372

434

status=0

fvn_c_det=(0.e0,0.e0)

373

435

fvn_c_det=(0.e0,0.e0)

deallocate(ipiv)

374

436

deallocate(ipiv)

deallocate(wc_a)

375

437

deallocate(wc_a)

return

376

438

return

end if

377

439

end if

fvn_c_det=(1.e0,0.e0)

378

440

fvn_c_det=(1.e0,0.e0)

do i=1,d

379

441

do i=1,d

if (ipiv(i)==i) then

380

442

if (ipiv(i)==i) then

fvn_c_det=fvn_c_det*wc_a(i,i)

381

443

fvn_c_det=fvn_c_det*wc_a(i,i)

else

382

444

else

fvn_c_det=-fvn_c_det*wc_a(i,i)

383

445

fvn_c_det=-fvn_c_det*wc_a(i,i)

end if

384

446

end if

end do

385

447

end do

deallocate(ipiv)

386

448

deallocate(ipiv)

deallocate(wc_a)

387

449

deallocate(wc_a)

388

450

end function

389

451

end function

390

452

function fvn_z_det(d,a,status)

391

453

function fvn_z_det(d,a,status)

!

392

454

!

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

393

455

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

!

394

456

!

! d (in) : matrix rank

395

457

! d (in) : matrix rank

! a (in) : The Matrix

396

458

! a (in) : The Matrix

! det (out) : determinant

397

459

! det (out) : determinant

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

398

460

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

!

399

461

!

implicit none

400

462

implicit none

integer, intent(in) :: d

401

463

integer, intent(in) :: d

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

402

464

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

integer, intent(out) :: status

403

465

integer, intent(out) :: status

double complex :: fvn_z_det

404

466

double complex :: fvn_z_det

405

467

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

406

468

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

integer, allocatable :: ipiv(:)

407

469

integer, allocatable :: ipiv(:)

integer :: info,i

408

470

integer :: info,i

409

471

status=1

410

472

status=1

allocate(wc_a(d,d))

411

473

allocate(wc_a(d,d))

allocate(ipiv(d))

412

474

allocate(ipiv(d))

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

413

475

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

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

414

476

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

if (info/= 0) then

415

477

if (info/= 0) then

status=0

416

478

status=0

fvn_z_det=(0.d0,0.d0)

417

479

fvn_z_det=(0.d0,0.d0)

deallocate(ipiv)

418

480

deallocate(ipiv)

deallocate(wc_a)

419

481

deallocate(wc_a)

return

420

482

return

end if

421

483

end if

fvn_z_det=(1.d0,0.d0)

422

484

fvn_z_det=(1.d0,0.d0)

do i=1,d

423

485

do i=1,d

if (ipiv(i)==i) then

424

486

if (ipiv(i)==i) then

fvn_z_det=fvn_z_det*wc_a(i,i)

425

487

fvn_z_det=fvn_z_det*wc_a(i,i)

else

426

488

else

fvn_z_det=-fvn_z_det*wc_a(i,i)

427

489

fvn_z_det=-fvn_z_det*wc_a(i,i)

end if

428

490

end if

end do

429

491

end do

deallocate(ipiv)

430

492

deallocate(ipiv)

deallocate(wc_a)

431

493

deallocate(wc_a)

432

494

end function

433

495

end function

434

496

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

435

497

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

!

436

498

!

! Condition test

437

499

! Condition test

!

438

500

!

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

439

501

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

! 1-norm

440

502

! 1-norm

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

441

503

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

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

442

504

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

!

443

505

!

subroutine fvn_s_matcon(d,a,rcond,status)

444

506

subroutine fvn_s_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

445

507

! Matrix condition (reciprocal of condition number)

!

446

508

!

! d (in) : matrix rank

447

509

! d (in) : matrix rank

! a (in) : The Matrix

448

510

! a (in) : The Matrix

! rcond (out) : guess what

449

511

! rcond (out) : guess what

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

450

512

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

!

451

513

!

implicit none

452

514

implicit none

integer, intent(in) :: d

453

515

integer, intent(in) :: d

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

454

516

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

real, intent(out) :: rcond

455

517

real, intent(out) :: rcond

integer, intent(out) :: status

456

518

integer, intent(out) :: status

457

519

real, allocatable :: work(:)

458

520

real, allocatable :: work(:)

integer, allocatable :: iwork(:)

459

521

integer, allocatable :: iwork(:)

real :: anorm

460

522

real :: anorm

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

461

523

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

integer :: info

462

524

integer :: info

integer, allocatable :: ipiv(:)

463

525

integer, allocatable :: ipiv(:)

464

526

real, external :: slange

465

527

real, external :: slange

466

528

467

529

status=1

468

530

status=1

469

531

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

470

532

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

471

533

allocate(wc_a(d,d))

472

534

allocate(wc_a(d,d))

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

473

535

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

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

474

536

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

475

537

allocate(ipiv(d))

476

538

allocate(ipiv(d))

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

477

539

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

if (info /= 0) then

478

540

if (info /= 0) then

status=0

479

541

status=0

deallocate(ipiv)

480

542

deallocate(ipiv)

deallocate(wc_a)

481

543

deallocate(wc_a)

return

482

544

return

end if

483

545

end if

allocate(work(4*d))

484

546

allocate(work(4*d))

allocate(iwork(d))

485

547

allocate(iwork(d))

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

486

548

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

if (info /= 0) then

487

549

if (info /= 0) then

status=0

488

550

status=0

end if

489

551

end if

deallocate(iwork)

490

552

deallocate(iwork)

deallocate(work)

491

553

deallocate(work)

deallocate(ipiv)

492

554

deallocate(ipiv)

deallocate(wc_a)

493

555

deallocate(wc_a)

494

556

end subroutine

495

557

end subroutine

496

558

subroutine fvn_d_matcon(d,a,rcond,status)

497

559

subroutine fvn_d_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

498

560

! Matrix condition (reciprocal of condition number)

!

499

561

!

! d (in) : matrix rank

500

562

! d (in) : matrix rank

! a (in) : The Matrix

501

563

! a (in) : The Matrix

! rcond (out) : guess what

502

564

! rcond (out) : guess what

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

503

565

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

!

504

566

!

implicit none

505

567

implicit none

integer, intent(in) :: d

506

568

integer, intent(in) :: d

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

507

569

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

double precision, intent(out) :: rcond

508

570

double precision, intent(out) :: rcond

integer, intent(out) :: status

509

571

integer, intent(out) :: status

510

572

double precision, allocatable :: work(:)

511

573

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

512

574

integer, allocatable :: iwork(:)

double precision :: anorm

513

575

double precision :: anorm

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

514

576

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

integer :: info

515

577

integer :: info

integer, allocatable :: ipiv(:)

516

578

integer, allocatable :: ipiv(:)

517

579

double precision, external :: dlange

518

580

double precision, external :: dlange

519

581

520

582

status=1

521

583

status=1

522

584

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

523

585

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

524

586

allocate(wc_a(d,d))

525

587

allocate(wc_a(d,d))

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

526

588

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

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

527

589

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

528

590

allocate(ipiv(d))

529

591

allocate(ipiv(d))

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

530

592

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

if (info /= 0) then

531

593

if (info /= 0) then

status=0

532

594

status=0

deallocate(ipiv)

533

595

deallocate(ipiv)

deallocate(wc_a)

534

596

deallocate(wc_a)

return

535

597

return

end if

536

598

end if

537

599

allocate(work(4*d))

538

600

allocate(work(4*d))

allocate(iwork(d))

539

601

allocate(iwork(d))

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

540

602

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

if (info /= 0) then

541

603

if (info /= 0) then

status=0

542

604

status=0

end if

543

605

end if

deallocate(iwork)

544

606

deallocate(iwork)

deallocate(work)

545

607

deallocate(work)

deallocate(ipiv)

546

608

deallocate(ipiv)

deallocate(wc_a)

547

609

deallocate(wc_a)

548

610

end subroutine

549

611

end subroutine

550

612

subroutine fvn_c_matcon(d,a,rcond,status)

551

613

subroutine fvn_c_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

552

614

! Matrix condition (reciprocal of condition number)

!

553

615

!

! d (in) : matrix rank

554

616

! d (in) : matrix rank

! a (in) : The Matrix

555

617

! a (in) : The Matrix

! rcond (out) : guess what

556

618

! rcond (out) : guess what

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

557

619

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

!

558

620

!

implicit none

559

621

implicit none

integer, intent(in) :: d

560

622

integer, intent(in) :: d

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

561

623

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

real, intent(out) :: rcond

562

624

real, intent(out) :: rcond

integer, intent(out) :: status

563

625

integer, intent(out) :: status

564

626

real, allocatable :: rwork(:)

565

627

real, allocatable :: rwork(:)

complex, allocatable :: work(:)

566

628

complex, allocatable :: work(:)

integer, allocatable :: iwork(:)

567

629

integer, allocatable :: iwork(:)

real :: anorm

568

630

real :: anorm

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

569

631

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

integer :: info

570

632

integer :: info

integer, allocatable :: ipiv(:)

571

633

integer, allocatable :: ipiv(:)

572

634

real, external :: clange

573

635

real, external :: clange

574

636

575

637

status=1

576

638

status=1

577

639

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

578

640

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

579

641

allocate(wc_a(d,d))

580

642

allocate(wc_a(d,d))

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

581

643

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

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

582

644

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

583

645

allocate(ipiv(d))

584

646

allocate(ipiv(d))

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

585

647

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

if (info /= 0) then

586

648

if (info /= 0) then

status=0

587

649

status=0

deallocate(ipiv)

588

650

deallocate(ipiv)

deallocate(wc_a)

589

651

deallocate(wc_a)

return

590

652

return

end if

591

653

end if

allocate(work(2*d))

592

654

allocate(work(2*d))

allocate(rwork(2*d))

593

655

allocate(rwork(2*d))

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

594

656

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

if (info /= 0) then

595

657

if (info /= 0) then

status=0

596

658

status=0

end if

597

659

end if

deallocate(rwork)

598

660

deallocate(rwork)

deallocate(work)

599

661

deallocate(work)

deallocate(ipiv)

600

662

deallocate(ipiv)

deallocate(wc_a)

601

663

deallocate(wc_a)

end subroutine

602

664

end subroutine

603

665

subroutine fvn_z_matcon(d,a,rcond,status)

604

666

subroutine fvn_z_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

605

667

! Matrix condition (reciprocal of condition number)

!

606

668

!

! d (in) : matrix rank

607

669

! d (in) : matrix rank

! a (in) : The Matrix

608

670

! a (in) : The Matrix

! rcond (out) : guess what

609

671

! rcond (out) : guess what

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

610

672

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

!

611

673

!

implicit none

612

674

implicit none

integer, intent(in) :: d

613

675

integer, intent(in) :: d

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

614

676

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

double precision, intent(out) :: rcond

615

677

double precision, intent(out) :: rcond

integer, intent(out) :: status

616

678

integer, intent(out) :: status

617

679

double complex, allocatable :: work(:)

618

680

double complex, allocatable :: work(:)

double precision, allocatable :: rwork(:)

619

681

double precision, allocatable :: rwork(:)

double precision :: anorm

620

682

double precision :: anorm

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

621

683

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

integer :: info

622

684

integer :: info

integer, allocatable :: ipiv(:)

623

685

integer, allocatable :: ipiv(:)

624

686

double precision, external :: zlange

625

687

double precision, external :: zlange

626

688

627

689

status=1

628

690

status=1

629

691

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

630

692

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

631

693

allocate(wc_a(d,d))

632

694

allocate(wc_a(d,d))

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

633

695

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

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

634

696

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

635

697

allocate(ipiv(d))

636

698

allocate(ipiv(d))

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

637

699

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

if (info /= 0) then

638

700

if (info /= 0) then

status=0

639

701

status=0

deallocate(ipiv)

640

702

deallocate(ipiv)

deallocate(wc_a)

641

703

deallocate(wc_a)

return

642

704

return

end if

643

705

end if

644

706

allocate(work(2*d))

645

707

allocate(work(2*d))

allocate(rwork(2*d))

646

708

allocate(rwork(2*d))

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

647

709

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

if (info /= 0) then

648

710

if (info /= 0) then

status=0

649

711

status=0

end if

650

712

end if

deallocate(rwork)

651

713

deallocate(rwork)

deallocate(work)

652

714

deallocate(work)

deallocate(ipiv)

653

715

deallocate(ipiv)

deallocate(wc_a)

654

716

deallocate(wc_a)

end subroutine

655

717

end subroutine

656

718

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

657

719

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

!

658

720

!

! Valeurs propres/ Vecteurs propre

659

721

! Valeurs propres/ Vecteurs propre

!

660

722

!

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

661

723

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

662

724

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

663

725

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

!

664

726

!

! integer d (in) : matrice rank

665

727

! integer d (in) : matrice rank

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

666

728

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

! complex evala(d) (out) : eigenvalues

667

729

! complex evala(d) (out) : eigenvalues

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

668

730

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

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

669

731

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

!

670

732

!

! interfacing Lapack routine SGEEV

671

733

! interfacing Lapack routine SGEEV

implicit none

672

734

implicit none

integer, intent(in) :: d

673

735

integer, intent(in) :: d

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

674

736

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

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

675

737

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

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

676

738

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

integer, intent(out) :: status

677

739

integer, intent(out) :: status

678

740

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

679

741

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

integer :: info

680

742

integer :: info

integer :: lwork

681

743

integer :: lwork

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

682

744

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

real :: vl ! unused but necessary for the call

683

745

real :: vl ! unused but necessary for the call

real, allocatable :: vr(:,:)

684

746

real, allocatable :: vr(:,:)

real, allocatable :: work(:)

685

747

real, allocatable :: work(:)

real :: twork(1)

686

748

real :: twork(1)

integer i

687

749

integer i

integer j

688

750

integer j

689

751

! making a working copy of a

690

752

! making a working copy of a

allocate(wc_a(d,d))

691

753

allocate(wc_a(d,d))

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

692

754

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

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

693

755

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

694

756

allocate(wr(d))

695

757

allocate(wr(d))

allocate(wi(d))

696

758

allocate(wi(d))

allocate(vr(d,d))

697

759

allocate(vr(d,d))

! query optimal work size

698

760

! query optimal work size

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

699

761

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

lwork=int(twork(1))

700

762

lwork=int(twork(1))

allocate(work(lwork))

701

763

allocate(work(lwork))

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

702

764

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

703

765

if (info /= 0) then

704

766

if (info /= 0) then

status=0

705

767

status=0

deallocate(work)

706

768

deallocate(work)

deallocate(vr)

707

769

deallocate(vr)

deallocate(wi)

708

770

deallocate(wi)

deallocate(wr)

709

771

deallocate(wr)

deallocate(wc_a)

710

772

deallocate(wc_a)

return

711

773

return

end if

712

774

end if

713

775

! now fill in the results

714

776

! now fill in the results

i=1

715

777

i=1

do while(i<=d)

716

778

do while(i<=d)

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

717

779

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

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

718

780

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

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

719

781

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

else ! eigenvalue is complex

720

782

else ! eigenvalue is complex

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

721

783

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

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

722

784

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

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

723

785

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

i=i+1

724

786

i=i+1

end if

725

787

end if

i=i+1

726

788

i=i+1

enddo

727

789

enddo

deallocate(work)

728

790

deallocate(work)

deallocate(vr)

729

791

deallocate(vr)

deallocate(wi)

730

792

deallocate(wi)

deallocate(wr)

731

793

deallocate(wr)

deallocate(wc_a)

732

794

deallocate(wc_a)

733

795

end subroutine

734

796

end subroutine

735

797

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

736

798

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

!

737

799

!

! integer d (in) : matrice rank

738

800

! integer d (in) : matrice rank

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

739

801

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

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

740

802

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

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

741

803

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

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

742

804

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

!

743

805

!

! interfacing Lapack routine DGEEV

744

806

! interfacing Lapack routine DGEEV

implicit none

745

807

implicit none

integer, intent(in) :: d

746

808

integer, intent(in) :: d

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

747

809

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

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

748

810

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

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

749

811

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

integer, intent(out) :: status

750

812

integer, intent(out) :: status

751

813

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

752

814

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

integer :: info

753

815

integer :: info

integer :: lwork

754

816

integer :: lwork

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

755

817

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

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

756

818

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

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

757

819

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

double precision, allocatable :: work(:)

758

820

double precision, allocatable :: work(:)

double precision :: twork(1)

759

821

double precision :: twork(1)

integer i

760

822

integer i

integer j

761

823

integer j

762

824

! making a working copy of a

763

825

! making a working copy of a

allocate(wc_a(d,d))

764

826

allocate(wc_a(d,d))

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

765

827

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

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

766

828

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

767

829

allocate(wr(d))

768

830

allocate(wr(d))

allocate(wi(d))

769

831

allocate(wi(d))

allocate(vr(d,d))

770

832

allocate(vr(d,d))

! query optimal work size

771

833

! query optimal work size

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

772

834

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

lwork=int(twork(1))

773

835

lwork=int(twork(1))

allocate(work(lwork))

774

836

allocate(work(lwork))

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

775

837

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

776

838

if (info /= 0) then

777

839

if (info /= 0) then

status=0

778

840

status=0

deallocate(work)

779

841

deallocate(work)

deallocate(vr)

780

842

deallocate(vr)

deallocate(wi)

781

843

deallocate(wi)

deallocate(wr)

782

844

deallocate(wr)

deallocate(wc_a)

783

845

deallocate(wc_a)

return

784

846

return

end if

785

847

end if

786

848

! now fill in the results

787

849

! now fill in the results

i=1

788

850

i=1

do while(i<=d)

789

851

do while(i<=d)

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

790

852

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

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

791

853

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

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

792

854

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

else ! eigenvalue is complex

793

855

else ! eigenvalue is complex

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

794

856

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

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

795

857

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

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

796

858

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

i=i+1

797

859

i=i+1

end if

798

860

end if

i=i+1

799

861

i=i+1

enddo

800

862

enddo

801

863

deallocate(work)

802

864

deallocate(work)

deallocate(vr)

803

865

deallocate(vr)

deallocate(wi)

804

866

deallocate(wi)

deallocate(wr)

805

867

deallocate(wr)

deallocate(wc_a)

806

868

deallocate(wc_a)

807

869

end subroutine

808

870

end subroutine

809

871

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

810

872

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

!

811

873

!

! integer d (in) : matrice rank

812

874

! integer d (in) : matrice rank

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

813

875

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

! complex evala(d) (out) : eigenvalues

814

876

! complex evala(d) (out) : eigenvalues

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

815

877

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

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

816

878

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

!

817

879

!

! interfacing Lapack routine CGEEV

818

880

! interfacing Lapack routine CGEEV

implicit none

819

881

implicit none

integer, intent(in) :: d

820

882

integer, intent(in) :: d

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

821

883

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

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

822

884

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

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

823

885

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

integer, intent(out) :: status

824

886

integer, intent(out) :: status

825

887

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

826

888

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

integer :: info

827

889

integer :: info

integer :: lwork

828

890

integer :: lwork

complex, allocatable :: work(:)

829

891

complex, allocatable :: work(:)

complex :: twork(1)

830

892

complex :: twork(1)

real, allocatable :: rwork(:)

831

893

real, allocatable :: rwork(:)

complex :: vl ! unused but necessary for the call

832

894

complex :: vl ! unused but necessary for the call

833

895

status=1

834

896

status=1

835

897

! making a working copy of a

836

898

! making a working copy of a

allocate(wc_a(d,d))

837

899

allocate(wc_a(d,d))

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

838

900

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

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

839

901

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

840

902

841

903

! query optimal work size

842

904

! query optimal work size

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

843

905

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

lwork=int(twork(1))

844

906

lwork=int(twork(1))

allocate(work(lwork))

845

907

allocate(work(lwork))

allocate(rwork(2*d))

846

908

allocate(rwork(2*d))

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

847

909

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

848

910

if (info /= 0) then

849

911

if (info /= 0) then

status=0

850

912

status=0

end if

851

913

end if

deallocate(rwork)

852

914

deallocate(rwork)

deallocate(work)

853

915

deallocate(work)

deallocate(wc_a)

854

916

deallocate(wc_a)

855

917

end subroutine

856

918

end subroutine

857

919

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

858

920

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

!

859

921

!

! integer d (in) : matrice rank

860

922

! integer d (in) : matrice rank

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

861

923

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

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

862

924

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

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

863

925

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

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

864

926

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

!

865

927

!

! interfacing Lapack routine ZGEEV

866

928

! interfacing Lapack routine ZGEEV

implicit none

867

929

implicit none

integer, intent(in) :: d

868

930

integer, intent(in) :: d

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

869

931

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

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

870

932

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

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

871

933

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

integer, intent(out) :: status

872

934

integer, intent(out) :: status

873

935

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

874

936

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

integer :: info

875

937

integer :: info

integer :: lwork

876

938

integer :: lwork

double complex, allocatable :: work(:)

877

939

double complex, allocatable :: work(:)

double complex :: twork(1)

878

940

double complex :: twork(1)

double precision, allocatable :: rwork(:)

879

941

double precision, allocatable :: rwork(:)

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

880

942

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

881

943

status=1

882

944

status=1

883

945

! making a working copy of a

884

946

! making a working copy of a

allocate(wc_a(d,d))

885

947

allocate(wc_a(d,d))

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

886

948

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

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

887

949

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

888

950

! query optimal work size

889

951

! query optimal work size

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

890

952

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

lwork=int(twork(1))

891

953

lwork=int(twork(1))

allocate(work(lwork))

892

954

allocate(work(lwork))

allocate(rwork(2*d))

893

955

allocate(rwork(2*d))

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

894

956

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

895

957

if (info /= 0) then

896

958

if (info /= 0) then

status=0

897

959

status=0

end if

898

960

end if

deallocate(rwork)

899

961

deallocate(rwork)

deallocate(work)

900

962

deallocate(work)

deallocate(wc_a)

901

963

deallocate(wc_a)

902

964

end subroutine

903

965

end subroutine

904

966

905

967

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

906

968

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

!

907

969

!

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

908

970

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

!

909

971

!

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

910

972

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

911

973

subroutine fvn_s_find_interval(x,i,xdata,n)

912

974

subroutine fvn_s_find_interval(x,i,xdata,n)

implicit none

913

975

implicit none

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

914

976

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

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

915

977

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

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

916

978

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

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

917

979

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

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

918

980

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

! we will not exclude the upper limit

919

981

! we will not exclude the upper limit

! a simple dichotomy method is used

920

982

! a simple dichotomy method is used

921

983

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

922

984

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

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

923

985

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

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

924

986

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

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

925

987

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

926

988

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

927

989

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

928

990

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

929

991

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

! we will not exclude the upper limit

930

992

! we will not exclude the upper limit

if (x == xdata(n)) then

931

993

if (x == xdata(n)) then

i=n-1

932

994

i=n-1

return

933

995

return

end if

934

996

end if

935

997

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

936

998

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

if (x < xdata(1)) then

937

999

if (x < xdata(1)) then

i=0

938

1000

i=0

return

939

1001

return

end if

940

1002

end if

941

1003

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

942

1004

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

if (x > xdata(n)) then

943

1005

if (x > xdata(n)) then

i=n

944

1006

i=n

return

945

1007

return

end if

946

1008

end if

947

1009

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

948

1010

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

imin=0

949

1011

imin=0

imax=n+1

950

1012

imax=n+1

951

1013

do while((imax-imin) > 1)

952

1014

do while((imax-imin) > 1)

imoyen=(imax+imin)/2

953

1015

imoyen=(imax+imin)/2

if (x >= xdata(imoyen)) then

954

1016

if (x >= xdata(imoyen)) then

imin=imoyen

955

1017

imin=imoyen

else

956

1018

else

imax=imoyen

957

1019

imax=imoyen

end if

958

1020

end if

end do

959

1021

end do

960

1022

i=imin

961

1023

i=imin

962

1024

end subroutine

963

1025

end subroutine

964

1026

965

1027

subroutine fvn_d_find_interval(x,i,xdata,n)

966

1028

subroutine fvn_d_find_interval(x,i,xdata,n)

implicit none

967

1029

implicit none

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

968

1030

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

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

969

1031

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

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

970

1032

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

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

971

1033

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

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

972

1034

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

! we will not exclude the upper limit

973

1035

! we will not exclude the upper limit

! a simple dichotomy method is used

974

1036

! a simple dichotomy method is used

975

1037

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

976

1038

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

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

977

1039

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

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

978

1040

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

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

979

1041

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

980

1042

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

981

1043

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

982

1044

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

983

1045

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

! we will not exclude the upper limit

984

1046

! we will not exclude the upper limit

if (x == xdata(n)) then

985

1047

if (x == xdata(n)) then

i=n-1

986

1048

i=n-1

return

987

1049

return

end if

988

1050

end if

989

1051

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

990

1052

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

if (x < xdata(1)) then

991

1053

if (x < xdata(1)) then

i=0

992

1054

i=0

return

993

1055

return

end if

994

1056

end if

995

1057

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

996

1058

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

if (x > xdata(n)) then

997

1059

if (x > xdata(n)) then

i=n

998

1060

i=n

return

999

1061

return

end if

1000

1062

end if

1001

1063

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

1002

1064

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

imin=0

1003

1065

imin=0

imax=n+1

1004

1066

imax=n+1

1005

1067

do while((imax-imin) > 1)

1006

1068

do while((imax-imin) > 1)

imoyen=(imax+imin)/2

1007

1069

imoyen=(imax+imin)/2

if (x >= xdata(imoyen)) then

1008

1070

if (x >= xdata(imoyen)) then

imin=imoyen

1009

1071

imin=imoyen

else

1010

1072

else

imax=imoyen

1011

1073

imax=imoyen

end if

1012

1074

end if

end do

1013

1075

end do

1014

1076

i=imin

1015

1077

i=imin

1016

1078

end subroutine

1017

1079

end subroutine

1018

1080

1019

1081

function fvn_s_quad_interpol(x,n,xdata,ydata)

1020

1082

function fvn_s_quad_interpol(x,n,xdata,ydata)

implicit none

1021

1083

implicit none

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

1022

1084

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

! and values, using a quadratic interpolation

1023

1085

! and values, using a quadratic interpolation

! xdata must be increasingly ordered

1024

1086

! xdata must be increasingly ordered

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

1025

1087

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

! otherwise extrapolation is done

1026

1088

! otherwise extrapolation is done

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

1027

1089

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

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

1028

1090

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

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

1029

1091

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

real(kind=4) :: fvn_s_quad_interpol

1030

1092

real(kind=4) :: fvn_s_quad_interpol

1031

1093

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

1032

1094

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

real(kind=4) :: p

1033

1095

real(kind=4) :: p

1034

1096

call fvn_s_find_interval(x,iinf,xdata,n)

1035

1097

call fvn_s_find_interval(x,iinf,xdata,n)

1036

1098

! Settings for extrapolation

1037

1099

! Settings for extrapolation

if (iinf==0) then

1038

1100

if (iinf==0) then

! TODO -> Lower bound extrapolation warning

1039

1101

! TODO -> Lower bound extrapolation warning

iinf=1

1040

1102

iinf=1

end if

1041

1103

end if

1042

1104

if (iinf==n) then

1043

1105

if (iinf==n) then

! TODO -> Higher bound extrapolation warning

1044

1106

! TODO -> Higher bound extrapolation warning

iinf=n-1

1045

1107

iinf=n-1

end if

1046

1108

end if

1047

1109

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

1048

1110

! 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

1049

1111

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

if (iinf==1) then

1050

1112

if (iinf==1) then

base=0

1051

1113

base=0

else

1052

1114

else

base=iinf-2

1053

1115

base=iinf-2

end if

1054

1116

end if

1055

1117

! The three points we will use are :

1056

1118

! The three points we will use are :

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

1057

1119

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

1058

1120

! Straight forward Lagrange polynomial

1059

1121

! Straight forward Lagrange polynomial

fvn_s_quad_interpol=0.

1060

1122

fvn_s_quad_interpol=0.

do i=1,3

1061

1123

do i=1,3

! polynome i

1062

1124

! polynome i

p=ydata(base+i)

1063

1125

p=ydata(base+i)

do j=1,3

1064

1126

do j=1,3

if (j /= i) then

1065

1127

if (j /= i) then

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

1066

1128

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

end if

1067

1129

end if

end do

1068

1130

end do

fvn_s_quad_interpol=fvn_s_quad_interpol+p

1069

1131

fvn_s_quad_interpol=fvn_s_quad_interpol+p

end do

1070

1132

end do

1071

1133

end function

1072

1134

end function

1073

1135

1074

1136

function fvn_d_quad_interpol(x,n,xdata,ydata)

1075

1137

function fvn_d_quad_interpol(x,n,xdata,ydata)

implicit none

1076

1138

implicit none

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

1077

1139

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

! and values, using a quadratic interpolation

1078

1140

! and values, using a quadratic interpolation

! xdata must be increasingly ordered

1079

1141

! xdata must be increasingly ordered

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

1080

1142

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

! otherwise extrapolation is done

1081

1143

! otherwise extrapolation is done

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

1082

1144

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

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

1083

1145

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

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

1084

1146

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

real(kind=8) :: fvn_d_quad_interpol

1085

1147

real(kind=8) :: fvn_d_quad_interpol

1086

1148

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

1087

1149

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

real(kind=8) :: p

1088

1150

real(kind=8) :: p

1089

1151

call fvn_d_find_interval(x,iinf,xdata,n)

1090

1152

call fvn_d_find_interval(x,iinf,xdata,n)

1091

1153

! Settings for extrapolation

1092

1154

! Settings for extrapolation

if (iinf==0) then

1093

1155

if (iinf==0) then

! TODO -> Lower bound extrapolation warning

1094

1156

! TODO -> Lower bound extrapolation warning

iinf=1

1095

1157

iinf=1

end if

1096

1158

end if

1097

1159

if (iinf==n) then

1098

1160

if (iinf==n) then

! TODO Higher bound extrapolation warning

1099

1161

! TODO Higher bound extrapolation warning

iinf=n-1

1100

1162

iinf=n-1

end if

1101

1163

end if

1102

1164

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

1103

1165

! 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

1104

1166

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

if (iinf==1) then

1105

1167

if (iinf==1) then

base=0

1106

1168

base=0

else

1107

1169

else

base=iinf-2

1108

1170

base=iinf-2

end if

1109

1171

end if

1110

1172

! The three points we will use are :

1111

1173

! The three points we will use are :

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

1112

1174

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

1113

1175

! Straight forward Lagrange polynomial

1114

1176

! Straight forward Lagrange polynomial

fvn_d_quad_interpol=0.

1115

1177

fvn_d_quad_interpol=0.

do i=1,3

1116

1178

do i=1,3

! polynome i

1117

1179

! polynome i

p=ydata(base+i)

1118

1180

p=ydata(base+i)

do j=1,3

1119

1181

do j=1,3

if (j /= i) then

1120

1182

if (j /= i) then

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

1121

1183

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

end if

1122

1184

end if

end do

1123

1185

end do

fvn_d_quad_interpol=fvn_d_quad_interpol+p

1124

1186

fvn_d_quad_interpol=fvn_d_quad_interpol+p

end do

1125

1187

end do

1126

1188

end function

1127

1189

end function

1128

1190

1129

1191

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

1130

1192

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

implicit none

1131

1193

implicit none

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

1132

1194

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

! set of points and values, using a quadratic interpolation

1133

1195

! set of points and values, using a quadratic interpolation

! xdata and ydata must be increasingly ordered

1134

1196

! xdata and ydata must be increasingly ordered

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

1135

1197

! 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

1136

1198

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

! otherwise extrapolation is done

1137

1199

! otherwise extrapolation is done

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

1138

1200

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

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

1139

1201

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

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

1140

1202

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

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

1141

1203

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

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

1142

1204

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

real(kind=4) :: fvn_s_quad_2d_interpol

1143

1205

real(kind=4) :: fvn_s_quad_2d_interpol

1144

1206

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

1145

1207

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

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

1146

1208

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

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

1147

1209

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

1148

1210

call fvn_s_find_interval(x,ixinf,xdata,nx)

1149

1211

call fvn_s_find_interval(x,ixinf,xdata,nx)

call fvn_s_find_interval(y,iyinf,ydata,ny)

1150

1212

call fvn_s_find_interval(y,iyinf,ydata,ny)

1151

1213

! Settings for extrapolation

1152

1214

! Settings for extrapolation

if (ixinf==0) then

1153

1215

if (ixinf==0) then

! TODO -> Lower x bound extrapolation warning

1154

1216

! TODO -> Lower x bound extrapolation warning

ixinf=1

1155

1217

ixinf=1

end if

1156

1218

end if

1157

1219

if (ixinf==nx) then

1158

1220

if (ixinf==nx) then

! TODO -> Higher x bound extrapolation warning

1159

1221

! TODO -> Higher x bound extrapolation warning

ixinf=nx-1

1160

1222

ixinf=nx-1

end if

1161

1223

end if

1162

1224

if (iyinf==0) then

1163

1225

if (iyinf==0) then

! TODO -> Lower y bound extrapolation warning

1164

1226

! TODO -> Lower y bound extrapolation warning

iyinf=1

1165

1227

iyinf=1

end if

1166

1228

end if

1167

1229

if (iyinf==ny) then

1168

1230

if (iyinf==ny) then

! TODO -> Higher y bound extrapolation warning

1169

1231

! TODO -> Higher y bound extrapolation warning

iyinf=ny-1

1170

1232

iyinf=ny-1

end if

1171

1233

end if

1172

1234

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

1173

1235

! 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

1174

1236

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

if (ixinf==1) then

1175

1237

if (ixinf==1) then

basex=0

1176

1238

basex=0

else

1177

1239

else

basex=ixinf-2

1178

1240

basex=ixinf-2

end if

1179

1241

end if

1180

1242

if (iyinf==1) then

1181

1243

if (iyinf==1) then

basey=0

1182

1244

basey=0

else

1183

1245

else

basey=iyinf-2

1184

1246

basey=iyinf-2

end if

1185

1247

end if

1186

1248

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

1187

1249

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

! stored in ztmp(1:3)

1188

1250

! stored in ztmp(1:3)

do i=1,3

1189

1251

do i=1,3

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

1190

1252

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

end do

1191

1253

end do

1192

1254

! Then we make an interpolation for y using previous interpolations

1193

1255

! 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)

1194

1256

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

end function

1195

1257

end function

1196

1258

1197

1259

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

1198

1260

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

implicit none

1199

1261

implicit none

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

1200

1262

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

! set of points and values, using a quadratic interpolation

1201

1263

! set of points and values, using a quadratic interpolation

! xdata and ydata must be increasingly ordered

1202

1264

! xdata and ydata must be increasingly ordered

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

1203

1265

! 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

1204

1266

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

! otherwise extrapolation is done

1205

1267

! otherwise extrapolation is done

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

1206

1268

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

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

1207

1269

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

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

1208

1270

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

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

1209

1271

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

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

1210

1272

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

real(kind=8) :: fvn_d_quad_2d_interpol

1211

1273

real(kind=8) :: fvn_d_quad_2d_interpol

1212

1274

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

1213

1275

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

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

1214

1276

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

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

1215

1277

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

1216

1278

call fvn_d_find_interval(x,ixinf,xdata,nx)

1217

1279

call fvn_d_find_interval(x,ixinf,xdata,nx)

call fvn_d_find_interval(y,iyinf,ydata,ny)

1218

1280

call fvn_d_find_interval(y,iyinf,ydata,ny)

1219

1281

! Settings for extrapolation

1220

1282

! Settings for extrapolation

if (ixinf==0) then

1221

1283

if (ixinf==0) then

! TODO -> Lower x bound extrapolation warning

1222

1284

! TODO -> Lower x bound extrapolation warning

ixinf=1

1223

1285

ixinf=1

end if

1224

1286

end if

1225

1287

if (ixinf==nx) then

1226

1288

if (ixinf==nx) then

! TODO -> Higher x bound extrapolation warning

1227

1289

! TODO -> Higher x bound extrapolation warning

ixinf=nx-1

1228

1290

ixinf=nx-1

end if

1229

1291

end if

1230

1292

if (iyinf==0) then

1231

1293

if (iyinf==0) then

! TODO -> Lower y bound extrapolation warning

1232

1294

! TODO -> Lower y bound extrapolation warning

iyinf=1

1233

1295

iyinf=1

end if

1234

1296

end if

1235

1297

if (iyinf==ny) then

1236

1298

if (iyinf==ny) then

! TODO -> Higher y bound extrapolation warning

1237

1299

! TODO -> Higher y bound extrapolation warning

iyinf=ny-1

1238

1300

iyinf=ny-1

end if

1239

1301

end if

1240

1302

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

1241

1303

! 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

1242

1304

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

if (ixinf==1) then

1243

1305

if (ixinf==1) then

basex=0

1244

1306

basex=0

else

1245

1307

else

basex=ixinf-2

1246

1308

basex=ixinf-2

end if

1247

1309

end if

1248

1310

if (iyinf==1) then

1249

1311

if (iyinf==1) then

basey=0

1250

1312

basey=0

else

1251

1313

else

basey=iyinf-2

1252

1314

basey=iyinf-2

end if

1253

1315

end if

1254

1316

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

1255

1317

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

! stored in ztmp(1:3)

1256

1318

! stored in ztmp(1:3)

do i=1,3

1257

1319

do i=1,3

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

1258

1320

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

end do

1259

1321

end do

1260

1322

! Then we make an interpolation for y using previous interpolations

1261

1323

! 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)

1262

1324

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

end function

1263

1325

end function

1264

1326

1265

1327

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

1266

1328

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

implicit none

1267

1329

implicit none

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

1268

1330

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

! set of points and values, using a quadratic interpolation

1269

1331

! set of points and values, using a quadratic interpolation

! xdata, ydata and zdata must be increasingly ordered

1270

1332

! xdata, ydata and zdata must be increasingly ordered

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

1271

1333

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

! perform an interpolation, otherwise extrapolation is done

1272

1334

! perform an interpolation, otherwise extrapolation is done

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

1273

1335

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

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

1274

1336

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

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

1275

1337

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

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

1276

1338

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

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

1277

1339

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

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

1278

1340

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

real(kind=4) :: fvn_s_quad_3d_interpol

1279

1341

real(kind=4) :: fvn_s_quad_3d_interpol

1280

1342

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

1281

1343

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

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

1282

1344

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

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

1283

1345

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

1284

1346

call fvn_s_find_interval(x,ixinf,xdata,nx)

1285

1347

call fvn_s_find_interval(x,ixinf,xdata,nx)

call fvn_s_find_interval(y,iyinf,ydata,ny)

1286

1348

call fvn_s_find_interval(y,iyinf,ydata,ny)

call fvn_s_find_interval(z,izinf,zdata,nz)

1287

1349

call fvn_s_find_interval(z,izinf,zdata,nz)

1288

1350

! Settings for extrapolation

1289

1351

! Settings for extrapolation

if (ixinf==0) then

1290

1352

if (ixinf==0) then

! TODO -> Lower x bound extrapolation warning

1291

1353

! TODO -> Lower x bound extrapolation warning

ixinf=1

1292

1354

ixinf=1

end if

1293

1355

end if

1294

1356

if (ixinf==nx) then

1295

1357

if (ixinf==nx) then

! TODO -> Higher x bound extrapolation warning

1296

1358

! TODO -> Higher x bound extrapolation warning

ixinf=nx-1

1297

1359

ixinf=nx-1

end if

1298

1360

end if

1299

1361

if (iyinf==0) then

1300

1362

if (iyinf==0) then

! TODO -> Lower y bound extrapolation warning

1301

1363

! TODO -> Lower y bound extrapolation warning

iyinf=1

1302

1364

iyinf=1

end if

1303

1365

end if

1304

1366

if (iyinf==ny) then

1305

1367

if (iyinf==ny) then

! TODO -> Higher y bound extrapolation warning

1306

1368

! TODO -> Higher y bound extrapolation warning

iyinf=ny-1

1307

1369

iyinf=ny-1

end if

1308

1370

end if

1309

1371

if (izinf==0) then

1310

1372

if (izinf==0) then

! TODO -> Lower z bound extrapolation warning

1311

1373

! TODO -> Lower z bound extrapolation warning

izinf=1

1312

1374

izinf=1

end if

1313

1375

end if

1314

1376

if (izinf==nz) then

1315

1377

if (izinf==nz) then

! TODO -> Higher z bound extrapolation warning

1316

1378

! TODO -> Higher z bound extrapolation warning

izinf=nz-1

1317

1379

izinf=nz-1

end if

1318

1380

end if

1319

1381

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

1320

1382

! 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

1321

1383

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

if (ixinf==1) then

1322

1384

if (ixinf==1) then

basex=0

1323

1385

basex=0

else

1324

1386

else

basex=ixinf-2

1325

1387

basex=ixinf-2

end if

1326

1388

end if

1327

1389

if (iyinf==1) then

1328

1390

if (iyinf==1) then

basey=0

1329

1391

basey=0

else

1330

1392

else

basey=iyinf-2

1331

1393

basey=iyinf-2

end if

1332

1394

end if

1333

1395

if (izinf==1) then

1334

1396

if (izinf==1) then

basez=0

1335

1397

basez=0

else

1336

1398

else

basez=izinf-2

1337

1399

basez=izinf-2

end if

1338

1400

end if

1339

1401

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

1340

1402

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

! results are stored in ttmp

1341

1403

! results are stored in ttmp

do i=1,3

1342

1404

do i=1,3

do j=1,3

1343

1405

do j=1,3

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

1344

1406

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

end do

1345

1407

end do

1346

1408

end do

1347

1409

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

1348

1410

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

fvn_s_quad_3d_interpol=fvn_s_quad_2d_interpol(y,z, &

1349

1411

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)

1350

1412

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

end function

1351

1413

end function

1352

1414

1353

1415

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

1354

1416

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

implicit none

1355

1417

implicit none

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

1356

1418

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

! set of points and values, using a quadratic interpolation

1357

1419

! set of points and values, using a quadratic interpolation

! xdata, ydata and zdata must be increasingly ordered

1358

1420

! xdata, ydata and zdata must be increasingly ordered

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

1359

1421

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

! perform an interpolation, otherwise extrapolation is done

1360

1422

! perform an interpolation, otherwise extrapolation is done

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

1361

1423

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

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

1362

1424

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

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

1363

1425

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

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

1364

1426

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

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

1365

1427

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

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

1366

1428

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

real(kind=8) :: fvn_d_quad_3d_interpol

1367

1429

real(kind=8) :: fvn_d_quad_3d_interpol

1368

1430

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

1369

1431

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

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

1370

1432

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

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

1371

1433

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

1372

1434

call fvn_d_find_interval(x,ixinf,xdata,nx)

1373

1435

call fvn_d_find_interval(x,ixinf,xdata,nx)

call fvn_d_find_interval(y,iyinf,ydata,ny)

1374

1436

call fvn_d_find_interval(y,iyinf,ydata,ny)

call fvn_d_find_interval(z,izinf,zdata,nz)

1375

1437

call fvn_d_find_interval(z,izinf,zdata,nz)

1376

1438

! Settings for extrapolation

1377

1439

! Settings for extrapolation

if (ixinf==0) then

1378

1440

if (ixinf==0) then

! TODO -> Lower x bound extrapolation warning

1379

1441

! TODO -> Lower x bound extrapolation warning

ixinf=1

1380

1442

ixinf=1

end if

1381

1443

end if

1382

1444

if (ixinf==nx) then

1383

1445

if (ixinf==nx) then

! TODO -> Higher x bound extrapolation warning

1384

1446

! TODO -> Higher x bound extrapolation warning

ixinf=nx-1

1385

1447

ixinf=nx-1

end if

1386

1448

end if

1387

1449

if (iyinf==0) then

1388

1450

if (iyinf==0) then

! TODO -> Lower y bound extrapolation warning

1389

1451

! TODO -> Lower y bound extrapolation warning

iyinf=1

1390

1452

iyinf=1

end if

1391

1453

end if

1392

1454

if (iyinf==ny) then

1393

1455

if (iyinf==ny) then

! TODO -> Higher y bound extrapolation warning

1394

1456

! TODO -> Higher y bound extrapolation warning

iyinf=ny-1

1395

1457

iyinf=ny-1

end if

1396

1458

end if

1397

1459

if (izinf==0) then

1398

1460

if (izinf==0) then

! TODO -> Lower z bound extrapolation warning

1399

1461

! TODO -> Lower z bound extrapolation warning

izinf=1

1400

1462

izinf=1

end if

1401

1463

end if

1402

1464

if (izinf==nz) then

1403

1465

if (izinf==nz) then

! TODO -> Higher z bound extrapolation warning

1404

1466

! TODO -> Higher z bound extrapolation warning

izinf=nz-1

1405

1467

izinf=nz-1

end if

1406

1468

end if

1407

1469

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

1408

1470

! 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

1409

1471

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

if (ixinf==1) then

1410

1472

if (ixinf==1) then

basex=0

1411

1473

basex=0

else

1412

1474

else

basex=ixinf-2

1413

1475

basex=ixinf-2

end if

1414

1476

end if

1415

1477

if (iyinf==1) then

1416

1478

if (iyinf==1) then

basey=0

1417

1479

basey=0

else

1418

1480

else

basey=iyinf-2

1419

1481

basey=iyinf-2

end if

1420

1482

end if

1421

1483

if (izinf==1) then

1422

1484

if (izinf==1) then

basez=0

1423

1485

basez=0

else

1424

1486

else

basez=izinf-2

1425

1487

basez=izinf-2

end if

1426

1488

end if

1427

1489

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

1428

1490

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

! results are stored in ttmp

1429

1491

! results are stored in ttmp

do i=1,3

1430

1492

do i=1,3

do j=1,3

1431

1493

do j=1,3

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

1432

1494

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

end do

1433

1495

end do

1434

1496

end do

1435

1497

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

1436

1498

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

fvn_d_quad_3d_interpol=fvn_d_quad_2d_interpol(y,z, &

1437

1499

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)

1438

1500

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

end function

1439

1501

end function

1440

1502

1441

1503

1442

1504

1443

1505

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

1444

1506

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

!

1445

1507

!

! Akima spline interpolation and spline evaluation

1446

1508

! Akima spline interpolation and spline evaluation

!

1447

1509

!

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

1448

1510

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

1449

1511

! Single precision

1450

1512

! Single precision

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

1451

1513

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

implicit none

1452

1514

implicit none

integer, intent(in) :: n

1453

1515

integer, intent(in) :: n

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

1454

1516

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

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

1455

1517

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

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

1456

1518

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

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

1457

1519

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

1458

1520

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

1459

1521

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

real :: wi_1,wi

1460

1522

real :: wi_1,wi

integer :: i

1461

1523

integer :: i

real :: dx,a,b

1462

1524

real :: dx,a,b

1463

1525

! br is just a copy of x

1464

1526

! br is just a copy of x

br(:)=x(:)

1465

1527

br(:)=x(:)

1466

1528

allocate(var(n))

1467

1529

allocate(var(n))

allocate(z(n))

1468

1530

allocate(z(n))

! evaluate the variations

1469

1531

! evaluate the variations

do i=1, n-1

1470

1532

do i=1, n-1

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

1471

1533

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

end do

1472

1534

end do

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

1473

1535

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

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

1474

1536

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

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

1475

1537

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

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

1476

1538

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

1477

1539

do i = 1, n

1478

1540

do i = 1, n

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

1479

1541

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

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

1480

1542

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

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

1481

1543

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

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

1482

1544

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

else

1483

1545

else

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

1484

1546

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

end if

1485

1547

end if

end do

1486

1548

end do

1487

1549

do i=1, n-1

1488

1550

do i=1, n-1

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

1489

1551

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

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

1490

1552

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

1491

1553

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

co(1,i)=y(i)

1492

1554

co(1,i)=y(i)

co(2,i)=z(i)

1493

1555

co(2,i)=z(i)

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

1494

1556

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

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

1495

1557

!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

1496

1558

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 !

1497

1559

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

1498

1560

! 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 ...

1499

1561

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

end do

1500

1562

end do

co(1,n)=y(n)

1501

1563

co(1,n)=y(n)

co(2,n)=z(n)

1502

1564

co(2,n)=z(n)

co(3,n)=0.e0

1503

1565

co(3,n)=0.e0

co(4,n)=0.e0

1504

1566

co(4,n)=0.e0

1505

1567

deallocate(z)

1506

1568

deallocate(z)

deallocate(var)

1507

1569

deallocate(var)

1508

1570

end subroutine

1509

1571

end subroutine

1510

1572

! Double precision

1511

1573

! Double precision

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

1512

1574

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

1513

1575

implicit none

1514

1576

implicit none

integer, intent(in) :: n

1515

1577

integer, intent(in) :: n

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

1516

1578

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

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

1517

1579

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

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

1518

1580

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

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

1519

1581

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

1520

1582

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

1521

1583

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

double precision :: wi_1,wi

1522

1584

double precision :: wi_1,wi

integer :: i

1523

1585

integer :: i

double precision :: dx,a,b

1524

1586

double precision :: dx,a,b

1525

1587

! br is just a copy of x

1526

1588

! br is just a copy of x

br(:)=x(:)

1527

1589

br(:)=x(:)

1528

1590

allocate(var(n))

1529

1591

allocate(var(n))

allocate(z(n))

1530

1592

allocate(z(n))

! evaluate the variations

1531

1593

! evaluate the variations

do i=1, n-1

1532

1594

do i=1, n-1

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

1533

1595

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

end do

1534

1596

end do

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

1535

1597

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

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

1536

1598

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

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

1537

1599

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

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

1538

1600

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

1539

1601

do i = 1, n

1540

1602

do i = 1, n

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

1541

1603

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

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

1542

1604

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

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

1543

1605

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

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

1544

1606

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

else

1545

1607

else

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

1546

1608

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

end if

1547

1609

end if

end do

1548

1610

end do

1549

1611

do i=1, n-1

1550

1612

do i=1, n-1

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

1551

1613

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

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

1552

1614

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

1553

1615

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

co(1,i)=y(i)

1554

1616

co(1,i)=y(i)

co(2,i)=z(i)

1555

1617

co(2,i)=z(i)

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

1556

1618

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

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

1557

1619

!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

1558

1620

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 !

1559

1621

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

1560

1622

! 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 ...

1561

1623

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

end do

1562

1624

end do

co(1,n)=y(n)

1563

1625

co(1,n)=y(n)

co(2,n)=z(n)

1564

1626

co(2,n)=z(n)

co(3,n)=0.d0

1565

1627

co(3,n)=0.d0

co(4,n)=0.d0

1566

1628

co(4,n)=0.d0

1567

1629

deallocate(z)

1568

1630

deallocate(z)

deallocate(var)

1569

1631

deallocate(var)

1570

1632

end subroutine

1571

1633

end subroutine

1572

1634

!

1573

1635

!

! Single precision spline evaluation

1574

1636

! Single precision spline evaluation

!

1575

1637

!

function fvn_s_spline_eval(x,n,br,co)

1576

1638

function fvn_s_spline_eval(x,n,br,co)

implicit none

1577

1639

implicit none

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

1578

1640

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

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

1579

1641

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

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

1580

1642

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

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

1581

1643

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

real :: fvn_s_spline_eval

1582

1644

real :: fvn_s_spline_eval

1583

1645

integer :: i

1584

1646

integer :: i

real :: dx

1585

1647

real :: dx

1586

1648

if (x<=br(1)) then

1587

1649

if (x<=br(1)) then

i=1

1588

1650

i=1

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

1589

1651

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

i=n

1590

1652

i=n

else

1591

1653

else

i=1

1592

1654

i=1

do while(x>=br(i))

1593

1655

do while(x>=br(i))

i=i+1

1594

1656

i=i+1

end do

1595

1657

end do

i=i-1

1596

1658

i=i-1

end if

1597

1659

end if

dx=x-br(i)

1598

1660

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

1599

1661

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

1600

1662

end function

1601

1663

end function

1602

1664

! Double precision spline evaluation

1603

1665

! Double precision spline evaluation

function fvn_d_spline_eval(x,n,br,co)

1604

1666

function fvn_d_spline_eval(x,n,br,co)

implicit none

1605

1667

implicit none

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

1606

1668

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

1607

1669

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

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

1608

1670

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

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

1609

1671

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

double precision :: fvn_d_spline_eval

1610

1672

double precision :: fvn_d_spline_eval

1611

1673

integer :: i

1612

1674

integer :: i

double precision :: dx

1613

1675

double precision :: dx

1614

1676

1615

1677

if (x<=br(1)) then

1616

1678

if (x<=br(1)) then

i=1

1617

1679

i=1

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

1618

1680

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

i=n

1619

1681

i=n

else

1620

1682

else

i=1

1621

1683

i=1

do while(x>=br(i))

1622

1684

do while(x>=br(i))

i=i+1

1623

1685

i=i+1

end do

1624

1686

end do

i=i-1

1625

1687

i=i-1

end if

1626

1688

end if

1627

1689

dx=x-br(i)

1628

1690

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

1629

1691

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

1630

1692

end function

1631

1693

end function

1632

1694

1633

1695

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

1634

1696

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

!

1635

1697

!

! Least square problem

1636

1698

! Least square problem

!

1637

1699

!

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

1638

1700

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

!

1639

1701

!

1640

1702

!

1641

1703

1642

1704

1643

1705

1644

1706

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

1645

1707

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

!

1646

1708

!

! Least square polynomial fitting

1647

1709

! Least square polynomial fitting

!

1648

1710

!

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

1649

1711

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

! for a set of coordinates.

1650

1712

! for a set of coordinates.

!

1651

1713

!

! The degree must be lower than the number of points

1652

1714

! The degree must be lower than the number of points

!

1653

1715

!

! np (in) : number of points

1654

1716

! np (in) : number of points

! x(np) (in) : x data

1655

1717

! x(np) (in) : x data

! y(np) (in) : y data

1656

1718

! y(np) (in) : y data

! deg (in) : polynomial's degree

1657

1719

! deg (in) : polynomial's degree

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

1658

1720

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

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

1659

1721

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

!

1660

1722

!

implicit none

1661

1723

implicit none

1662

1724

integer, intent(in) :: np,deg

1663

1725

integer, intent(in) :: np,deg

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

1664

1726

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

real(kind=8), intent(out), dimension(deg+1) :: coeff

1665

1727

real(kind=8), intent(out), dimension(deg+1) :: coeff

integer, intent(out) :: status

1666

1728

integer, intent(out) :: status

1667

1729

real(kind=8), allocatable, dimension(:,:) :: mat,bmat

1668

1730

real(kind=8), allocatable, dimension(:,:) :: mat,bmat

real(kind=8),dimension(:),allocatable :: work

1669

1731

real(kind=8),dimension(:),allocatable :: work

real(kind=8),dimension(1) :: twork

1670

1732

real(kind=8),dimension(1) :: twork

integer :: lwork,info

1671

1733

integer :: lwork,info

1672

1734

integer :: i,j

1673

1735

integer :: i,j

1674

1736

status=1

1675

1737

status=1

allocate(mat(np,deg+1),bmat(np,1))

1676

1738

allocate(mat(np,deg+1),bmat(np,1))

1677

1739

! Design matrix valorisation

1678

1740

! Design matrix valorisation

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1679

1741

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1680

1742

! second member valorisation

1681

1743

! second member valorisation

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1682

1744

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1683

1745

! query workspace size

1684

1746

! query workspace size

call dgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)

1685

1747

call dgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)

lwork=twork(1)

1686

1748

lwork=twork(1)

allocate(work(int(lwork)))

1687

1749

allocate(work(int(lwork)))

! real call

1688

1750

! real call

call dgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)

1689

1751

call dgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)

1690

1752

if (info /= 0) then

1691

1753

if (info /= 0) then

status=0

1692

1754

status=0

end if

1693

1755

end if

1694

1756

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1695

1757

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1696

1758

deallocate(work)

1697

1759

deallocate(work)

deallocate(mat,bmat)

1698

1760

deallocate(mat,bmat)

end subroutine

1699

1761

end subroutine

1700

1762

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

1701

1763

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

!

1702

1764

!

! Least square polynomial fitting

1703

1765

! Least square polynomial fitting

!

1704

1766

!

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

1705

1767

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

! for a set of coordinates.

1706

1768

! for a set of coordinates.

!

1707

1769

!

! The degree must be lower than the number of points

1708

1770

! The degree must be lower than the number of points

!

1709

1771

!

! np (in) : number of points

1710

1772

! np (in) : number of points

! x(np) (in) : x data

1711

1773

! x(np) (in) : x data

! y(np) (in) : y data

1712

1774

! y(np) (in) : y data

! deg (in) : polynomial's degree

1713

1775

! deg (in) : polynomial's degree

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

1714

1776

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

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

1715

1777

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

!

1716

1778

!

implicit none

1717

1779

implicit none

1718

1780

integer, intent(in) :: np,deg

1719

1781

integer, intent(in) :: np,deg

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

1720

1782

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

real(kind=4), intent(out), dimension(deg+1) :: coeff

1721

1783

real(kind=4), intent(out), dimension(deg+1) :: coeff

integer, intent(out) :: status

1722

1784

integer, intent(out) :: status

1723

1785

real(kind=4), allocatable, dimension(:,:) :: mat,bmat

1724

1786

real(kind=4), allocatable, dimension(:,:) :: mat,bmat

real(kind=4),dimension(:),allocatable :: work

1725

1787

real(kind=4),dimension(:),allocatable :: work

real(kind=4),dimension(1) :: twork

1726

1788

real(kind=4),dimension(1) :: twork

integer :: lwork,info

1727

1789

integer :: lwork,info

1728

1790

integer :: i,j

1729

1791

integer :: i,j

1730

1792

status=1

1731

1793

status=1

allocate(mat(np,deg+1),bmat(np,1))

1732

1794

allocate(mat(np,deg+1),bmat(np,1))

1733

1795

! Design matrix valorisation

1734

1796

! Design matrix valorisation

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1735

1797

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1736

1798

! second member valorisation

1737

1799

! second member valorisation

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1738

1800

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1739

1801

! query workspace size

1740

1802

! query workspace size

call sgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)

1741

1803

call sgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)

lwork=twork(1)

1742

1804

lwork=twork(1)

allocate(work(int(lwork)))

1743

1805

allocate(work(int(lwork)))

! real call

1744

1806

! real call

call sgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)

1745

1807

call sgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)

1746

1808

if (info /= 0) then

1747

1809

if (info /= 0) then

status=0

1748

1810

status=0

end if

1749

1811

end if

1750

1812

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1751

1813

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1752

1814

deallocate(work)

1753

1815

deallocate(work)

deallocate(mat,bmat)

1754

1816

deallocate(mat,bmat)

end subroutine

1755

1817

end subroutine

1756

1818

1757

1819

1758

1820

1759

1821

1760

1822

1761

1823

1762

1824

1763

1825

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

1764

1826

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

!

1765

1827

!

! Least square polynomial fitting

1766

1828

! Least square polynomial fitting

!

1767

1829

!

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

1768

1830

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

! for a set of coordinates.

1769

1831

! for a set of coordinates.

!

1770

1832

!

! The degree must be lower than the number of points

1771

1833

! The degree must be lower than the number of points

!

1772

1834

!

! np (in) : number of points

1773

1835

! np (in) : number of points

! x(np) (in) : x data

1774

1836

! x(np) (in) : x data

! y(np) (in) : y data

1775

1837

! y(np) (in) : y data

! deg (in) : polynomial's degree

1776

1838

! deg (in) : polynomial's degree

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

1777

1839

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

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

1778

1840

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

!

1779

1841

!

implicit none

1780

1842

implicit none

1781

1843

integer, intent(in) :: np,deg

1782

1844

integer, intent(in) :: np,deg

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

1783

1845

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

real(kind=8), intent(out), dimension(deg+1) :: coeff

1784

1846

real(kind=8), intent(out), dimension(deg+1) :: coeff

integer, intent(out) :: status

1785

1847

integer, intent(out) :: status

1786

1848

real(kind=8), allocatable, dimension(:,:) :: mat,bmat

1787

1849

real(kind=8), allocatable, dimension(:,:) :: mat,bmat

real(kind=8),dimension(:),allocatable :: work,singval

1788

1850

real(kind=8),dimension(:),allocatable :: work,singval

real(kind=8),dimension(1) :: twork

1789

1851

real(kind=8),dimension(1) :: twork

integer :: lwork,info,rank

1790

1852

integer :: lwork,info,rank

1791

1853

integer :: i,j

1792

1854

integer :: i,j

1793

1855

status=1

1794

1856

status=1

allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))

1795

1857

allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))

1796

1858

! Design matrix valorisation

1797

1859

! Design matrix valorisation

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1798

1860

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1799

1861

! second member valorisation

1800

1862

! second member valorisation

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1801

1863

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1802

1864

! query workspace size

1803

1865

! query workspace size

call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)

1804

1866

call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)

lwork=twork(1)

1805

1867

lwork=twork(1)

allocate(work(int(lwork)))

1806

1868

allocate(work(int(lwork)))

! real call

1807

1869

! real call

call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)

1808

1870

call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)

1809

1871

if (info /= 0) then

1810

1872

if (info /= 0) then

status=0

1811

1873

status=0

end if

1812

1874

end if

1813

1875

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1814

1876

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1815

1877

deallocate(work)

1816

1878

deallocate(work)

deallocate(mat,bmat,singval)

1817

1879

deallocate(mat,bmat,singval)

end subroutine

1818

1880

end subroutine

1819

1881

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

1820

1882

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

!

1821

1883

!

! Least square polynomial fitting

1822

1884

! Least square polynomial fitting

!

1823

1885

!

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

1824

1886

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

! for a set of coordinates.

1825

1887

! for a set of coordinates.

!

1826

1888

!

! The degree must be lower than the number of points

1827

1889

! The degree must be lower than the number of points

!

1828

1890

!

! np (in) : number of points

1829

1891

! np (in) : number of points

! x(np) (in) : x data

1830

1892

! x(np) (in) : x data

! y(np) (in) : y data

1831

1893

! y(np) (in) : y data

! deg (in) : polynomial's degree

1832

1894

! deg (in) : polynomial's degree

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

1833

1895

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

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

1834

1896

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

!

1835

1897

!

implicit none

1836

1898

implicit none

1837

1899

integer, intent(in) :: np,deg

1838

1900

integer, intent(in) :: np,deg

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

1839

1901

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

real(kind=4), intent(out), dimension(deg+1) :: coeff

1840

1902

real(kind=4), intent(out), dimension(deg+1) :: coeff

integer, intent(out) :: status

1841

1903

integer, intent(out) :: status

1842

1904

real(kind=4), allocatable, dimension(:,:) :: mat,bmat

1843

1905

real(kind=4), allocatable, dimension(:,:) :: mat,bmat

real(kind=4),dimension(:),allocatable :: work,singval

1844

1906

real(kind=4),dimension(:),allocatable :: work,singval

real(kind=4),dimension(1) :: twork

1845

1907

real(kind=4),dimension(1) :: twork

integer :: lwork,info,rank

1846

1908

integer :: lwork,info,rank

1847

1909

integer :: i,j

1848

1910

integer :: i,j

1849

1911

status=1

1850

1912

status=1

allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))

1851

1913

allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))

1852

1914

! Design matrix valorisation

1853

1915

! Design matrix valorisation

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1854

1916

mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )

1855

1917

! second member valorisation

1856

1918

! second member valorisation

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1857

1919

bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))

1858

1920

! query workspace size

1859

1921

! query workspace size

call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)

1860

1922

call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)

lwork=twork(1)

1861

1923

lwork=twork(1)

allocate(work(int(lwork)))

1862

1924

allocate(work(int(lwork)))

! real call

1863

1925

! real call

call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)

1864

1926

call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)

1865

1927

if (info /= 0) then

1866

1928

if (info /= 0) then

status=0

1867

1929

status=0

end if

1868

1930

end if

1869

1931

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1870

1932

coeff = (/ (bmat(i,1),i=1,deg+1) /)

1871

1933

deallocate(work)

1872

1934

deallocate(work)

deallocate(mat,bmat,singval)

1873

1935

deallocate(mat,bmat,singval)

end subroutine

1874

1936

end subroutine

1875

1937

1876

1938

!

1877

1939

!

! Muller

1878

1940

! Muller

!

1879

1941

!

1880

1942

!

1881

1943

!

! William Daniau 2007

1882

1944

! William Daniau 2007

!

1883

1945

!

! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f

1884

1946

! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f

! http://plato.asu.edu/ftp/other_software/muller.f

1885

1947

! http://plato.asu.edu/ftp/other_software/muller.f

!

1886

1948

!

! it can be used as a replacement for imsl routine dzanly with minor changes

1887

1949

! it can be used as a replacement for imsl routine dzanly with minor changes

!

1888

1950

!

!-----------------------------------------------------------------------

1889

1951

!-----------------------------------------------------------------------

!

1890

1952

!

! purpose - zeros of an analytic complex function

1891

1953

! purpose - zeros of an analytic complex function

! using the muller method with deflation

1892

1954

! using the muller method with deflation

!

1893

1955

!

! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,

1894

1956

! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,

! infer,ier)

1895

1957

! infer,ier)

!

1896

1958

!

! arguments f - a complex function subprogram, f(z), written

1897

1959

! arguments f - a complex function subprogram, f(z), written

! by the user specifying the equation whose

1898

1960

! by the user specifying the equation whose

! roots are to be found. f must appear in

1899

1961

! roots are to be found. f must appear in

! an external statement in the calling pro-

1900

1962

! an external statement in the calling pro-

! gram.

1901

1963

! gram.

! eps - 1st stopping criterion. let fp(z)=f(z)/p

1902

1964

! eps - 1st stopping criterion. let fp(z)=f(z)/p

! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))

1903

1965

! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))

! and z(1),...,z(k-1) are previously found

1904

1966

! and z(1),...,z(k-1) are previously found

! roots. if ((cdabs(f(z)).le.eps) .and.

1905

1967

! roots. if ((cdabs(f(z)).le.eps) .and.

! (cdabs(fp(z)).le.eps)), then z is accepted

1906

1968

! (cdabs(fp(z)).le.eps)), then z is accepted

! as a root. (input)

1907

1969

! as a root. (input)

! eps1 - 2nd stopping criterion. a root is accepted

1908

1970

! eps1 - 2nd stopping criterion. a root is accepted

! if two successive approximations to a given

1909

1971

! if two successive approximations to a given

! root agree within eps1. (input)

1910

1972

! root agree within eps1. (input)

! note. if either or both of the stopping

1911

1973

! note. if either or both of the stopping

! criteria are fulfilled, the root is

1912

1974

! criteria are fulfilled, the root is

! accepted.

1913

1975

! accepted.

! kn - the number of known roots which must be stored

1914

1976

! kn - the number of known roots which must be stored

! in x(1),...,x(kn), prior to entry to muller

1915

1977

! in x(1),...,x(kn), prior to entry to muller

! nguess - the number of initial guesses provided. these

1916

1978

! nguess - the number of initial guesses provided. these

! guesses must be stored in x(kn+1),...,

1917

1979

! guesses must be stored in x(kn+1),...,

! x(kn+nguess). nguess must be set equal

1918

1980

! x(kn+nguess). nguess must be set equal

! to zero if no guesses are provided. (input)

1919

1981

! to zero if no guesses are provided. (input)

! n - the number of new roots to be found by

1920

1982

! n - the number of new roots to be found by

! muller (input)

1921

1983

! muller (input)

! x - a complex vector of length kn+n. x(1),...,

1922

1984

! x - a complex vector of length kn+n. x(1),...,

! x(kn) on input must contain any known

1923

1985

! x(kn) on input must contain any known

! roots. x(kn+1),..., x(kn+n) on input may,

1924

1986

! roots. x(kn+1),..., x(kn+n) on input may,

! on user option, contain initial guesses for

1925

1987

! on user option, contain initial guesses for

! the n new roots which are to be computed.

1926

1988

! the n new roots which are to be computed.

! if the user does not provide an initial

1927

1989

! if the user does not provide an initial

! guess, zero is used.

1928

1990

! guess, zero is used.

! on output, x(kn+1),...,x(kn+n) contain the

1929

1991

! on output, x(kn+1),...,x(kn+n) contain the

! approximate roots found by muller.

1930

1992

! approximate roots found by muller.

! itmax - the maximum allowable number of iterations

1931

1993

! itmax - the maximum allowable number of iterations

! per root (input)

1932

1994

! per root (input)

! infer - an integer vector of length kn+n. on

1933

1995

! infer - an integer vector of length kn+n. on

! output infer(j) contains the number of

1934

1996

! output infer(j) contains the number of

! iterations used in finding the j-th root

1935

1997

! iterations used in finding the j-th root

! when convergence was achieved. if

1936

1998

! when convergence was achieved. if

! convergence was not obtained in itmax

1937

1999

! convergence was not obtained in itmax

! iterations, infer(j) will be greater than

1938

2000

! iterations, infer(j) will be greater than

! itmax (output).

1939

2001

! itmax (output).

! ier - error parameter (output)

1940

2002

! ier - error parameter (output)

! warning error

1941

2003

! warning error

! ier = 33 indicates failure to converge with-

1942

2004

! ier = 33 indicates failure to converge with-

! in itmax iterations for at least one of

1943

2005

! in itmax iterations for at least one of

! the (n) new roots.

1944

2006

! the (n) new roots.

!

1945

2007

!

1946

2008

!

! remarks muller always returns the last approximation for root j

1947

2009

! remarks muller always returns the last approximation for root j

! in x(j). if the convergence criterion is satisfied,

1948

2010

! in x(j). if the convergence criterion is satisfied,

! then infer(j) is less than or equal to itmax. if the

1949

2011

! then infer(j) is less than or equal to itmax. if the

! convergence criterion is not satisified, then infer(j)

1950

2012

! convergence criterion is not satisified, then infer(j)

! is set to either itmax+1 or itmax+k, with k greater

1951

2013

! is set to either itmax+1 or itmax+k, with k greater

! than 1. infer(j) = itmax+1 indicates that muller did

1952

2014

! than 1. infer(j) = itmax+1 indicates that muller did

! not obtain convergence in the allowed number of iter-

1953

2015

! not obtain convergence in the allowed number of iter-

! ations. in this case, the user may wish to set itmax

1954

2016

! ations. in this case, the user may wish to set itmax

! to a larger value. infer(j) = itmax+k means that con-

1955

2017

! to a larger value. infer(j) = itmax+k means that con-

! vergence was obtained (on iteration k) for the defla-

1956

2018

! vergence was obtained (on iteration k) for the defla-

! ted function

1957

2019

! ted function

! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))

1958

2020

! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))

!

1959

2021

!

! but failed for f(z). in this case, better initial

1960

2022

! but failed for f(z). in this case, better initial

! guesses might help or, it might be necessary to relax

1961

2023

! guesses might help or, it might be necessary to relax

! the convergence criterion.

1962

2024

! the convergence criterion.

!

1963

2025

!

!-----------------------------------------------------------------------

1964

2026

!-----------------------------------------------------------------------

!

1965

2027

!

subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)

1966

2028

subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)

implicit none

1967

2029

implicit none

double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq

1968

2030

double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq

double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &

1969

2031

double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &

tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &

1970

2032

tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &

zero,p1,one,four,p5

1971

2033

zero,p1,one,four,p5

1972

2034

double complex, external :: f

1973

2035

double complex, external :: f

integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &

1974

2036

integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &

knpng,jk,ick,nn,lm1,errcode

1975

2037

knpng,jk,ick,nn,lm1,errcode

double complex :: x(kn+n)

1976

2038

double complex :: x(kn+n)

integer :: infer(kn+n)

1977

2039

integer :: infer(kn+n)

1978

2040

1979

2041

data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &

1980

2042

data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &

one/(1.0,0.0)/,four/(4.0,0.0)/, &

1981

2043

one/(1.0,0.0)/,four/(4.0,0.0)/, &

p5/(0.5,0.0)/, &

1982

2044

p5/(0.5,0.0)/, &

rzero/0.0/,rten/10.0/,rhun/100.0/, &

1983

2045

rzero/0.0/,rten/10.0/,rhun/100.0/, &

ax/0.1/,ickmax/3/,rp01/0.01/

1984

2046

ax/0.1/,ickmax/3/,rp01/0.01/

1985

2047

ier = 0

1986

2048

ier = 0

if (n .lt. 1) then ! What the hell are doing here then ...

1987

2049

if (n .lt. 1) then ! What the hell are doing here then ...

return

1988

2050

return

end if

1989

2051

end if

!eps1 = rten **(-nsig)

1990

2052

!eps1 = rten **(-nsig)

eps1 = min(eps1,rp01)

1991

2053

eps1 = min(eps1,rp01)

1992

2054

knp1 = kn+1

1993

2055

knp1 = kn+1

knpn = kn+n

1994

2056

knpn = kn+n

knpng = kn+nguess

1995

2057

knpng = kn+nguess

do i=1,knpn

1996

2058

do i=1,knpn

infer(i) = 0

1997

2059

infer(i) = 0

if (i .gt. knpng) x(i) = zero

1998

2060

if (i .gt. knpng) x(i) = zero

end do

1999

2061

end do

l= knp1

2000

2062

l= knp1

2001

2063

ic=0

2002

2064

ic=0

rloop: do while (l<=knpn) ! Main loop over new roots

2003

2065

rloop: do while (l<=knpn) ! Main loop over new roots

jk = 0

2004

2066

jk = 0

ick = 0

2005

2067

ick = 0

xl = x(l)

2006

2068

xl = x(l)

icloop: do

2007

2069

icloop: do

ic = 0

2008

2070

ic = 0

h = ax

2009

2071

h = ax

h = p1*h

2010

2072

h = p1*h

if (cdabs(xl) .gt. ax) h = p1*xl

2011

2073

if (cdabs(xl) .gt. ax) h = p1*xl

! first three points are

2012

2074

! first three points are

! xl+h, xl-h, xl

2013

2075

! xl+h, xl-h, xl

rt = xl+h

2014

2076

rt = xl+h

call deflated_work(errcode)

2015

2077

call deflated_work(errcode)

if (errcode == 1) then

2016

2078

if (errcode == 1) then

exit icloop

2017

2079

exit icloop

end if

2018

2080

end if

2019

2081

z0 = fprt

2020

2082

z0 = fprt

y0 = frt

2021

2083

y0 = frt

x0 = rt

2022

2084

x0 = rt

rt = xl-h

2023

2085

rt = xl-h

call deflated_work(errcode)

2024

2086

call deflated_work(errcode)

if (errcode == 1) then

2025

2087

if (errcode == 1) then

exit icloop

2026

2088

exit icloop

end if

2027

2089

end if

2028

2090

z1 = fprt

2029

2091

z1 = fprt

y1 = frt

2030

2092

y1 = frt

h = xl-rt

2031

2093

h = xl-rt

d = h/(rt-x0)

2032

2094

d = h/(rt-x0)

rt = xl

2033

2095

rt = xl

2034

2096

call deflated_work(errcode)

2035

2097

call deflated_work(errcode)

if (errcode == 1) then

2036

2098

if (errcode == 1) then

exit icloop

2037

2099

exit icloop

end if

2038

2100

end if

2039

2101

2040

2102

z2 = fprt

2041

2103

z2 = fprt

y2 = frt

2042

2104

y2 = frt

! begin main algorithm

2043

2105

! begin main algorithm

iloop: do

2044

2106

iloop: do

dd = one + d

2045

2107

dd = one + d

t1 = z0*d*d

2046

2108

t1 = z0*d*d

t2 = z1*dd*dd

2047

2109

t2 = z1*dd*dd

xx = z2*dd

2048

2110

xx = z2*dd

t3 = z2*d

2049

2111

t3 = z2*d

bi = t1-t2+xx+t3

2050

2112

bi = t1-t2+xx+t3

den = bi*bi-four*(xx*t1-t3*(t2-xx))

2051

2113

den = bi*bi-four*(xx*t1-t3*(t2-xx))

! use denominator of maximum amplitude

2052

2114

! use denominator of maximum amplitude

t1 = cdsqrt(den)

2053

2115

t1 = cdsqrt(den)

qz = rhun*max(cdabs(bi),cdabs(t1))

2054

2116

qz = rhun*max(cdabs(bi),cdabs(t1))

t2 = bi + t1

2055

2117

t2 = bi + t1

tpq = cdabs(t2)+qz

2056

2118

tpq = cdabs(t2)+qz

if (tpq .eq. qz) t2 = zero

2057

2119

if (tpq .eq. qz) t2 = zero

t3 = bi - t1

2058

2120

t3 = bi - t1

tpq = cdabs(t3) + qz

2059

2121

tpq = cdabs(t3) + qz

if (tpq .eq. qz) t3 = zero

2060

2122

if (tpq .eq. qz) t3 = zero

den = t2

2061

2123

den = t2

qz = cdabs(t3)-cdabs(t2)

2062

2124

qz = cdabs(t3)-cdabs(t2)

if (qz .gt. rzero) den = t3

2063

2125

if (qz .gt. rzero) den = t3

! test for zero denominator

2064

2126

! test for zero denominator

if (cdabs(den) .eq. rzero) then

2065

2127

if (cdabs(den) .eq. rzero) then

call trans_rt()

2066

2128

call trans_rt()

call deflated_work(errcode)

2067

2129

call deflated_work(errcode)

if (errcode == 1) then

2068

2130

if (errcode == 1) then

exit icloop

2069

2131

exit icloop

end if

2070

2132

end if

z2 = fprt

2071

2133

z2 = fprt

y2 = frt

2072

2134

y2 = frt

cycle iloop

2073

2135

cycle iloop

end if

2074

2136

end if

2075

2137

2076

2138

d = -xx/den

2077

2139

d = -xx/den

d = d+d

2078

2140

d = d+d

h = d*h

2079

2141

h = d*h

rt = rt + h

2080

2142

rt = rt + h

! check convergence of the first kind

2081

2143

! check convergence of the first kind

if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then

2082

2144

if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then

if (ic .ne. 0) then

2083

2145

if (ic .ne. 0) then

exit icloop

2084

2146

exit icloop

end if

2085

2147

end if

ic = 1

2086

2148

ic = 1

z0 = y1

2087

2149

z0 = y1

z1 = y2

2088

2150

z1 = y2

z2 = f(rt)

2089

2151

z2 = f(rt)

xl = rt

2090

2152

xl = rt

ick = ick+1

2091

2153

ick = ick+1

if (ick .le. ickmax) then

2092

2154

if (ick .le. ickmax) then

cycle iloop

2093

2155

cycle iloop

end if

2094

2156

end if

! warning error, itmax = maximum

2095

2157

! warning error, itmax = maximum

jk = itmax + jk

2096

2158

jk = itmax + jk

ier = 33

2097

2159

ier = 33

end if

2098

2160

end if

if (ic .ne. 0) then

2099

2161

if (ic .ne. 0) then

cycle icloop

2100

2162

cycle icloop

end if

2101

2163

end if

call deflated_work(errcode)

2102

2164

call deflated_work(errcode)

if (errcode == 1) then

2103

2165

if (errcode == 1) then

exit icloop

2104

2166

exit icloop

end if

2105

2167

end if

2106

2168

do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)

2107

2169

do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)

! take remedial action to induce

2108

2170

! take remedial action to induce

! convergence

2109

2171

! convergence

d = d*p5

2110

2172

d = d*p5

h = h*p5

2111

2173

h = h*p5

rt = rt-h

2112

2174

rt = rt-h

call deflated_work(errcode)

2113

2175

call deflated_work(errcode)

if (errcode == 1) then

2114

2176

if (errcode == 1) then

exit icloop

2115

2177

exit icloop

end if

2116

2178

end if

end do

2117

2179

end do

z0 = z1

2118

2180

z0 = z1

z1 = z2

2119

2181

z1 = z2

z2 = fprt

2120

2182

z2 = fprt

y0 = y1

2121

2183

y0 = y1

y1 = y2

2122

2184

y1 = y2

y2 = frt

2123

2185

y2 = frt

end do iloop

2124

2186

end do iloop

end do icloop

2125

2187

end do icloop

x(l) = rt

2126

2188

x(l) = rt

infer(l) = jk

2127

2189

infer(l) = jk

l = l+1

2128

2190

l = l+1

end do rloop

2129

2191

end do rloop

2130

2192

contains

2131

2193

contains

subroutine trans_rt()

2132

2194

subroutine trans_rt()

tem = rten*eps1

2133

2195

tem = rten*eps1

if (cdabs(rt) .gt. ax) tem = tem*rt

2134

2196

if (cdabs(rt) .gt. ax) tem = tem*rt

rt = rt+tem

2135

2197

rt = rt+tem

d = (h+tem)*d/h

2136

2198

d = (h+tem)*d/h

h = h+tem

2137

2199

h = h+tem

end subroutine trans_rt

2138

2200

end subroutine trans_rt

2139

2201

subroutine deflated_work(errcode)

2140

2202

subroutine deflated_work(errcode)

! errcode=0 => no errors

2141

2203

! errcode=0 => no errors

! errcode=1 => jk>itmax or convergence of second kind achieved

2142

2204

! errcode=1 => jk>itmax or convergence of second kind achieved

integer :: errcode,flag

2143

2205

integer :: errcode,flag

2144

2206

flag=1

2145

2207

flag=1

loop1: do while(flag==1)

2146

2208

loop1: do while(flag==1)

errcode=0

2147

2209

errcode=0

jk = jk+1

2148

2210

jk = jk+1

if (jk .gt. itmax) then

2149

2211

if (jk .gt. itmax) then

ier=33

2150

2212

ier=33

errcode=1

2151

2213

errcode=1

return

2152

2214

return

end if

2153

2215

end if

frt = f(rt)

2154

2216

frt = f(rt)

fprt = frt

2155

2217

fprt = frt

if (l /= 1) then

2156

2218

if (l /= 1) then

lm1 = l-1

2157

2219

lm1 = l-1

do i=1,lm1

2158

2220

do i=1,lm1

tem = rt - x(i)

2159

2221

tem = rt - x(i)

if (cdabs(tem) .eq. rzero) then

2160

2222

if (cdabs(tem) .eq. rzero) then

!if (ic .ne. 0) go to 15 !! ?? possible?

2161

2223

!if (ic .ne. 0) go to 15 !! ?? possible?

call trans_rt()

2162

2224

call trans_rt()

cycle loop1

2163

2225

cycle loop1

end if

2164

2226

end if

fprt = fprt/tem

2165

2227

fprt = fprt/tem

end do

2166

2228

end do

end if

2167

2229

end if

flag=0

2168

2230

flag=0

end do loop1

2169

2231

end do loop1

2170

2232

if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then

2171

2233

if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then

errcode=1

2172

2234

errcode=1

return

2173

2235

return

end if

2174

2236

end if

2175

2237

end subroutine deflated_work

2176

2238

end subroutine deflated_work

2177

2239

end subroutine

2178

2240

end subroutine

2179

2241

2180

2242

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

2181

2243

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

!

2182

2244

!

! Integration

2183

2245

! Integration

!

2184

2246

!

! Only double precision coded atm

2185

2247

! Only double precision coded atm

!

2186

2248

!

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

2187

2249

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

2188

2250

2189

2251

subroutine fvn_d_gauss_legendre(n,qx,qw)

2190

2252

subroutine fvn_d_gauss_legendre(n,qx,qw)

!

2191

2253

!

! This routine compute the n Gauss Legendre abscissas and weights

2192

2254

! This routine compute the n Gauss Legendre abscissas and weights

! Adapted from Numerical Recipes routine gauleg

2193

2255

! Adapted from Numerical Recipes routine gauleg

!

2194

2256

!

! n (in) : number of points

2195

2257

! n (in) : number of points

! qx(out) : abscissas

2196

2258

! qx(out) : abscissas

! qw(out) : weights

2197

2259

! qw(out) : weights

!

2198

2260

!

implicit none

2199

2261

implicit none

double precision,parameter :: pi=3.141592653589793d0

2200

2262

double precision,parameter :: pi=3.141592653589793d0

integer, intent(in) :: n

2201

2263

integer, intent(in) :: n

double precision, intent(out) :: qx(n),qw(n)

2202

2264

double precision, intent(out) :: qx(n),qw(n)

2203

2265

integer :: m,i,j

2204

2266

integer :: m,i,j

double precision :: z,z1,p1,p2,p3,pp

2205

2267

double precision :: z,z1,p1,p2,p3,pp

m=(n+1)/2

2206

2268

m=(n+1)/2

do i=1,m

2207

2269

do i=1,m

z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))

2208

2270

z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))

iloop: do

2209

2271

iloop: do

p1=1.d0

2210

2272

p1=1.d0

p2=0.d0

2211

2273

p2=0.d0

do j=1,n

2212

2274

do j=1,n

p3=p2

2213

2275

p3=p2

p2=p1

2214

2276

p2=p1

p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)

2215

2277

p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)

end do

2216

2278

end do

pp=dble(n)*(z*p1-p2)/(z*z-1.d0)

2217

2279

pp=dble(n)*(z*p1-p2)/(z*z-1.d0)

z1=z

2218

2280

z1=z

z=z1-p1/pp

2219

2281

z=z1-p1/pp

if (dabs(z-z1)<=epsilon(z)) then

2220

2282

if (dabs(z-z1)<=epsilon(z)) then

exit iloop

2221

2283

exit iloop

end if

2222

2284

end if

end do iloop

2223

2285

end do iloop

qx(i)=-z

2224

2286

qx(i)=-z

qx(n+1-i)=z

2225

2287

qx(n+1-i)=z

qw(i)=2.d0/((1.d0-z*z)*pp*pp)

2226

2288

qw(i)=2.d0/((1.d0-z*z)*pp*pp)

qw(n+1-i)=qw(i)

2227

2289

qw(n+1-i)=qw(i)

end do

2228

2290

end do

end subroutine

2229

2291

end subroutine

2230

2292

2231

2293

2232

2294

subroutine fvn_d_gl_integ(f,a,b,n,res)

2233

2295

subroutine fvn_d_gl_integ(f,a,b,n,res)

!

2234

2296

!

! This is a simple non adaptative integration routine

2235

2297

! This is a simple non adaptative integration routine

! using n gauss legendre abscissas and weights

2236

2298

! using n gauss legendre abscissas and weights

!

2237

2299

!

! f(in) : the function to integrate

2238

2300

! f(in) : the function to integrate

! a(in) : lower bound

2239

2301

! a(in) : lower bound

! b(in) : higher bound

2240

2302

! b(in) : higher bound

! n(in) : number of gauss legendre pairs

2241

2303

! n(in) : number of gauss legendre pairs

! res(out): the evaluation of the integral

2242

2304

! res(out): the evaluation of the integral

!

2243

2305

!

double precision,external :: f

2244

2306

double precision,external :: f

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

2245

2307

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

integer, intent(in):: n

2246

2308

integer, intent(in):: n

double precision, intent(out) :: res

2247

2309

double precision, intent(out) :: res

2248

2310

double precision, allocatable :: qx(:),qw(:)

2249

2311

double precision, allocatable :: qx(:),qw(:)

double precision :: xm,xr

2250

2312

double precision :: xm,xr

integer :: i

2251

2313

integer :: i

2252

2314

! First compute n gauss legendre abs and weight

2253

2315

! First compute n gauss legendre abs and weight

allocate(qx(n))

2254

2316

allocate(qx(n))

allocate(qw(n))

2255

2317

allocate(qw(n))

call fvn_d_gauss_legendre(n,qx,qw)

2256

2318

call fvn_d_gauss_legendre(n,qx,qw)

2257

2319

xm=0.5d0*(b+a)

2258

2320

xm=0.5d0*(b+a)

xr=0.5d0*(b-a)

2259

2321

xr=0.5d0*(b-a)

2260

2322

res=0.d0

2261

2323

res=0.d0

2262

2324

do i=1,n

2263

2325

do i=1,n

res=res+qw(i)*f(xm+xr*qx(i))

2264

2326

res=res+qw(i)*f(xm+xr*qx(i))

end do

2265

2327

end do

2266

2328

res=xr*res

2267

2329

res=xr*res

2268

2330

deallocate(qw)

2269

2331

deallocate(qw)

deallocate(qx)

2270

2332

deallocate(qx)

2271

2333

end subroutine

2272

2334

end subroutine

2273

2335

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

2274

2336

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

!

2275

2337

!

! Simple and double adaptative Gauss Kronrod integration based on

2276

2338

! Simple and double adaptative Gauss Kronrod integration based on

! a modified version of quadpack ( http://www.netlib.org/quadpack

2277

2339

! a modified version of quadpack ( http://www.netlib.org/quadpack

!

2278

2340

!

! Common parameters :

2279

2341

! Common parameters :

!

2280

2342

!

! key (in)

2281

2343

! key (in)

! epsabs

2282

2344

! epsabs

! epsrel

2283

2345

! epsrel

!

2284

2346

!

2285

2347

!

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

2286

2348

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

2287

2349

subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

2288

2350

subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

!

2289

2351

!

! Evaluate the integral of function f(x) between a and b

2290

2352

! Evaluate the integral of function f(x) between a and b

!

2291

2353

!

! f(in) : the function

2292

2354

! f(in) : the function

! a(in) : lower bound

2293

2355

! a(in) : lower bound

! b(in) : higher bound

2294

2356

! b(in) : higher bound

! epsabs(in) : desired absolute error

2295

2357

! epsabs(in) : desired absolute error

! epsrel(in) : desired relative error

2296

2358

! epsrel(in) : desired relative error

! key(in) : gauss kronrod rule

2297

2359

! key(in) : gauss kronrod rule

! 1: 7 - 15 points

2298

2360

! 1: 7 - 15 points

! 2: 10 - 21 points

2299

2361

! 2: 10 - 21 points

! 3: 15 - 31 points

2300

2362

! 3: 15 - 31 points

! 4: 20 - 41 points

2301

2363

! 4: 20 - 41 points

! 5: 25 - 51 points

2302

2364

! 5: 25 - 51 points

! 6: 30 - 61 points

2303

2365

! 6: 30 - 61 points

!

2304

2366

!

! limit(in) : maximum number of subintervals in the partition of the

2305

2367

! limit(in) : maximum number of subintervals in the partition of the

! given integration interval (a,b). A value of 500 will give the same

2306

2368

! given integration interval (a,b). A value of 500 will give the same

! behaviour as the imsl routine dqdag

2307

2369

! behaviour as the imsl routine dqdag

!

2308

2370

!

! res(out) : estimated integral value

2309

2371

! res(out) : estimated integral value

! abserr(out) : estimated absolute error

2310

2372

! abserr(out) : estimated absolute error

! ier(out) : error flag from quadpack routines

2311

2373

! ier(out) : error flag from quadpack routines

! 0 : no error

2312

2374

! 0 : no error

! 1 : maximum number of subdivisions allowed

2313

2375

! 1 : maximum number of subdivisions allowed

! has been achieved. one can allow more

2314

2376

! has been achieved. one can allow more

! subdivisions by increasing the value of

2315

2377

! subdivisions by increasing the value of

! limit (and taking the according dimension

2316

2378

! limit (and taking the according dimension

! adjustments into account). however, if

2317

2379

! adjustments into account). however, if

! this yield no improvement it is advised

2318

2380

! this yield no improvement it is advised

! to analyze the integrand in order to

2319

2381

! to analyze the integrand in order to

! determine the integration difficulaties.

2320

2382

! determine the integration difficulaties.

! if the position of a local difficulty can

2321

2383

! if the position of a local difficulty can

! be determined (i.e.singularity,

2322

2384

! be determined (i.e.singularity,

! discontinuity within the interval) one

2323

2385

! discontinuity within the interval) one

! will probably gain from splitting up the

2324

2386

! will probably gain from splitting up the

! interval at this point and calling the

2325

2387

! interval at this point and calling the

! integrator on the subranges. if possible,

2326

2388

! integrator on the subranges. if possible,

! an appropriate special-purpose integrator

2327

2389

! an appropriate special-purpose integrator

! should be used which is designed for

2328

2390

! should be used which is designed for

! handling the type of difficulty involved.

2329

2391

! handling the type of difficulty involved.

! 2 : the occurrence of roundoff error is

2330

2392

! 2 : the occurrence of roundoff error is

! detected, which prevents the requested

2331

2393

! detected, which prevents the requested

! tolerance from being achieved.

2332

2394

! tolerance from being achieved.

! 3 : extremely bad integrand behaviour occurs

2333

2395

! 3 : extremely bad integrand behaviour occurs

! at some points of the integration

2334

2396

! at some points of the integration

! interval.

2335

2397

! interval.

! 6 : the input is invalid, because

2336

2398

! 6 : the input is invalid, because

! (epsabs.le.0 and

2337

2399

! (epsabs.le.0 and

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

2338

2400

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

! or limit.lt.1 or lenw.lt.limit*4.

2339

2401

! or limit.lt.1 or lenw.lt.limit*4.

! result, abserr, neval, last are set

2340

2402

! result, abserr, neval, last are set

! to zero.

2341

2403

! to zero.

! except when lenw is invalid, iwork(1),

2342

2404

! except when lenw is invalid, iwork(1),

! work(limit*2+1) and work(limit*3+1) are

2343

2405

! work(limit*2+1) and work(limit*3+1) are

! set to zero, work(1) is set to a and

2344

2406

! set to zero, work(1) is set to a and

! work(limit+1) to b.

2345

2407

! work(limit+1) to b.

2346

2408

implicit none

2347

2409

implicit none

double precision, external :: f

2348

2410

double precision, external :: f

double precision, intent(in) :: a,b,epsabs,epsrel

2349

2411

double precision, intent(in) :: a,b,epsabs,epsrel

integer, intent(in) :: key

2350

2412

integer, intent(in) :: key

integer, intent(in) :: limit

2351

2413

integer, intent(in) :: limit

double precision, intent(out) :: res,abserr

2352

2414

double precision, intent(out) :: res,abserr

integer, intent(out) :: ier

2353

2415

integer, intent(out) :: ier

2354

2416

double precision, allocatable :: work(:)

2355

2417

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

2356

2418

integer, allocatable :: iwork(:)

integer :: lenw,neval,last

2357

2419

integer :: lenw,neval,last

2358

2420

! imsl value for limit is 500

2359

2421

! imsl value for limit is 500

lenw=limit*4

2360

2422

lenw=limit*4

2361

2423

allocate(iwork(limit))

2362

2424

allocate(iwork(limit))

allocate(work(lenw))

2363

2425

allocate(work(lenw))

2364

2426

call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

2365

2427

call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

2366

2428

deallocate(work)

2367

2429

deallocate(work)

deallocate(iwork)

2368

2430

deallocate(iwork)

2369

2431

end subroutine

2370

2432

end subroutine

2371

2433

2372

2434

2373

2435

subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)

2374

2436

subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)

!

2375

2437

!

! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)

2376

2438

! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)

!

2377

2439

!

! f(in) : the function

2378

2440

! f(in) : the function

! a(in) : lower bound

2379

2441

! a(in) : lower bound

! b(in) : higher bound

2380

2442

! b(in) : higher bound

! g(in) : external function describing lower bound for y

2381

2443

! g(in) : external function describing lower bound for y

! h(in) : external function describing higher bound for y

2382

2444

! h(in) : external function describing higher bound for y

! epsabs(in) : desired absolute error

2383

2445

! epsabs(in) : desired absolute error

! epsrel(in) : desired relative error

2384

2446

! epsrel(in) : desired relative error

! key(in) : gauss kronrod rule

2385

2447

! key(in) : gauss kronrod rule

! 1: 7 - 15 points

2386

2448

! 1: 7 - 15 points

! 2: 10 - 21 points

2387

2449

! 2: 10 - 21 points

! 3: 15 - 31 points

2388

2450

! 3: 15 - 31 points

! 4: 20 - 41 points

2389

2451

! 4: 20 - 41 points

! 5: 25 - 51 points

2390

2452

! 5: 25 - 51 points

! 6: 30 - 61 points

2391

2453

! 6: 30 - 61 points

!

2392

2454

!

! limit(in) : maximum number of subintervals in the partition of the

2393

2455

! limit(in) : maximum number of subintervals in the partition of the

! given integration interval (a,b). A value of 500 will give the same

2394

2456

! given integration interval (a,b). A value of 500 will give the same

! behaviour as the imsl routine dqdag

2395

2457

! behaviour as the imsl routine dqdag

!

2396

2458

!

! res(out) : estimated integral value

2397

2459

! res(out) : estimated integral value

! abserr(out) : estimated absolute error

2398

2460

! abserr(out) : estimated absolute error

! ier(out) : error flag from quadpack routines

2399

2461

! ier(out) : error flag from quadpack routines

! 0 : no error

2400

2462

! 0 : no error

! 1 : maximum number of subdivisions allowed

2401

2463

! 1 : maximum number of subdivisions allowed

! has been achieved. one can allow more

2402

2464

! has been achieved. one can allow more

! subdivisions by increasing the value of

2403

2465

! subdivisions by increasing the value of

! limit (and taking the according dimension

2404

2466

! limit (and taking the according dimension

! adjustments into account). however, if

2405

2467

! adjustments into account). however, if

! this yield no improvement it is advised

2406

2468

! this yield no improvement it is advised

! to analyze the integrand in order to

2407

2469

! to analyze the integrand in order to

! determine the integration difficulaties.

2408

2470

! determine the integration difficulaties.

! if the position of a local difficulty can

2409

2471

! if the position of a local difficulty can

! be determined (i.e.singularity,

2410

2472

! be determined (i.e.singularity,

! discontinuity within the interval) one

2411

2473

! discontinuity within the interval) one

! will probably gain from splitting up the

2412

2474

! will probably gain from splitting up the

! interval at this point and calling the

2413

2475

! interval at this point and calling the

! integrator on the subranges. if possible,

2414

2476

! integrator on the subranges. if possible,

! an appropriate special-purpose integrator

2415

2477

! an appropriate special-purpose integrator

! should be used which is designed for

2416

2478

! should be used which is designed for

! handling the type of difficulty involved.

2417

2479

! handling the type of difficulty involved.

! 2 : the occurrence of roundoff error is

2418

2480

! 2 : the occurrence of roundoff error is

! detected, which prevents the requested

2419

2481

! detected, which prevents the requested

! tolerance from being achieved.

2420

2482

! tolerance from being achieved.

! 3 : extremely bad integrand behaviour occurs

2421

2483

! 3 : extremely bad integrand behaviour occurs

! at some points of the integration

2422

2484

! at some points of the integration

! interval.

2423

2485

! interval.

! 6 : the input is invalid, because

2424

2486

! 6 : the input is invalid, because

! (epsabs.le.0 and

2425

2487

! (epsabs.le.0 and

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

2426

2488

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

! or limit.lt.1 or lenw.lt.limit*4.

2427

2489

! or limit.lt.1 or lenw.lt.limit*4.

! result, abserr, neval, last are set

2428

2490

! result, abserr, neval, last are set

! to zero.

2429

2491

! to zero.

! except when lenw is invalid, iwork(1),

2430

2492

! except when lenw is invalid, iwork(1),

! work(limit*2+1) and work(limit*3+1) are

2431

2493

! work(limit*2+1) and work(limit*3+1) are

! set to zero, work(1) is set to a and

2432

2494

! set to zero, work(1) is set to a and

! work(limit+1) to b.

2433

2495

! work(limit+1) to b.

2434

2496

implicit none

2435

2497

implicit none

double precision, external:: f,g,h

2436

2498

double precision, external:: f,g,h

double precision, intent(in) :: a,b,epsabs,epsrel

2437

2499

double precision, intent(in) :: a,b,epsabs,epsrel

integer, intent(in) :: key,limit

2438

2500

integer, intent(in) :: key,limit

integer, intent(out) :: ier

2439

2501

integer, intent(out) :: ier

double precision, intent(out) :: res,abserr

2440

2502

double precision, intent(out) :: res,abserr

2441

2503

2442

2504

double precision, allocatable :: work(:)

2443

2505

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

2444

2506

integer, allocatable :: iwork(:)

integer :: lenw,neval,last

2445

2507

integer :: lenw,neval,last

2446

2508

! imsl value for limit is 500

2447

2509

! imsl value for limit is 500

lenw=limit*4

2448

2510

lenw=limit*4

allocate(work(lenw))

2449

2511

allocate(work(lenw))

allocate(iwork(limit))

2450

2512

allocate(iwork(limit))

2451

2513

call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

2452

2514

call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

2453

2515

deallocate(iwork)

2454

2516

deallocate(iwork)

deallocate(work)

2455

2517

deallocate(work)

end subroutine

2456

2518

end subroutine

2457

2519

2458

2520

2459

2521

subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

2460

2522

subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

!

2461

2523

!

! Evaluate the single integral of function f(x,y) for y between a and b with a

2462

2524

! Evaluate the single integral of function f(x,y) for y between a and b with a

! given x value

2463

2525

! given x value

!

2464

2526

!

! This function is used for the evaluation of the double integral fvn_d_integ_2_gk

2465

2527

! This function is used for the evaluation of the double integral fvn_d_integ_2_gk

!

2466

2528

!

! f(in) : the function

2467

2529

! f(in) : the function

! x(in) : x

2468

2530

! x(in) : x

! a(in) : lower bound

2469

2531

! a(in) : lower bound

! b(in) : higher bound

2470

2532

! b(in) : higher bound

! epsabs(in) : desired absolute error

2471

2533

! epsabs(in) : desired absolute error

! epsrel(in) : desired relative error

2472

2534

! epsrel(in) : desired relative error

! key(in) : gauss kronrod rule

2473

2535

! key(in) : gauss kronrod rule

! 1: 7 - 15 points

2474

2536

! 1: 7 - 15 points

! 2: 10 - 21 points

2475

2537

! 2: 10 - 21 points

! 3: 15 - 31 points

2476

2538

! 3: 15 - 31 points

! 4: 20 - 41 points

2477

2539

! 4: 20 - 41 points

! 5: 25 - 51 points

2478

2540

! 5: 25 - 51 points

! 6: 30 - 61 points

2479

2541

! 6: 30 - 61 points

!

2480

2542

!

! limit(in) : maximum number of subintervals in the partition of the

2481

2543

! limit(in) : maximum number of subintervals in the partition of the

! given integration interval (a,b). A value of 500 will give the same

2482

2544

! given integration interval (a,b). A value of 500 will give the same

! behaviour as the imsl routine dqdag

2483

2545

! behaviour as the imsl routine dqdag

!

2484

2546

!

! res(out) : estimated integral value

2485

2547

! res(out) : estimated integral value

! abserr(out) : estimated absolute error

2486

2548

! abserr(out) : estimated absolute error

! ier(out) : error flag from quadpack routines

2487

2549

! ier(out) : error flag from quadpack routines

! 0 : no error

2488

2550

! 0 : no error

! 1 : maximum number of subdivisions allowed

2489

2551

! 1 : maximum number of subdivisions allowed

! has been achieved. one can allow more

2490

2552

! has been achieved. one can allow more

! subdivisions by increasing the value of

2491

2553

! subdivisions by increasing the value of

! limit (and taking the according dimension

2492

2554

! limit (and taking the according dimension

! adjustments into account). however, if

2493

2555

! adjustments into account). however, if

! this yield no improvement it is advised

2494

2556

! this yield no improvement it is advised

! to analyze the integrand in order to

2495

2557

! to analyze the integrand in order to

! determine the integration difficulaties.

2496

2558

! determine the integration difficulaties.

! if the position of a local difficulty can

2497

2559

! if the position of a local difficulty can

! be determined (i.e.singularity,

2498

2560

! be determined (i.e.singularity,

! discontinuity within the interval) one

2499

2561

! discontinuity within the interval) one

! will probably gain from splitting up the

2500

2562

! will probably gain from splitting up the

! interval at this point and calling the

2501

2563

! interval at this point and calling the

! integrator on the subranges. if possible,

2502

2564

! integrator on the subranges. if possible,

! an appropriate special-purpose integrator

2503

2565

! an appropriate special-purpose integrator

! should be used which is designed for

2504

2566

! should be used which is designed for

! handling the type of difficulty involved.

2505

2567

! handling the type of difficulty involved.

! 2 : the occurrence of roundoff error is

2506

2568

! 2 : the occurrence of roundoff error is

! detected, which prevents the requested

2507

2569

! detected, which prevents the requested

! tolerance from being achieved.

2508

2570

! tolerance from being achieved.

! 3 : extremely bad integrand behaviour occurs

2509

2571

! 3 : extremely bad integrand behaviour occurs

! at some points of the integration

2510

2572

! at some points of the integration

! interval.

2511

2573

! interval.

! 6 : the input is invalid, because

2512

2574

! 6 : the input is invalid, because

! (epsabs.le.0 and

2513

2575

! (epsabs.le.0 and

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

2514

2576

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

! or limit.lt.1 or lenw.lt.limit*4.

2515

2577

! or limit.lt.1 or lenw.lt.limit*4.

! result, abserr, neval, last are set

2516

2578

! result, abserr, neval, last are set

! to zero.

2517

2579

! to zero.

! except when lenw is invalid, iwork(1),

2518

2580

! except when lenw is invalid, iwork(1),

! work(limit*2+1) and work(limit*3+1) are

2519

2581

! work(limit*2+1) and work(limit*3+1) are

! set to zero, work(1) is set to a and

2520

2582

! set to zero, work(1) is set to a and

! work(limit+1) to b.

2521

2583

! work(limit+1) to b.

2522

2584

implicit none

2523

2585

implicit none

double precision, external:: f

2524

2586

double precision, external:: f

double precision, intent(in) :: x,a,b,epsabs,epsrel

2525

2587

double precision, intent(in) :: x,a,b,epsabs,epsrel

integer, intent(in) :: key,limit

2526

2588

integer, intent(in) :: key,limit

integer, intent(out) :: ier

2527

2589

integer, intent(out) :: ier

double precision, intent(out) :: res,abserr

2528

2590

double precision, intent(out) :: res,abserr

2529

2591

2530

2592

double precision, allocatable :: work(:)

2531

2593

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

2532

2594

integer, allocatable :: iwork(:)

integer :: lenw,neval,last

2533

2595

integer :: lenw,neval,last

2534

2596

! imsl value for limit is 500

2535

2597

! imsl value for limit is 500

lenw=limit*4

2536

2598

lenw=limit*4

allocate(work(lenw))

2537

2599

allocate(work(lenw))

allocate(iwork(limit))

2538

2600

allocate(iwork(limit))

2539

2601

call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

2540

2602

call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

2541

2603

deallocate(iwork)

2542

2604

deallocate(iwork)

deallocate(work)

2543

2605

deallocate(work)

end subroutine

2544

2606

end subroutine

2545

2607

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

2546

2608

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

! Include the modified quadpack files

2547

2609

! Include the modified quadpack files

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

2548

2610

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

include "fvn_quadpack/dqag_2d_inner.f"

2549

2611

include "fvn_quadpack/dqag_2d_inner.f"

include "fvn_quadpack/dqk15_2d_inner.f"

2550

2612

include "fvn_quadpack/dqk15_2d_inner.f"

include "fvn_quadpack/dqk31_2d_outer.f"

2551

2613

include "fvn_quadpack/dqk31_2d_outer.f"

include "fvn_quadpack/d1mach.f"

2552

2614

include "fvn_quadpack/d1mach.f"

include "fvn_quadpack/dqk31_2d_inner.f"

2553

2615

include "fvn_quadpack/dqk31_2d_inner.f"

include "fvn_quadpack/dqage.f"

2554

2616

include "fvn_quadpack/dqage.f"

include "fvn_quadpack/dqk15.f"

2555

2617

include "fvn_quadpack/dqk15.f"

include "fvn_quadpack/dqk21.f"

2556

2618

include "fvn_quadpack/dqk21.f"

include "fvn_quadpack/dqk31.f"

2557

2619

include "fvn_quadpack/dqk31.f"

include "fvn_quadpack/dqk41.f"

2558

2620

include "fvn_quadpack/dqk41.f"

include "fvn_quadpack/dqk51.f"

2559

2621

include "fvn_quadpack/dqk51.f"

include "fvn_quadpack/dqk61.f"

2560

2622

include "fvn_quadpack/dqk61.f"

include "fvn_quadpack/dqk41_2d_outer.f"

2561

2623

include "fvn_quadpack/dqk41_2d_outer.f"

include "fvn_quadpack/dqk41_2d_inner.f"

2562

2624

include "fvn_quadpack/dqk41_2d_inner.f"

include "fvn_quadpack/dqag_2d_outer.f"

2563

2625

include "fvn_quadpack/dqag_2d_outer.f"

include "fvn_quadpack/dqpsrt.f"

2564

2626

include "fvn_quadpack/dqpsrt.f"

include "fvn_quadpack/dqag.f"

2565

2627

include "fvn_quadpack/dqag.f"

include "fvn_quadpack/dqage_2d_outer.f"

2566

2628

include "fvn_quadpack/dqage_2d_outer.f"

include "fvn_quadpack/dqage_2d_inner.f"

2567

2629

include "fvn_quadpack/dqage_2d_inner.f"

include "fvn_quadpack/dqk51_2d_outer.f"

2568

2630

include "fvn_quadpack/dqk51_2d_outer.f"

include "fvn_quadpack/dqk51_2d_inner.f"

2569

2631

include "fvn_quadpack/dqk51_2d_inner.f"

include "fvn_quadpack/dqk61_2d_outer.f"

2570

2632

include "fvn_quadpack/dqk61_2d_outer.f"

include "fvn_quadpack/dqk21_2d_outer.f"

2571

2633

include "fvn_quadpack/dqk21_2d_outer.f"

include "fvn_quadpack/dqk61_2d_inner.f"

2572

2634

include "fvn_quadpack/dqk61_2d_inner.f"

include "fvn_quadpack/dqk21_2d_inner.f"

2573

2635

include "fvn_quadpack/dqk21_2d_inner.f"

include "fvn_quadpack/dqk15_2d_outer.f"

2574

2636

include "fvn_quadpack/dqk15_2d_outer.f"

2575

2637

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

2576

2638

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

!

2577

2639

!

! Trigonometric functions

2578

2640

! Trigonometric functions

!

2579

2641

!

! fvn_z_acos, fvn_z_asin : complex arc cosine and sine

2580

2642

! fvn_z_acos, fvn_z_asin : complex arc cosine and sine

! fvn_d_acosh : arc cosinus hyperbolic

2581

2643

! fvn_d_acosh : arc cosinus hyperbolic

! fvn_d_asinh : arc sinus hyperbolic

2582

2644

! fvn_d_asinh : arc sinus hyperbolic

!

2583

2645

!

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

2584

2646

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

function fvn_z_acos(z)

2585

2647

function fvn_z_acos(z)

! double complex arccos function adapted from

2586

2648

! double complex arccos function adapted from

! the c gsl library

2587

2649

! the c gsl library

! http://www.gnu.org/software/gsl/

2588

2650

! http://www.gnu.org/software/gsl/

implicit none

2589

2651

implicit none

complex(kind=8) :: fvn_z_acos

2590

2652

complex(kind=8) :: fvn_z_acos

complex(kind=8) :: z

2591

2653

complex(kind=8) :: z

real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1

2592

2654

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

2593

2655

real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8

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

2594

2656

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

real(kind=8) :: r_res,i_res

2595

2657

real(kind=8) :: r_res,i_res

2596

2658

rz=dreal(z)

2597

2659

rz=dreal(z)

iz=dimag(z)

2598

2660

iz=dimag(z)

if ( iz == 0._8 ) then

2599

2661

if ( iz == 0._8 ) then

fvn_z_acos=fvn_z_acos_real(rz)

2600

2662

fvn_z_acos=fvn_z_acos_real(rz)

return

2601

2663

return

end if

2602

2664

end if

2603

2665

x=dabs(rz)

2604

2666

x=dabs(rz)

y=dabs(iz)

2605

2667

y=dabs(iz)

r=fvn_d_hypot(x+1.,y)

2606

2668

r=fvn_d_hypot(x+1.,y)

s=fvn_d_hypot(x-1.,y)

2607

2669

s=fvn_d_hypot(x-1.,y)

a=0.5*(r + s)

2608

2670

a=0.5*(r + s)

b=x/a

2609

2671

b=x/a

y2=y*y

2610

2672

y2=y*y

2611

2673

if (b <= b_crossover) then

2612

2674

if (b <= b_crossover) then

r_res=dacos(b)

2613

2675

r_res=dacos(b)

else

2614

2676

else

if (x <= 1.) then

2615

2677

if (x <= 1.) then

d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))

2616

2678

d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))

r_res=datan(dsqrt(d)/x)

2617

2679

r_res=datan(dsqrt(d)/x)

else

2618

2680

else

apx=a+x

2619

2681

apx=a+x

d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))

2620

2682

d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))

r_res=datan((y*dsqrt(d))/x);

2621

2683

r_res=datan((y*dsqrt(d))/x);

end if

2622

2684

end if

2623

2685

end if

2624

2686

if (a <= a_crossover) then

2625

2687

if (a <= a_crossover) then

if (x < 1.) then

2626

2688

if (x < 1.) then

am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))

2627

2689

am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))

else

2628

2690

else

am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))

2629

2691

am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))

end if

2630

2692

end if

i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));

2631

2693

i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));

else

2632

2694

else

i_res = dlog (a + dsqrt (a*a - 1.));

2633

2695

i_res = dlog (a + dsqrt (a*a - 1.));

end if

2634

2696

end if

if (rz <0.) then

2635

2697

if (rz <0.) then

r_res=fvn_pi-r_res

2636

2698

r_res=fvn_pi-r_res

end if

2637

2699

end if

i_res=-sign(1._8,iz)*i_res

2638

2700

i_res=-sign(1._8,iz)*i_res

fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res)

2639

2701

fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res)

2640

2702

end function fvn_z_acos

2641

2703

end function fvn_z_acos

2642

2704

function fvn_z_acos_real(r)

2643

2705

function fvn_z_acos_real(r)

! return the double complex arc cosinus for a

2644

2706

! return the double complex arc cosinus for a

! double precision argument

2645

2707

! double precision argument

implicit none

2646

2708

implicit none

real(kind=8) :: r

2647

2709

real(kind=8) :: r

complex(kind=8) :: fvn_z_acos_real

2648

2710

complex(kind=8) :: fvn_z_acos_real

2649

2711

if (dabs(r)<=1._8) then

2650

2712

if (dabs(r)<=1._8) then

fvn_z_acos_real=dcmplx(dacos(r))

2651

2713

fvn_z_acos_real=dcmplx(dacos(r))

return

2652

2714

return

end if

2653

2715

end if

if (r < 0._8) then

2654

2716

if (r < 0._8) then

fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r))

2655

2717

fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r))

else

2656

2718

else

fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r))

2657

2719

fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r))

end if

2658

2720

end if

end function

2659

2721

end function

2660

2722

2661

2723

function fvn_z_asin(z)

2662

2724

function fvn_z_asin(z)

! double complex arcsin function derived from

2663

2725

! double complex arcsin function derived from

! the c gsl library

2664

2726

! the c gsl library

! http://www.gnu.org/software/gsl/

2665

2727

! http://www.gnu.org/software/gsl/

implicit none

2666

2728

implicit none

complex(kind=8) :: fvn_z_asin

2667

2729

complex(kind=8) :: fvn_z_asin

complex(kind=8) :: z

2668

2730

complex(kind=8) :: z

real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1

2669

2731

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

2670

2732

real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8

real(kind=8) :: r_res,i_res

2671

2733

real(kind=8) :: r_res,i_res

2672

2734

rz=dreal(z)

2673

2735

rz=dreal(z)

iz=dimag(z)

2674

2736

iz=dimag(z)

if ( iz == 0._8 ) then

2675

2737

if ( iz == 0._8 ) then

! z is real

2676

2738

! z is real

fvn_z_asin=fvn_z_asin_real(rz)

2677

2739

fvn_z_asin=fvn_z_asin_real(rz)

return

2678

2740

return

end if

2679

2741

end if

2680

2742

x=dabs(rz)

2681

2743

x=dabs(rz)

y=dabs(iz)

2682

2744

y=dabs(iz)

r=fvn_d_hypot(x+1.,y)

2683

2745

r=fvn_d_hypot(x+1.,y)

s=fvn_d_hypot(x-1.,y)

2684

2746

s=fvn_d_hypot(x-1.,y)

a=0.5*(r + s)

2685

2747

a=0.5*(r + s)

b=x/a

2686

2748

b=x/a

y2=y*y

2687

2749

y2=y*y

2688

2750

if (b <= b_crossover) then

2689

2751

if (b <= b_crossover) then

r_res=dasin(b)

2690

2752

r_res=dasin(b)

else

2691

2753

else

if (x <= 1.) then

2692

2754

if (x <= 1.) then

d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))

2693

2755

d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))

r_res=datan(x/dsqrt(d))

2694

2756

r_res=datan(x/dsqrt(d))

else

2695

2757

else

apx=a+x

2696

2758

apx=a+x

d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))

2697

2759

d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))

r_res=datan(x/(y*dsqrt(d)));

2698

2760

r_res=datan(x/(y*dsqrt(d)));

end if

2699

2761

end if

2700

2762

end if

2701

2763

if (a <= a_crossover) then

2702

2764

if (a <= a_crossover) then

if (x < 1.) then

2703

2765

if (x < 1.) then

am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))

2704

2766

am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))

else

2705

2767

else

am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))

2706

2768

am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))

end if

2707

2769

end if

i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));

2708

2770

i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));

else

2709

2771

else

i_res = dlog (a + dsqrt (a*a - 1.));

2710

2772

i_res = dlog (a + dsqrt (a*a - 1.));

end if

2711

2773

end if

r_res=sign(1._8,rz)*r_res

2712

2774

r_res=sign(1._8,rz)*r_res

i_res=sign(1._8,iz)*i_res

2713

2775

i_res=sign(1._8,iz)*i_res

fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res)

2714

2776

fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res)

2715

2777

end function fvn_z_asin

2716

2778

end function fvn_z_asin

2717

2779

function fvn_z_asin_real(r)

2718

2780

function fvn_z_asin_real(r)

! return the double complex arc sinus for a

2719

2781

! return the double complex arc sinus for a

! double precision argument

2720

2782

! double precision argument

implicit none

2721

2783

implicit none

real(kind=8) :: r

2722

2784

real(kind=8) :: r

complex(kind=8) :: fvn_z_asin_real

2723

2785

complex(kind=8) :: fvn_z_asin_real

2724

2786

if (dabs(r)<=1._8) then

2725

2787

if (dabs(r)<=1._8) then

fvn_z_asin_real=dcmplx(dasin(r))

2726

2788

fvn_z_asin_real=dcmplx(dasin(r))

return

2727

2789

return

end if

2728

2790

end if

if (r < 0._8) then

2729

2791

if (r < 0._8) then

fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r))

2730

2792

fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r))

else

2731

2793

else

fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r))

2732

2794

fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r))

end if

2733

2795

end if

end function fvn_z_asin_real

2734

2796

end function fvn_z_asin_real

2735

2797

function fvn_d_acosh(r)

2736

2798

function fvn_d_acosh(r)

! return the arc hyperbolic cosine

2737

2799

! return the arc hyperbolic cosine

implicit none

2738

2800

implicit none

real(kind=8) :: r

2739

2801

real(kind=8) :: r

real(kind=8) :: fvn_d_acosh

2740

2802

real(kind=8) :: fvn_d_acosh

if (r >=1) then

2741

2803

if (r >=1) then

fvn_d_acosh=dlog(r+dsqrt(r*r-1))

2742

2804

fvn_d_acosh=dlog(r+dsqrt(r*r-1))

else

2743

2805

else

!! TODO : Better error handling!!!!!!

2744

2806

!! TODO : Better error handling!!!!!!

stop "Argument to fvn_d_acosh lesser than 1"

2745

2807

stop "Argument to fvn_d_acosh lesser than 1"

end if

2746

2808

end if

end function fvn_d_acosh

2747

2809

end function fvn_d_acosh

2748

2810

function fvn_d_asinh(r)

2749

2811

function fvn_d_asinh(r)

! return the arc hyperbolic sine

2750

2812

! return the arc hyperbolic sine

implicit none

2751

2813

implicit none

real(kind=8) :: r

2752

2814

real(kind=8) :: r

real(kind=8) :: fvn_d_asinh

2753

2815

real(kind=8) :: fvn_d_asinh

fvn_d_asinh=dlog(r+dsqrt(r*r+1))

2754

2816

fvn_d_asinh=dlog(r+dsqrt(r*r+1))

end function fvn_d_asinh

2755

2817

end function fvn_d_asinh

2756

2818

function fvn_d_hypot(a,b)

2757

2819

function fvn_d_hypot(a,b)

implicit none

2758

2820

implicit none

! return the euclidian norm of vector(a,b)

2759

2821

! return the euclidian norm of vector(a,b)

real(kind=8) :: a,b

2760

2822

real(kind=8) :: a,b

real(kind=8) :: fvn_d_hypot

2761

2823

real(kind=8) :: fvn_d_hypot

fvn_d_hypot=dsqrt(a*a+b*b)

2762

2824

fvn_d_hypot=dsqrt(a*a+b*b)

end function

2763

2825

end function

2764

2826

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

2765

2827

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

!

2766

2828

!

! SPARSE RESOLUTION

2767

2829

! SPARSE RESOLUTION

!

2768

2830

!

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

2769

2831

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

!

2770

2832

!

! Sparse resolution is done by interfaçing Tim Davi's UMFPACK

2771

2833

! Sparse resolution is done by interfaçing Tim Davi's UMFPACK

! http://www.cise.ufl.edu/research/sparse/SuiteSparse/

2772

2834

! http://www.cise.ufl.edu/research/sparse/SuiteSparse/

! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig

2773

2835

! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig

!

2774

2836

!

! Solve Ax=B using UMFPACK

2775

2837

! Solve Ax=B using UMFPACK

!

2776

2838

!

! Where A is a sparse matrix given in its triplet form

2777

2839

! Where A is a sparse matrix given in its triplet form

! T -> non zero elements

2778

2840

! T -> non zero elements

! Ti,Tj -> row and column index (1-based) of the given elt

2779

2841

! Ti,Tj -> row and column index (1-based) of the given elt

! n : rank of matrix A

2780

2842

! n : rank of matrix A

! nz : number of non zero elts

2781

2843

! nz : number of non zero elts

!

2782

2844

!

! fvn_*_sparse_solve

2783

2845

! fvn_*_sparse_solve

! * = zl : double complex + integer(8)

2784

2846

! * = zl : double complex + integer(8)

! * = zi : double complex + integer(4)

2785

2847

! * = zi : double complex + integer(4)

!

2786

2848

!

subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status)

2787

2849

subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status)

implicit none

2788

2850

implicit none

integer(8), intent(in) :: n,nz

2789

2851

integer(8), intent(in) :: n,nz

complex(8),dimension(nz),intent(in) :: T

2790

2852

complex(8),dimension(nz),intent(in) :: T

integer(8),dimension(nz),intent(in) :: Ti,Tj

2791

2853

integer(8),dimension(nz),intent(in) :: Ti,Tj

complex(8),dimension(n),intent(in) :: B

2792

2854

complex(8),dimension(n),intent(in) :: B

complex(8),dimension(n),intent(out) :: x

2793

2855

complex(8),dimension(n),intent(out) :: x

integer(8), intent(out) :: status

2794

2856

integer(8), intent(out) :: status

2795

2857

integer(8),dimension(:),allocatable :: wTi,wTj

2796

2858

integer(8),dimension(:),allocatable :: wTi,wTj

real(8),dimension(:),allocatable :: Tx,Tz

2797

2859

real(8),dimension(:),allocatable :: Tx,Tz

GITLAB

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