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

! 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 ;-) Nevertheless, you

22

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

! may give credits to quadpack authors.

23

! may give credits to quadpack authors.

!

24

!

! Version 1.1

25

! svn version

26

! June 2007 : added some complex trigonometric functions

!

26

27

!

! TO DO LIST :

27

28

! TO DO LIST :

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

28

29

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

! eigenvalues are given with no particular order.

29

30

! eigenvalues are given with no particular order.

! + Generic interface for fvn_x_name family -> fvn_name

30

31

! + Generic interface for fvn_x_name family -> fvn_name

! + Make some parameters optional, status for example

31

32

! + Make some parameters optional, status for example

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

32

33

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

! + unify quadpack routines

33

34

! + unify quadpack routines

! + ...

34

35

! + ...

!

35

36

!

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

36

37

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

37

38

implicit none

38

39

implicit none

40

! We define pi and i for the module

41

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

42

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

43

44

! All quadpack routines are private to the module

39

45

! All quadpack routines are private to the module

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

40

46

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

41

47

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

42

48

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

43

49

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

44

50

45

51

contains

46

52

contains

47

53

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

48

54

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

!

49

55

!

! Matrix inversion subroutines

50

56

! Matrix inversion subroutines

!

51

57

!

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

52

58

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

subroutine fvn_s_matinv(d,a,inva,status)

53

59

subroutine fvn_s_matinv(d,a,inva,status)

!

54

60

!

! Matrix inversion of a real matrix using BLAS and LAPACK

55

61

! Matrix inversion of a real matrix using BLAS and LAPACK

!

56

62

!

! d (in) : matrix rank

57

63

! d (in) : matrix rank

! a (in) : input matrix

58

64

! a (in) : input matrix

! inva (out) : inversed matrix

59

65

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

60

66

! status (ou) : =0 if something failed

!

61

67

!

integer, intent(in) :: d

62

68

integer, intent(in) :: d

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

63

69

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

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

64

70

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

integer, intent(out) :: status

65

71

integer, intent(out) :: status

66

72

integer, allocatable :: ipiv(:)

67

73

integer, allocatable :: ipiv(:)

real, allocatable :: work(:)

68

74

real, allocatable :: work(:)

real twork(1)

69

75

real twork(1)

integer :: info

70

76

integer :: info

integer :: lwork

71

77

integer :: lwork

72

78

status=1

73

79

status=1

74

80

allocate(ipiv(d))

75

81

allocate(ipiv(d))

! copy a into inva using BLAS

76

82

! copy a into inva using BLAS

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

77

83

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

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

78

84

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

! LU factorization using LAPACK

79

85

! LU factorization using LAPACK

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

80

86

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

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

81

87

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

if (info /= 0) then

82

88

if (info /= 0) then

status=0

83

89

status=0

deallocate(ipiv)

84

90

deallocate(ipiv)

return

85

91

return

end if

86

92

end if

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

87

93

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

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

88

94

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

lwork=int(twork(1))

89

95

lwork=int(twork(1))

allocate(work(lwork))

90

96

allocate(work(lwork))

! Matrix inversion using LAPACK

91

97

! Matrix inversion using LAPACK

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

92

98

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

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

93

99

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

if (info /= 0) then

94

100

if (info /= 0) then

status=0

95

101

status=0

end if

96

102

end if

deallocate(work)

97

103

deallocate(work)

deallocate(ipiv)

98

104

deallocate(ipiv)

end subroutine

99

105

end subroutine

100

106

subroutine fvn_d_matinv(d,a,inva,status)

101

107

subroutine fvn_d_matinv(d,a,inva,status)

!

102

108

!

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

103

109

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

!

104

110

!

! d (in) : matrix rank

105

111

! d (in) : matrix rank

! a (in) : input matrix

106

112

! a (in) : input matrix

! inva (out) : inversed matrix

107

113

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

108

114

! status (ou) : =0 if something failed

!

109

115

!

integer, intent(in) :: d

110

116

integer, intent(in) :: d

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

111

117

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

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

112

118

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

integer, intent(out) :: status

113

119

integer, intent(out) :: status

114

120

integer, allocatable :: ipiv(:)

115

121

integer, allocatable :: ipiv(:)

double precision, allocatable :: work(:)

116

122

double precision, allocatable :: work(:)

double precision :: twork(1)

117

123

double precision :: twork(1)

integer :: info

118

124

integer :: info

integer :: lwork

119

125

integer :: lwork

120

126

status=1

121

127

status=1

122

128

allocate(ipiv(d))

123

129

allocate(ipiv(d))

! copy a into inva using BLAS

124

130

! copy a into inva using BLAS

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

125

131

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

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

126

132

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

! LU factorization using LAPACK

127

133

! LU factorization using LAPACK

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

128

134

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

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

129

135

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

if (info /= 0) then

130

136

if (info /= 0) then

status=0

131

137

status=0

deallocate(ipiv)

132

138

deallocate(ipiv)

return

133

139

return

end if

134

140

end if

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

135

141

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

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

136

142

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

lwork=int(twork(1))

137

143

lwork=int(twork(1))

allocate(work(lwork))

138

144

allocate(work(lwork))

! Matrix inversion using LAPACK

139

145

! Matrix inversion using LAPACK

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

140

146

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

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

141

147

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

if (info /= 0) then

142

148

if (info /= 0) then

status=0

143

149

status=0

end if

144

150

end if

deallocate(work)

145

151

deallocate(work)

deallocate(ipiv)

146

152

deallocate(ipiv)

end subroutine

147

153

end subroutine

148

154

subroutine fvn_c_matinv(d,a,inva,status)

149

155

subroutine fvn_c_matinv(d,a,inva,status)

!

150

156

!

! Matrix inversion of a complex matrix using BLAS and LAPACK

151

157

! Matrix inversion of a complex matrix using BLAS and LAPACK

!

152

158

!

! d (in) : matrix rank

153

159

! d (in) : matrix rank

! a (in) : input matrix

154

160

! a (in) : input matrix

! inva (out) : inversed matrix

155

161

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

156

162

! status (ou) : =0 if something failed

!

157

163

!

integer, intent(in) :: d

158

164

integer, intent(in) :: d

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

159

165

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

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

160

166

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

integer, intent(out) :: status

161

167

integer, intent(out) :: status

162

168

integer, allocatable :: ipiv(:)

163

169

integer, allocatable :: ipiv(:)

complex, allocatable :: work(:)

164

170

complex, allocatable :: work(:)

complex :: twork(1)

165

171

complex :: twork(1)

integer :: info

166

172

integer :: info

integer :: lwork

167

173

integer :: lwork

168

174

status=1

169

175

status=1

170

176

allocate(ipiv(d))

171

177

allocate(ipiv(d))

! copy a into inva using BLAS

172

178

! copy a into inva using BLAS

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

173

179

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

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

174

180

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

175

181

! LU factorization using LAPACK

176

182

! LU factorization using LAPACK

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

177

183

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

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

178

184

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

if (info /= 0) then

179

185

if (info /= 0) then

status=0

180

186

status=0

deallocate(ipiv)

181

187

deallocate(ipiv)

return

182

188

return

end if

183

189

end if

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

184

190

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

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

185

191

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

lwork=int(twork(1))

186

192

lwork=int(twork(1))

allocate(work(lwork))

187

193

allocate(work(lwork))

! Matrix inversion using LAPACK

188

194

! Matrix inversion using LAPACK

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

189

195

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

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

190

196

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

if (info /= 0) then

191

197

if (info /= 0) then

status=0

192

198

status=0

end if

193

199

end if

deallocate(work)

194

200

deallocate(work)

deallocate(ipiv)

195

201

deallocate(ipiv)

end subroutine

196

202

end subroutine

197

203

subroutine fvn_z_matinv(d,a,inva,status)

198

204

subroutine fvn_z_matinv(d,a,inva,status)

!

199

205

!

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

200

206

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

!

201

207

!

! d (in) : matrix rank

202

208

! d (in) : matrix rank

! a (in) : input matrix

203

209

! a (in) : input matrix

! inva (out) : inversed matrix

204

210

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

205

211

! status (ou) : =0 if something failed

!

206

212

!

integer, intent(in) :: d

207

213

integer, intent(in) :: d

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

208

214

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

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

209

215

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

integer, intent(out) :: status

210

216

integer, intent(out) :: status

211

217

integer, allocatable :: ipiv(:)

212

218

integer, allocatable :: ipiv(:)

double complex, allocatable :: work(:)

213

219

double complex, allocatable :: work(:)

double complex :: twork(1)

214

220

double complex :: twork(1)

integer :: info

215

221

integer :: info

integer :: lwork

216

222

integer :: lwork

217

223

status=1

218

224

status=1

219

225

allocate(ipiv(d))

220

226

allocate(ipiv(d))

! copy a into inva using BLAS

221

227

! copy a into inva using BLAS

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

222

228

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

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

223

229

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

224

230

! LU factorization using LAPACK

225

231

! LU factorization using LAPACK

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

226

232

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

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

227

233

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

if (info /= 0) then

228

234

if (info /= 0) then

status=0

229

235

status=0

deallocate(ipiv)

230

236

deallocate(ipiv)

return

231

237

return

end if

232

238

end if

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

233

239

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

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

234

240

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

lwork=int(twork(1))

235

241

lwork=int(twork(1))

allocate(work(lwork))

236

242

allocate(work(lwork))

! Matrix inversion using LAPACK

237

243

! Matrix inversion using LAPACK

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

238

244

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

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

239

245

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

if (info /= 0) then

240

246

if (info /= 0) then

status=0

241

247

status=0

end if

242

248

end if

deallocate(work)

243

249

deallocate(work)

deallocate(ipiv)

244

250

deallocate(ipiv)

end subroutine

245

251

end subroutine

246

252

247

253

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

248

254

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

!

249

255

!

! Determinants

250

256

! Determinants

!

251

257

!

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

252

258

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

function fvn_s_det(d,a,status)

253

259

function fvn_s_det(d,a,status)

!

254

260

!

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

255

261

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

!

256

262

!

! d (in) : matrix rank

257

263

! d (in) : matrix rank

! a (in) : The Matrix

258

264

! a (in) : The Matrix

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

259

265

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

!

260

266

!

integer, intent(in) :: d

261

267

integer, intent(in) :: d

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

262

268

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

integer, intent(out) :: status

263

269

integer, intent(out) :: status

real :: fvn_s_det

264

270

real :: fvn_s_det

265

271

real, allocatable :: wc_a(:,:)

266

272

real, allocatable :: wc_a(:,:)

integer, allocatable :: ipiv(:)

267

273

integer, allocatable :: ipiv(:)

integer :: info,i

268

274

integer :: info,i

269

275

status=1

270

276

status=1

allocate(wc_a(d,d))

271

277

allocate(wc_a(d,d))

allocate(ipiv(d))

272

278

allocate(ipiv(d))

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

273

279

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

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

274

280

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

if (info/= 0) then

275

281

if (info/= 0) then

status=0

276

282

status=0

fvn_s_det=0.e0

277

283

fvn_s_det=0.e0

deallocate(ipiv)

278

284

deallocate(ipiv)

deallocate(wc_a)

279

285

deallocate(wc_a)

return

280

286

return

end if

281

287

end if

fvn_s_det=1.e0

282

288

fvn_s_det=1.e0

do i=1,d

283

289

do i=1,d

if (ipiv(i)==i) then

284

290

if (ipiv(i)==i) then

fvn_s_det=fvn_s_det*wc_a(i,i)

285

291

fvn_s_det=fvn_s_det*wc_a(i,i)

else

286

292

else

fvn_s_det=-fvn_s_det*wc_a(i,i)

287

293

fvn_s_det=-fvn_s_det*wc_a(i,i)

end if

288

294

end if

end do

289

295

end do

deallocate(ipiv)

290

296

deallocate(ipiv)

deallocate(wc_a)

291

297

deallocate(wc_a)

292

298

end function

293

299

end function

294

300

function fvn_d_det(d,a,status)

295

301

function fvn_d_det(d,a,status)

!

296

302

!

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

297

303

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

!

298

304

!

! d (in) : matrix rank

299

305

! d (in) : matrix rank

! a (in) : The Matrix

300

306

! a (in) : The Matrix

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

301

307

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

!

302

308

!

integer, intent(in) :: d

303

309

integer, intent(in) :: d

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

304

310

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

integer, intent(out) :: status

305

311

integer, intent(out) :: status

double precision :: fvn_d_det

306

312

double precision :: fvn_d_det

307

313

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

308

314

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

integer, allocatable :: ipiv(:)

309

315

integer, allocatable :: ipiv(:)

integer :: info,i

310

316

integer :: info,i

311

317

status=1

312

318

status=1

allocate(wc_a(d,d))

313

319

allocate(wc_a(d,d))

allocate(ipiv(d))

314

320

allocate(ipiv(d))

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

315

321

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

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

316

322

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

if (info/= 0) then

317

323

if (info/= 0) then

status=0

318

324

status=0

fvn_d_det=0.d0

319

325

fvn_d_det=0.d0

deallocate(ipiv)

320

326

deallocate(ipiv)

deallocate(wc_a)

321

327

deallocate(wc_a)

return

322

328

return

end if

323

329

end if

fvn_d_det=1.d0

324

330

fvn_d_det=1.d0

do i=1,d

325

331

do i=1,d

if (ipiv(i)==i) then

326

332

if (ipiv(i)==i) then

fvn_d_det=fvn_d_det*wc_a(i,i)

327

333

fvn_d_det=fvn_d_det*wc_a(i,i)

else

328

334

else

fvn_d_det=-fvn_d_det*wc_a(i,i)

329

335

fvn_d_det=-fvn_d_det*wc_a(i,i)

end if

330

336

end if

end do

331

337

end do

deallocate(ipiv)

332

338

deallocate(ipiv)

deallocate(wc_a)

333

339

deallocate(wc_a)

334

340

end function

335

341

end function

336

342

function fvn_c_det(d,a,status) !

337

343

function fvn_c_det(d,a,status) !

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

338

344

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

!

339

345

!

! d (in) : matrix rank

340

346

! d (in) : matrix rank

! a (in) : The Matrix

341

347

! a (in) : The Matrix

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

342

348

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

!

343

349

!

integer, intent(in) :: d

344

350

integer, intent(in) :: d

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

345

351

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

integer, intent(out) :: status

346

352

integer, intent(out) :: status

complex :: fvn_c_det

347

353

complex :: fvn_c_det

348

354

complex, allocatable :: wc_a(:,:)

349

355

complex, allocatable :: wc_a(:,:)

integer, allocatable :: ipiv(:)

350

356

integer, allocatable :: ipiv(:)

integer :: info,i

351

357

integer :: info,i

352

358

status=1

353

359

status=1

allocate(wc_a(d,d))

354

360

allocate(wc_a(d,d))

allocate(ipiv(d))

355

361

allocate(ipiv(d))

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

356

362

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

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

357

363

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

if (info/= 0) then

358

364

if (info/= 0) then

status=0

359

365

status=0

fvn_c_det=(0.e0,0.e0)

360

366

fvn_c_det=(0.e0,0.e0)

deallocate(ipiv)

361

367

deallocate(ipiv)

deallocate(wc_a)

362

368

deallocate(wc_a)

return

363

369

return

end if

364

370

end if

fvn_c_det=(1.e0,0.e0)

365

371

fvn_c_det=(1.e0,0.e0)

do i=1,d

366

372

do i=1,d

if (ipiv(i)==i) then

367

373

if (ipiv(i)==i) then

fvn_c_det=fvn_c_det*wc_a(i,i)

368

374

fvn_c_det=fvn_c_det*wc_a(i,i)

else

369

375

else

fvn_c_det=-fvn_c_det*wc_a(i,i)

370

376

fvn_c_det=-fvn_c_det*wc_a(i,i)

end if

371

377

end if

end do

372

378

end do

deallocate(ipiv)

373

379

deallocate(ipiv)

deallocate(wc_a)

374

380

deallocate(wc_a)

375

381

end function

376

382

end function

377

383

function fvn_z_det(d,a,status)

378

384

function fvn_z_det(d,a,status)

!

379

385

!

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

380

386

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

!

381

387

!

! d (in) : matrix rank

382

388

! d (in) : matrix rank

! a (in) : The Matrix

383

389

! a (in) : The Matrix

! det (out) : determinant

384

390

! det (out) : determinant

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

385

391

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

!

386

392

!

integer, intent(in) :: d

387

393

integer, intent(in) :: d

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

388

394

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

integer, intent(out) :: status

389

395

integer, intent(out) :: status

double complex :: fvn_z_det

390

396

double complex :: fvn_z_det

391

397

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

392

398

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

integer, allocatable :: ipiv(:)

393

399

integer, allocatable :: ipiv(:)

integer :: info,i

394

400

integer :: info,i

395

401

status=1

396

402

status=1

allocate(wc_a(d,d))

397

403

allocate(wc_a(d,d))

allocate(ipiv(d))

398

404

allocate(ipiv(d))

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

399

405

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

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

400

406

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

if (info/= 0) then

401

407

if (info/= 0) then

status=0

402

408

status=0

fvn_z_det=(0.d0,0.d0)

403

409

fvn_z_det=(0.d0,0.d0)

deallocate(ipiv)

404

410

deallocate(ipiv)

deallocate(wc_a)

405

411

deallocate(wc_a)

return

406

412

return

end if

407

413

end if

fvn_z_det=(1.d0,0.d0)

408

414

fvn_z_det=(1.d0,0.d0)

do i=1,d

409

415

do i=1,d

if (ipiv(i)==i) then

410

416

if (ipiv(i)==i) then

fvn_z_det=fvn_z_det*wc_a(i,i)

411

417

fvn_z_det=fvn_z_det*wc_a(i,i)

else

412

418

else

fvn_z_det=-fvn_z_det*wc_a(i,i)

413

419

fvn_z_det=-fvn_z_det*wc_a(i,i)

end if

414

420

end if

end do

415

421

end do

deallocate(ipiv)

416

422

deallocate(ipiv)

deallocate(wc_a)

417

423

deallocate(wc_a)

418

424

end function

419

425

end function

420

426

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

421

427

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

!

422

428

!

! Condition test

423

429

! Condition test

!

424

430

!

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

425

431

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

! 1-norm

426

432

! 1-norm

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

427

433

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

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

428

434

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

!

429

435

!

subroutine fvn_s_matcon(d,a,rcond,status)

430

436

subroutine fvn_s_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

431

437

! Matrix condition (reciprocal of condition number)

!

432

438

!

! d (in) : matrix rank

433

439

! d (in) : matrix rank

! a (in) : The Matrix

434

440

! a (in) : The Matrix

! rcond (out) : guess what

435

441

! rcond (out) : guess what

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

436

442

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

!

437

443

!

integer, intent(in) :: d

438

444

integer, intent(in) :: d

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

439

445

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

real, intent(out) :: rcond

440

446

real, intent(out) :: rcond

integer, intent(out) :: status

441

447

integer, intent(out) :: status

442

448

real, allocatable :: work(:)

443

449

real, allocatable :: work(:)

integer, allocatable :: iwork(:)

444

450

integer, allocatable :: iwork(:)

real :: anorm

445

451

real :: anorm

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

446

452

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

integer :: info

447

453

integer :: info

integer, allocatable :: ipiv(:)

448

454

integer, allocatable :: ipiv(:)

449

455

real, external :: slange

450

456

real, external :: slange

451

457

452

458

status=1

453

459

status=1

454

460

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

455

461

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

456

462

allocate(wc_a(d,d))

457

463

allocate(wc_a(d,d))

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

458

464

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

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

459

465

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

460

466

allocate(ipiv(d))

461

467

allocate(ipiv(d))

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

462

468

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

if (info /= 0) then

463

469

if (info /= 0) then

status=0

464

470

status=0

deallocate(ipiv)

465

471

deallocate(ipiv)

deallocate(wc_a)

466

472

deallocate(wc_a)

return

467

473

return

end if

468

474

end if

allocate(work(4*d))

469

475

allocate(work(4*d))

allocate(iwork(d))

470

476

allocate(iwork(d))

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

471

477

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

if (info /= 0) then

472

478

if (info /= 0) then

status=0

473

479

status=0

end if

474

480

end if

deallocate(iwork)

475

481

deallocate(iwork)

deallocate(work)

476

482

deallocate(work)

deallocate(ipiv)

477

483

deallocate(ipiv)

deallocate(wc_a)

478

484

deallocate(wc_a)

479

485

end subroutine

480

486

end subroutine

481

487

subroutine fvn_d_matcon(d,a,rcond,status)

482

488

subroutine fvn_d_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

483

489

! Matrix condition (reciprocal of condition number)

!

484

490

!

! d (in) : matrix rank

485

491

! d (in) : matrix rank

! a (in) : The Matrix

486

492

! a (in) : The Matrix

! rcond (out) : guess what

487

493

! rcond (out) : guess what

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

488

494

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

!

489

495

!

integer, intent(in) :: d

490

496

integer, intent(in) :: d

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

491

497

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

double precision, intent(out) :: rcond

492

498

double precision, intent(out) :: rcond

integer, intent(out) :: status

493

499

integer, intent(out) :: status

494

500

double precision, allocatable :: work(:)

495

501

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

496

502

integer, allocatable :: iwork(:)

double precision :: anorm

497

503

double precision :: anorm

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

498

504

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

integer :: info

499

505

integer :: info

integer, allocatable :: ipiv(:)

500

506

integer, allocatable :: ipiv(:)

501

507

double precision, external :: dlange

502

508

double precision, external :: dlange

503

509

504

510

status=1

505

511

status=1

506

512

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

507

513

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

508

514

allocate(wc_a(d,d))

509

515

allocate(wc_a(d,d))

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

510

516

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

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

511

517

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

512

518

allocate(ipiv(d))

513

519

allocate(ipiv(d))

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

514

520

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

if (info /= 0) then

515

521

if (info /= 0) then

status=0

516

522

status=0

deallocate(ipiv)

517

523

deallocate(ipiv)

deallocate(wc_a)

518

524

deallocate(wc_a)

return

519

525

return

end if

520

526

end if

521

527

allocate(work(4*d))

522

528

allocate(work(4*d))

allocate(iwork(d))

523

529

allocate(iwork(d))

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

524

530

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

if (info /= 0) then

525

531

if (info /= 0) then

status=0

526

532

status=0

end if

527

533

end if

deallocate(iwork)

528

534

deallocate(iwork)

deallocate(work)

529

535

deallocate(work)

deallocate(ipiv)

530

536

deallocate(ipiv)

deallocate(wc_a)

531

537

deallocate(wc_a)

532

538

end subroutine

533

539

end subroutine

534

540

subroutine fvn_c_matcon(d,a,rcond,status)

535

541

subroutine fvn_c_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

536

542

! Matrix condition (reciprocal of condition number)

!

537

543

!

! d (in) : matrix rank

538

544

! d (in) : matrix rank

! a (in) : The Matrix

539

545

! a (in) : The Matrix

! rcond (out) : guess what

540

546

! rcond (out) : guess what

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

541

547

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

!

542

548

!

integer, intent(in) :: d

543

549

integer, intent(in) :: d

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

544

550

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

real, intent(out) :: rcond

545

551

real, intent(out) :: rcond

integer, intent(out) :: status

546

552

integer, intent(out) :: status

547

553

real, allocatable :: rwork(:)

548

554

real, allocatable :: rwork(:)

complex, allocatable :: work(:)

549

555

complex, allocatable :: work(:)

integer, allocatable :: iwork(:)

550

556

integer, allocatable :: iwork(:)

real :: anorm

551

557

real :: anorm

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

552

558

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

integer :: info

553

559

integer :: info

integer, allocatable :: ipiv(:)

554

560

integer, allocatable :: ipiv(:)

555

561

real, external :: clange

556

562

real, external :: clange

557

563

558

564

status=1

559

565

status=1

560

566

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

561

567

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

562

568

allocate(wc_a(d,d))

563

569

allocate(wc_a(d,d))

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

564

570

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

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

565

571

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

566

572

allocate(ipiv(d))

567

573

allocate(ipiv(d))

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

568

574

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

if (info /= 0) then

569

575

if (info /= 0) then

status=0

570

576

status=0

deallocate(ipiv)

571

577

deallocate(ipiv)

deallocate(wc_a)

572

578

deallocate(wc_a)

return

573

579

return

end if

574

580

end if

allocate(work(2*d))

575

581

allocate(work(2*d))

allocate(rwork(2*d))

576

582

allocate(rwork(2*d))

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

577

583

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

if (info /= 0) then

578

584

if (info /= 0) then

status=0

579

585

status=0

end if

580

586

end if

deallocate(rwork)

581

587

deallocate(rwork)

deallocate(work)

582

588

deallocate(work)

deallocate(ipiv)

583

589

deallocate(ipiv)

deallocate(wc_a)

584

590

deallocate(wc_a)

end subroutine

585

591

end subroutine

586

592

subroutine fvn_z_matcon(d,a,rcond,status)

587

593

subroutine fvn_z_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

588

594

! Matrix condition (reciprocal of condition number)

!

589

595

!

! d (in) : matrix rank

590

596

! d (in) : matrix rank

! a (in) : The Matrix

591

597

! a (in) : The Matrix

! rcond (out) : guess what

592

598

! rcond (out) : guess what

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

593

599

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

!

594

600

!

integer, intent(in) :: d

595

601

integer, intent(in) :: d

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

596

602

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

double precision, intent(out) :: rcond

597

603

double precision, intent(out) :: rcond

integer, intent(out) :: status

598

604

integer, intent(out) :: status

599

605

double complex, allocatable :: work(:)

600

606

double complex, allocatable :: work(:)

double precision, allocatable :: rwork(:)

601

607

double precision, allocatable :: rwork(:)

double precision :: anorm

602

608

double precision :: anorm

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

603

609

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

integer :: info

604

610

integer :: info

integer, allocatable :: ipiv(:)

605

611

integer, allocatable :: ipiv(:)

606

612

double precision, external :: zlange

607

613

double precision, external :: zlange

608

614

609

615

status=1

610

616

status=1

611

617

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

612

618

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

613

619

allocate(wc_a(d,d))

614

620

allocate(wc_a(d,d))

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

615

621

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

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

616

622

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

617

623

allocate(ipiv(d))

618

624

allocate(ipiv(d))

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

619

625

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

if (info /= 0) then

620

626

if (info /= 0) then

status=0

621

627

status=0

deallocate(ipiv)

622

628

deallocate(ipiv)

deallocate(wc_a)

623

629

deallocate(wc_a)

return

624

630

return

end if

625

631

end if

626

632

allocate(work(2*d))

627

633

allocate(work(2*d))

allocate(rwork(2*d))

628

634

allocate(rwork(2*d))

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

629

635

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

if (info /= 0) then

630

636

if (info /= 0) then

status=0

631

637

status=0

end if

632

638

end if

deallocate(rwork)

633

639

deallocate(rwork)

deallocate(work)

634

640

deallocate(work)

deallocate(ipiv)

635

641

deallocate(ipiv)

deallocate(wc_a)

636

642

deallocate(wc_a)

end subroutine

637

643

end subroutine

638

644

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

639

645

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

!

640

646

!

! Valeurs propres/ Vecteurs propre

641

647

! Valeurs propres/ Vecteurs propre

!

642

648

!

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

643

649

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

644

650

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

645

651

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

!

646

652

!

! integer d (in) : matrice rank

647

653

! integer d (in) : matrice rank

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

648

654

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

! complex evala(d) (out) : eigenvalues

649

655

! complex evala(d) (out) : eigenvalues

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

650

656

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

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

651

657

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

!

652

658

!

! interfacing Lapack routine SGEEV

653

659

! interfacing Lapack routine SGEEV

654

660

integer, intent(in) :: d

655

661

integer, intent(in) :: d

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

656

662

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

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

657

663

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

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

658

664

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

integer, intent(out) :: status

659

665

integer, intent(out) :: status

660

666

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

661

667

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

integer :: info

662

668

integer :: info

integer :: lwork

663

669

integer :: lwork

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

664

670

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

real :: vl ! unused but necessary for the call

665

671

real :: vl ! unused but necessary for the call

real, allocatable :: vr(:,:)

666

672

real, allocatable :: vr(:,:)

real, allocatable :: work(:)

667

673

real, allocatable :: work(:)

real :: twork(1)

668

674

real :: twork(1)

integer i

669

675

integer i

integer j

670

676

integer j

671

677

! making a working copy of a

672

678

! making a working copy of a

allocate(wc_a(d,d))

673

679

allocate(wc_a(d,d))

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

674

680

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

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

675

681

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

676

682

allocate(wr(d))

677

683

allocate(wr(d))

allocate(wi(d))

678

684

allocate(wi(d))

allocate(vr(d,d))

679

685

allocate(vr(d,d))

! query optimal work size

680

686

! query optimal work size

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

681

687

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

lwork=int(twork(1))

682

688

lwork=int(twork(1))

allocate(work(lwork))

683

689

allocate(work(lwork))

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

684

690

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

685

691

if (info /= 0) then

686

692

if (info /= 0) then

status=0

687

693

status=0

deallocate(work)

688

694

deallocate(work)

deallocate(vr)

689

695

deallocate(vr)

deallocate(wi)

690

696

deallocate(wi)

deallocate(wr)

691

697

deallocate(wr)

deallocate(wc_a)

692

698

deallocate(wc_a)

return

693

699

return

end if

694

700

end if

695

701

! now fill in the results

696

702

! now fill in the results

i=1

697

703

i=1

do while(i<=d)

698

704

do while(i<=d)

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

699

705

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

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

700

706

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

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

701

707

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

else ! eigenvalue is complex

702

708

else ! eigenvalue is complex

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

703

709

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

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

704

710

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

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

705

711

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

i=i+1

706

712

i=i+1

end if

707

713

end if

i=i+1

708

714

i=i+1

enddo

709

715

enddo

deallocate(work)

710

716

deallocate(work)

deallocate(vr)

711

717

deallocate(vr)

deallocate(wi)

712

718

deallocate(wi)

deallocate(wr)

713

719

deallocate(wr)

deallocate(wc_a)

714

720

deallocate(wc_a)

715

721

end subroutine

716

722

end subroutine

717

723

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

718

724

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

!

719

725

!

! integer d (in) : matrice rank

720

726

! integer d (in) : matrice rank

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

721

727

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

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

722

728

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

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

723

729

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

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

724

730

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

!

725

731

!

! interfacing Lapack routine DGEEV

726

732

! interfacing Lapack routine DGEEV

integer, intent(in) :: d

727

733

integer, intent(in) :: d

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

728

734

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

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

729

735

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

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

730

736

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

integer, intent(out) :: status

731

737

integer, intent(out) :: status

732

738

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

733

739

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

integer :: info

734

740

integer :: info

integer :: lwork

735

741

integer :: lwork

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

736

742

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

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

737

743

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

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

738

744

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

double precision, allocatable :: work(:)

739

745

double precision, allocatable :: work(:)

double precision :: twork(1)

740

746

double precision :: twork(1)

integer i

741

747

integer i

integer j

742

748

integer j

743

749

! making a working copy of a

744

750

! making a working copy of a

allocate(wc_a(d,d))

745

751

allocate(wc_a(d,d))

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

746

752

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

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

747

753

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

748

754

allocate(wr(d))

749

755

allocate(wr(d))

allocate(wi(d))

750

756

allocate(wi(d))

allocate(vr(d,d))

751

757

allocate(vr(d,d))

! query optimal work size

752

758

! query optimal work size

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

753

759

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

lwork=int(twork(1))

754

760

lwork=int(twork(1))

allocate(work(lwork))

755

761

allocate(work(lwork))

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

756

762

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

757

763

if (info /= 0) then

758

764

if (info /= 0) then

status=0

759

765

status=0

deallocate(work)

760

766

deallocate(work)

deallocate(vr)

761

767

deallocate(vr)

deallocate(wi)

762

768

deallocate(wi)

deallocate(wr)

763

769

deallocate(wr)

deallocate(wc_a)

764

770

deallocate(wc_a)

return

765

771

return

end if

766

772

end if

767

773

! now fill in the results

768

774

! now fill in the results

i=1

769

775

i=1

do while(i<=d)

770

776

do while(i<=d)

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

771

777

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

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

772

778

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

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

773

779

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

else ! eigenvalue is complex

774

780

else ! eigenvalue is complex

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

775

781

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

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

776

782

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

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

777

783

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

i=i+1

778

784

i=i+1

end if

779

785

end if

i=i+1

780

786

i=i+1

enddo

781

787

enddo

782

788

deallocate(work)

783

789

deallocate(work)

deallocate(vr)

784

790

deallocate(vr)

deallocate(wi)

785

791

deallocate(wi)

deallocate(wr)

786

792

deallocate(wr)

deallocate(wc_a)

787

793

deallocate(wc_a)

788

794

end subroutine

789

795

end subroutine

790

796

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

791

797

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

!

792

798

!

! integer d (in) : matrice rank

793

799

! integer d (in) : matrice rank

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

794

800

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

! complex evala(d) (out) : eigenvalues

795

801

! complex evala(d) (out) : eigenvalues

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

796

802

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

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

797

803

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

!

798

804

!

! interfacing Lapack routine CGEEV

799

805

! interfacing Lapack routine CGEEV

800

806

integer, intent(in) :: d

801

807

integer, intent(in) :: d

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

802

808

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

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

803

809

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

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

804

810

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

integer, intent(out) :: status

805

811

integer, intent(out) :: status

806

812

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

807

813

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

integer :: info

808

814

integer :: info

integer :: lwork

809

815

integer :: lwork

complex, allocatable :: work(:)

810

816

complex, allocatable :: work(:)

complex :: twork(1)

811

817

complex :: twork(1)

real, allocatable :: rwork(:)

812

818

real, allocatable :: rwork(:)

complex :: vl ! unused but necessary for the call

813

819

complex :: vl ! unused but necessary for the call

814

820

status=1

815

821

status=1

816

822

! making a working copy of a

817

823

! making a working copy of a

allocate(wc_a(d,d))

818

824

allocate(wc_a(d,d))

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

819

825

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

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

820

826

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

821

827

822

828

! query optimal work size

823

829

! query optimal work size

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

824

830

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

lwork=int(twork(1))

825

831

lwork=int(twork(1))

allocate(work(lwork))

826

832

allocate(work(lwork))

allocate(rwork(2*d))

827

833

allocate(rwork(2*d))

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

828

834

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

829

835

if (info /= 0) then

830

836

if (info /= 0) then

status=0

831

837

status=0

end if

832

838

end if

deallocate(rwork)

833

839

deallocate(rwork)

deallocate(work)

834

840

deallocate(work)

deallocate(wc_a)

835

841

deallocate(wc_a)

836

842

end subroutine

837

843

end subroutine

838

844

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

839

845

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

!

840

846

!

! integer d (in) : matrice rank

841

847

! integer d (in) : matrice rank

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

842

848

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

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

843

849

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

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

844

850

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

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

845

851

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

!

846

852

!

! interfacing Lapack routine ZGEEV

847

853

! interfacing Lapack routine ZGEEV

848

854

integer, intent(in) :: d

849

855

integer, intent(in) :: d

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

850

856

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

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

851

857

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

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

852

858

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

integer, intent(out) :: status

853

859

integer, intent(out) :: status

854

860

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

855

861

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

integer :: info

856

862

integer :: info

integer :: lwork

857

863

integer :: lwork

double complex, allocatable :: work(:)

858

864

double complex, allocatable :: work(:)

double complex :: twork(1)

859

865

double complex :: twork(1)

double precision, allocatable :: rwork(:)

860

866

double precision, allocatable :: rwork(:)

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

861

867

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

862

868

status=1

863

869

status=1

864

870

! making a working copy of a

865

871

! making a working copy of a

allocate(wc_a(d,d))

866

872

allocate(wc_a(d,d))

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

867

873

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

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

868

874

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

869

875

! query optimal work size

870

876

! query optimal work size

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

871

877

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

lwork=int(twork(1))

872

878

lwork=int(twork(1))

allocate(work(lwork))

873

879

allocate(work(lwork))

allocate(rwork(2*d))

874

880

allocate(rwork(2*d))

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

875

881

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

876

882

if (info /= 0) then

877

883

if (info /= 0) then

status=0

878

884

status=0

end if

879

885

end if

deallocate(rwork)

880

886

deallocate(rwork)

deallocate(work)

881

887

deallocate(work)

deallocate(wc_a)

882

888

deallocate(wc_a)

883

889

end subroutine

884

890

end subroutine

885

891

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

886

892

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

!

887

893

!

! Akima spline interpolation and spline evaluation

888

894

! Akima spline interpolation and spline evaluation

!

889

895

!

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

890

896

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

891

897

! Single precision

892

898

! Single precision

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

893

899

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

implicit none

894

900

implicit none

integer, intent(in) :: n

895

901

integer, intent(in) :: n

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

896

902

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

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

897

903

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

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

898

904

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

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

899

905

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

900

906

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

901

907

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

real :: wi_1,wi

902

908

real :: wi_1,wi

integer :: i

903

909

integer :: i

real :: dx,a,b

904

910

real :: dx,a,b

905

911

! br is just a copy of x

906

912

! br is just a copy of x

br(:)=x(:)

907

913

br(:)=x(:)

908

914

allocate(var(n))

909

915

allocate(var(n))

allocate(z(n))

910

916

allocate(z(n))

! evaluate the variations

911

917

! evaluate the variations

do i=1, n-1

912

918

do i=1, n-1

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

913

919

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

end do

914

920

end do

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

915

921

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

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

916

922

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

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

917

923

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

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

918

924

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

919

925

do i = 1, n

920

926

do i = 1, n

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

921

927

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

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

922

928

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

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

923

929

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

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

924

930

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

else

925

931

else

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

926

932

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

end if

927

933

end if

end do

928

934

end do

929

935

do i=1, n-1

930

936

do i=1, n-1

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

931

937

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

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

932

938

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

933

939

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

co(1,i)=y(i)

934

940

co(1,i)=y(i)

co(2,i)=z(i)

935

941

co(2,i)=z(i)

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

936

942

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

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

937

943

!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

938

944

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 !

939

945

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

940

946

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

941

947

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

end do

942

948

end do

co(1,n)=y(n)

943

949

co(1,n)=y(n)

co(2,n)=z(n)

944

950

co(2,n)=z(n)

co(3,n)=0.e0

945

951

co(3,n)=0.e0

co(4,n)=0.e0

946

952

co(4,n)=0.e0

947

953

deallocate(z)

948

954

deallocate(z)

deallocate(var)

949

955

deallocate(var)

950

956

end subroutine

951

957

end subroutine

952

958

! Double precision

953

959

! Double precision

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

954

960

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

955

961

implicit none

956

962

implicit none

integer, intent(in) :: n

957

963

integer, intent(in) :: n

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

958

964

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

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

959

965

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

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

960

966

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

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

961

967

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

962

968

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

963

969

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

double precision :: wi_1,wi

964

970

double precision :: wi_1,wi

integer :: i

965

971

integer :: i

double precision :: dx,a,b

966

972

double precision :: dx,a,b

967

973

! br is just a copy of x

968

974

! br is just a copy of x

br(:)=x(:)

969

975

br(:)=x(:)

970

976

allocate(var(n))

971

977

allocate(var(n))

allocate(z(n))

972

978

allocate(z(n))

! evaluate the variations

973

979

! evaluate the variations

do i=1, n-1

974

980

do i=1, n-1

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

975

981

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

end do

976

982

end do

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

977

983

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

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

978

984

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

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

979

985

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

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

980

986

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

981

987

do i = 1, n

982

988

do i = 1, n

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

983

989

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

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

984

990

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

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

985

991

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

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

986

992

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

else

987

993

else

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

988

994

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

end if

989

995

end if

end do

990

996

end do

991

997

do i=1, n-1

992

998

do i=1, n-1

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

993

999

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

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

994

1000

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

995

1001

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

co(1,i)=y(i)

996

1002

co(1,i)=y(i)

co(2,i)=z(i)

997

1003

co(2,i)=z(i)

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

998

1004

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

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

999

1005

!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

1000

1006

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 !

1001

1007

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

1002

1008

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

1003

1009

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

end do

1004

1010

end do

co(1,n)=y(n)

1005

1011

co(1,n)=y(n)

co(2,n)=z(n)

1006

1012

co(2,n)=z(n)

co(3,n)=0.d0

1007

1013

co(3,n)=0.d0

co(4,n)=0.d0

1008

1014

co(4,n)=0.d0

1009

1015

deallocate(z)

1010

1016

deallocate(z)

deallocate(var)

1011

1017

deallocate(var)

1012

1018

end subroutine

1013

1019

end subroutine

1014

1020

!

1015

1021

!

! Single precision spline evaluation

1016

1022

! Single precision spline evaluation

!

1017

1023

!

function fvn_s_spline_eval(x,n,br,co)

1018

1024

function fvn_s_spline_eval(x,n,br,co)

implicit none

1019

1025

implicit none

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

1020

1026

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

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

1021

1027

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

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

1022

1028

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

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

1023

1029

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

real :: fvn_s_spline_eval

1024

1030

real :: fvn_s_spline_eval

1025

1031

integer :: i

1026

1032

integer :: i

real :: dx

1027

1033

real :: dx

1028

1034

if (x<=br(1)) then

1029

1035

if (x<=br(1)) then

i=1

1030

1036

i=1

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

1031

1037

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

i=n

1032

1038

i=n

else

1033

1039

else

i=1

1034

1040

i=1

do while(x>=br(i))

1035

1041

do while(x>=br(i))

i=i+1

1036

1042

i=i+1

end do

1037

1043

end do

i=i-1

1038

1044

i=i-1

end if

1039

1045

end if

dx=x-br(i)

1040

1046

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

1041

1047

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

1042

1048

end function

1043

1049

end function

1044

1050

! Double precision spline evaluation

1045

1051

! Double precision spline evaluation

function fvn_d_spline_eval(x,n,br,co)

1046

1052

function fvn_d_spline_eval(x,n,br,co)

implicit none

1047

1053

implicit none

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

1048

1054

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

1049

1055

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

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

1050

1056

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

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

1051

1057

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

double precision :: fvn_d_spline_eval

1052

1058

double precision :: fvn_d_spline_eval

1053

1059

integer :: i

1054

1060

integer :: i

double precision :: dx

1055

1061

double precision :: dx

1056

1062

1057

1063

if (x<=br(1)) then

1058

1064

if (x<=br(1)) then

i=1

1059

1065

i=1

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

1060

1066

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

i=n

1061

1067

i=n

else

1062

1068

else

i=1

1063

1069

i=1

do while(x>=br(i))

1064

1070

do while(x>=br(i))

i=i+1

1065

1071

i=i+1

end do

1066

1072

end do

i=i-1

1067

1073

i=i-1

end if

1068

1074

end if

1069

1075

dx=x-br(i)

1070

1076

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

1071

1077

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

1072

1078

end function

1073

1079

end function

1074

1080

1075

1081

!

1076

1082

!

! Muller

1077

1083

! Muller

!

1078

1084

!

1079

1085

!

1080

1086

!

! William Daniau 2007

1081

1087

! William Daniau 2007

!

1082

1088

!

! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f

1083

1089

! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f

! http://plato.asu.edu/ftp/other_software/muller.f

1084

1090

! http://plato.asu.edu/ftp/other_software/muller.f

!

1085

1091

!

! it can be used as a replacement for imsl routine dzanly with minor changes

1086

1092

! it can be used as a replacement for imsl routine dzanly with minor changes

!

1087

1093

!

!-----------------------------------------------------------------------

1088

1094

!-----------------------------------------------------------------------

!

1089

1095

!

! purpose - zeros of an analytic complex function

1090

1096

! purpose - zeros of an analytic complex function

! using the muller method with deflation

1091

1097

! using the muller method with deflation

!

1092

1098

!

! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,

1093

1099

! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,

! infer,ier)

1094

1100

! infer,ier)

!

1095

1101

!

! arguments f - a complex function subprogram, f(z), written

1096

1102

! arguments f - a complex function subprogram, f(z), written

! by the user specifying the equation whose

1097

1103

! by the user specifying the equation whose

! roots are to be found. f must appear in

1098

1104

! roots are to be found. f must appear in

! an external statement in the calling pro-

1099

1105

! an external statement in the calling pro-

! gram.

1100

1106

! gram.

! eps - 1st stopping criterion. let fp(z)=f(z)/p

1101

1107

! eps - 1st stopping criterion. let fp(z)=f(z)/p

! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))

1102

1108

! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))

! and z(1),...,z(k-1) are previously found

1103

1109

! and z(1),...,z(k-1) are previously found

! roots. if ((cdabs(f(z)).le.eps) .and.

1104

1110

! roots. if ((cdabs(f(z)).le.eps) .and.

! (cdabs(fp(z)).le.eps)), then z is accepted

1105

1111

! (cdabs(fp(z)).le.eps)), then z is accepted

! as a root. (input)

1106

1112

! as a root. (input)

! eps1 - 2nd stopping criterion. a root is accepted

1107

1113

! eps1 - 2nd stopping criterion. a root is accepted

! if two successive approximations to a given

1108

1114

! if two successive approximations to a given

! root agree within eps1. (input)

1109

1115

! root agree within eps1. (input)

! note. if either or both of the stopping

1110

1116

! note. if either or both of the stopping

! criteria are fulfilled, the root is

1111

1117

! criteria are fulfilled, the root is

! accepted.

1112

1118

! accepted.

! kn - the number of known roots which must be stored

1113

1119

! kn - the number of known roots which must be stored

! in x(1),...,x(kn), prior to entry to muller

1114

1120

! in x(1),...,x(kn), prior to entry to muller

! nguess - the number of initial guesses provided. these

1115

1121

! nguess - the number of initial guesses provided. these

! guesses must be stored in x(kn+1),...,

1116

1122

! guesses must be stored in x(kn+1),...,

! x(kn+nguess). nguess must be set equal

1117

1123

! x(kn+nguess). nguess must be set equal

! to zero if no guesses are provided. (input)

1118

1124

! to zero if no guesses are provided. (input)

! n - the number of new roots to be found by

1119

1125

! n - the number of new roots to be found by

! muller (input)

1120

1126

! muller (input)

! x - a complex vector of length kn+n. x(1),...,

1121

1127

! x - a complex vector of length kn+n. x(1),...,

! x(kn) on input must contain any known

1122

1128

! x(kn) on input must contain any known

! roots. x(kn+1),..., x(kn+n) on input may,

1123

1129

! roots. x(kn+1),..., x(kn+n) on input may,

! on user option, contain initial guesses for

1124

1130

! on user option, contain initial guesses for

! the n new roots which are to be computed.

1125

1131

! the n new roots which are to be computed.

! if the user does not provide an initial

1126

1132

! if the user does not provide an initial

! guess, zero is used.

1127

1133

! guess, zero is used.

! on output, x(kn+1),...,x(kn+n) contain the

1128

1134

! on output, x(kn+1),...,x(kn+n) contain the

! approximate roots found by muller.

1129

1135

! approximate roots found by muller.

! itmax - the maximum allowable number of iterations

1130

1136

! itmax - the maximum allowable number of iterations

! per root (input)

1131

1137

! per root (input)

! infer - an integer vector of length kn+n. on

1132

1138

! infer - an integer vector of length kn+n. on

! output infer(j) contains the number of

1133

1139

! output infer(j) contains the number of

! iterations used in finding the j-th root

1134

1140

! iterations used in finding the j-th root

! when convergence was achieved. if

1135

1141

! when convergence was achieved. if

! convergence was not obtained in itmax

1136

1142

! convergence was not obtained in itmax

! iterations, infer(j) will be greater than

1137

1143

! iterations, infer(j) will be greater than

! itmax (output).

1138

1144

! itmax (output).

! ier - error parameter (output)

1139

1145

! ier - error parameter (output)

! warning error

1140

1146

! warning error

! ier = 33 indicates failure to converge with-

1141

1147

! ier = 33 indicates failure to converge with-

! in itmax iterations for at least one of

1142

1148

! in itmax iterations for at least one of

! the (n) new roots.

1143

1149

! the (n) new roots.

!

1144

1150

!

1145

1151

!

! remarks muller always returns the last approximation for root j

1146

1152

! remarks muller always returns the last approximation for root j

! in x(j). if the convergence criterion is satisfied,

1147

1153

! in x(j). if the convergence criterion is satisfied,

! then infer(j) is less than or equal to itmax. if the

1148

1154

! then infer(j) is less than or equal to itmax. if the

! convergence criterion is not satisified, then infer(j)

1149

1155

! convergence criterion is not satisified, then infer(j)

! is set to either itmax+1 or itmax+k, with k greater

1150

1156

! is set to either itmax+1 or itmax+k, with k greater

! than 1. infer(j) = itmax+1 indicates that muller did

1151

1157

! than 1. infer(j) = itmax+1 indicates that muller did

! not obtain convergence in the allowed number of iter-

1152

1158

! not obtain convergence in the allowed number of iter-

! ations. in this case, the user may wish to set itmax

1153

1159

! ations. in this case, the user may wish to set itmax

! to a larger value. infer(j) = itmax+k means that con-

1154

1160

! to a larger value. infer(j) = itmax+k means that con-

! vergence was obtained (on iteration k) for the defla-

1155

1161

! vergence was obtained (on iteration k) for the defla-

! ted function

1156

1162

! ted function

! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))

1157

1163

! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))

!

1158

1164

!

! but failed for f(z). in this case, better initial

1159

1165

! but failed for f(z). in this case, better initial

! guesses might help or, it might be necessary to relax

1160

1166

! guesses might help or, it might be necessary to relax

! the convergence criterion.

1161

1167

! the convergence criterion.

!

1162

1168

!

!-----------------------------------------------------------------------

1163

1169

!-----------------------------------------------------------------------

!

1164

1170

!

subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)

1165

1171

subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)

implicit none

1166

1172

implicit none

double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq

1167

1173

double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq

double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &

1168

1174

double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &

tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &

1169

1175

tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &

zero,p1,one,four,p5

1170

1176

zero,p1,one,four,p5

1171

1177

double complex, external :: f

1172

1178

double complex, external :: f

integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &

1173

1179

integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &

knpng,jk,ick,nn,lm1,errcode

1174

1180

knpng,jk,ick,nn,lm1,errcode

double complex :: x(kn+n)

1175

1181

double complex :: x(kn+n)

integer :: infer(kn+n)

1176

1182

integer :: infer(kn+n)

1177

1183

1178

1184

data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &

1179

1185

data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &

one/(1.0,0.0)/,four/(4.0,0.0)/, &

1180

1186

one/(1.0,0.0)/,four/(4.0,0.0)/, &

p5/(0.5,0.0)/, &

1181

1187

p5/(0.5,0.0)/, &

rzero/0.0/,rten/10.0/,rhun/100.0/, &

1182

1188

rzero/0.0/,rten/10.0/,rhun/100.0/, &

ax/0.1/,ickmax/3/,rp01/0.01/

1183

1189

ax/0.1/,ickmax/3/,rp01/0.01/

1184

1190

ier = 0

1185

1191

ier = 0

if (n .lt. 1) then ! What the hell are doing here then ...

1186

1192

if (n .lt. 1) then ! What the hell are doing here then ...

return

1187

1193

return

end if

1188

1194

end if

!eps1 = rten **(-nsig)

1189

1195

!eps1 = rten **(-nsig)

eps1 = min(eps1,rp01)

1190

1196

eps1 = min(eps1,rp01)

1191

1197

knp1 = kn+1

1192

1198

knp1 = kn+1

knpn = kn+n

1193

1199

knpn = kn+n

knpng = kn+nguess

1194

1200

knpng = kn+nguess

do i=1,knpn

1195

1201

do i=1,knpn

infer(i) = 0

1196

1202

infer(i) = 0

if (i .gt. knpng) x(i) = zero

1197

1203

if (i .gt. knpng) x(i) = zero

end do

1198

1204

end do

l= knp1

1199

1205

l= knp1

1200

1206

ic=0

1201

1207

ic=0

rloop: do while (l<=knpn) ! Main loop over new roots

1202

1208

rloop: do while (l<=knpn) ! Main loop over new roots

jk = 0

1203

1209

jk = 0

ick = 0

1204

1210

ick = 0

xl = x(l)

1205

1211

xl = x(l)

icloop: do

1206

1212

icloop: do

ic = 0

1207

1213

ic = 0

h = ax

1208

1214

h = ax

h = p1*h

1209

1215

h = p1*h

if (cdabs(xl) .gt. ax) h = p1*xl

1210

1216

if (cdabs(xl) .gt. ax) h = p1*xl

! first three points are

1211

1217

! first three points are

! xl+h, xl-h, xl

1212

1218

! xl+h, xl-h, xl

rt = xl+h

1213

1219

rt = xl+h

call deflated_work(errcode)

1214

1220

call deflated_work(errcode)

if (errcode == 1) then

1215

1221

if (errcode == 1) then

exit icloop

1216

1222

exit icloop

end if

1217

1223

end if

1218

1224

z0 = fprt

1219

1225

z0 = fprt

y0 = frt

1220

1226

y0 = frt

x0 = rt

1221

1227

x0 = rt

rt = xl-h

1222

1228

rt = xl-h

call deflated_work(errcode)

1223

1229

call deflated_work(errcode)

if (errcode == 1) then

1224

1230

if (errcode == 1) then

exit icloop

1225

1231

exit icloop

end if

1226

1232

end if

1227

1233

z1 = fprt

1228

1234

z1 = fprt

y1 = frt

1229

1235

y1 = frt

h = xl-rt

1230

1236

h = xl-rt

d = h/(rt-x0)

1231

1237

d = h/(rt-x0)

rt = xl

1232

1238

rt = xl

1233

1239

call deflated_work(errcode)

1234

1240

call deflated_work(errcode)

if (errcode == 1) then

1235

1241

if (errcode == 1) then

exit icloop

1236

1242

exit icloop

end if

1237

1243

end if

1238

1244

1239

1245

z2 = fprt

1240

1246

z2 = fprt

y2 = frt

1241

1247

y2 = frt

! begin main algorithm

1242

1248

! begin main algorithm

iloop: do

1243

1249

iloop: do

dd = one + d

1244

1250

dd = one + d

t1 = z0*d*d

1245

1251

t1 = z0*d*d

t2 = z1*dd*dd

1246

1252

t2 = z1*dd*dd

xx = z2*dd

1247

1253

xx = z2*dd

t3 = z2*d

1248

1254

t3 = z2*d

bi = t1-t2+xx+t3

1249

1255

bi = t1-t2+xx+t3

den = bi*bi-four*(xx*t1-t3*(t2-xx))

1250

1256

den = bi*bi-four*(xx*t1-t3*(t2-xx))

! use denominator of maximum amplitude

1251

1257

! use denominator of maximum amplitude

t1 = cdsqrt(den)

1252

1258

t1 = cdsqrt(den)

qz = rhun*max(cdabs(bi),cdabs(t1))

1253

1259

qz = rhun*max(cdabs(bi),cdabs(t1))

t2 = bi + t1

1254

1260

t2 = bi + t1

tpq = cdabs(t2)+qz

1255

1261

tpq = cdabs(t2)+qz

if (tpq .eq. qz) t2 = zero

1256

1262

if (tpq .eq. qz) t2 = zero

t3 = bi - t1

1257

1263

t3 = bi - t1

tpq = cdabs(t3) + qz

1258

1264

tpq = cdabs(t3) + qz

if (tpq .eq. qz) t3 = zero

1259

1265

if (tpq .eq. qz) t3 = zero

den = t2

1260

1266

den = t2

qz = cdabs(t3)-cdabs(t2)

1261

1267

qz = cdabs(t3)-cdabs(t2)

if (qz .gt. rzero) den = t3

1262

1268

if (qz .gt. rzero) den = t3

! test for zero denominator

1263

1269

! test for zero denominator

if (cdabs(den) .eq. rzero) then

1264

1270

if (cdabs(den) .eq. rzero) then

call trans_rt()

1265

1271

call trans_rt()

call deflated_work(errcode)

1266

1272

call deflated_work(errcode)

if (errcode == 1) then

1267

1273

if (errcode == 1) then

exit icloop

1268

1274

exit icloop

end if

1269

1275

end if

z2 = fprt

1270

1276

z2 = fprt

y2 = frt

1271

1277

y2 = frt

cycle iloop

1272

1278

cycle iloop

end if

1273

1279

end if

1274

1280

1275

1281

d = -xx/den

1276

1282

d = -xx/den

d = d+d

1277

1283

d = d+d

h = d*h

1278

1284

h = d*h

rt = rt + h

1279

1285

rt = rt + h

! check convergence of the first kind

1280

1286

! check convergence of the first kind

if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then

1281

1287

if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then

if (ic .ne. 0) then

1282

1288

if (ic .ne. 0) then

exit icloop

1283

1289

exit icloop

end if

1284

1290

end if

ic = 1

1285

1291

ic = 1

z0 = y1

1286

1292

z0 = y1

z1 = y2

1287

1293

z1 = y2

z2 = f(rt)

1288

1294

z2 = f(rt)

xl = rt

1289

1295

xl = rt

ick = ick+1

1290

1296

ick = ick+1

if (ick .le. ickmax) then

1291

1297

if (ick .le. ickmax) then

cycle iloop

1292

1298

cycle iloop

end if

1293

1299

end if

! warning error, itmax = maximum

1294

1300

! warning error, itmax = maximum

jk = itmax + jk

1295

1301

jk = itmax + jk

ier = 33

1296

1302

ier = 33

end if

1297

1303

end if

if (ic .ne. 0) then

1298

1304

if (ic .ne. 0) then

cycle icloop

1299

1305

cycle icloop

end if

1300

1306

end if

call deflated_work(errcode)

1301

1307

call deflated_work(errcode)

if (errcode == 1) then

1302

1308

if (errcode == 1) then

exit icloop

1303

1309

exit icloop

end if

1304

1310

end if

1305

1311

do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)

1306

1312

do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)

! take remedial action to induce

1307

1313

! take remedial action to induce

! convergence

1308

1314

! convergence

d = d*p5

1309

1315

d = d*p5

h = h*p5

1310

1316

h = h*p5

rt = rt-h

1311

1317

rt = rt-h

call deflated_work(errcode)

1312

1318

call deflated_work(errcode)

if (errcode == 1) then

1313

1319

if (errcode == 1) then

exit icloop

1314

1320

exit icloop

end if

1315

1321

end if

end do

1316

1322

end do

z0 = z1

1317

1323

z0 = z1

z1 = z2

1318

1324

z1 = z2

z2 = fprt

1319

1325

z2 = fprt

y0 = y1

1320

1326

y0 = y1

y1 = y2

1321

1327

y1 = y2

y2 = frt

1322

1328

y2 = frt

end do iloop

1323

1329

end do iloop

end do icloop

1324

1330

end do icloop

x(l) = rt

1325

1331

x(l) = rt

infer(l) = jk

1326

1332

infer(l) = jk

l = l+1

1327

1333

l = l+1

end do rloop

1328

1334

end do rloop

1329

1335

contains

1330

1336

contains

subroutine trans_rt()

1331

1337

subroutine trans_rt()

tem = rten*eps1

1332

1338

tem = rten*eps1

if (cdabs(rt) .gt. ax) tem = tem*rt

1333

1339

if (cdabs(rt) .gt. ax) tem = tem*rt

rt = rt+tem

1334

1340

rt = rt+tem

d = (h+tem)*d/h

1335

1341

d = (h+tem)*d/h

h = h+tem

1336

1342

h = h+tem

end subroutine trans_rt

1337

1343

end subroutine trans_rt

1338

1344

subroutine deflated_work(errcode)

1339

1345

subroutine deflated_work(errcode)

! errcode=0 => no errors

1340

1346

! errcode=0 => no errors

! errcode=1 => jk>itmax or convergence of second kind achieved

1341

1347

! errcode=1 => jk>itmax or convergence of second kind achieved

integer :: errcode,flag

1342

1348

integer :: errcode,flag

1343

1349

flag=1

1344

1350

flag=1

loop1: do while(flag==1)

1345

1351

loop1: do while(flag==1)

errcode=0

1346

1352

errcode=0

jk = jk+1

1347

1353

jk = jk+1

if (jk .gt. itmax) then

1348

1354

if (jk .gt. itmax) then

ier=33

1349

1355

ier=33

errcode=1

1350

1356

errcode=1

return

1351

1357

return

end if

1352

1358

end if

frt = f(rt)

1353

1359

frt = f(rt)

fprt = frt

1354

1360

fprt = frt

if (l /= 1) then

1355

1361

if (l /= 1) then

lm1 = l-1

1356

1362

lm1 = l-1

do i=1,lm1

1357

1363

do i=1,lm1

tem = rt - x(i)

1358

1364

tem = rt - x(i)

if (cdabs(tem) .eq. rzero) then

1359

1365

if (cdabs(tem) .eq. rzero) then

!if (ic .ne. 0) go to 15 !! ?? possible?

1360

1366

!if (ic .ne. 0) go to 15 !! ?? possible?

call trans_rt()

1361

1367

call trans_rt()

cycle loop1

1362

1368

cycle loop1

end if

1363

1369

end if

fprt = fprt/tem

1364

1370

fprt = fprt/tem

end do

1365

1371

end do

end if

1366

1372

end if

flag=0

1367

1373

flag=0

end do loop1

1368

1374

end do loop1

1369

1375

if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then

1370

1376

if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then

errcode=1

1371

1377

errcode=1

return

1372

1378

return

end if

1373

1379

end if

1374

1380

end subroutine deflated_work

1375

1381

end subroutine deflated_work

1376

1382

end subroutine

1377

1383

end subroutine

1378

1384

1379

1385

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

1380

1386

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

!

1381

1387

!

! Integration

1382

1388

! Integration

!

1383

1389

!

! Only double precision coded atm

1384

1390

! Only double precision coded atm

!

1385

1391

!

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

1386

1392

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

1387

1393

1388

1394

subroutine fvn_d_gauss_legendre(n,qx,qw)

1389

1395

subroutine fvn_d_gauss_legendre(n,qx,qw)

!

1390

1396

!

! This routine compute the n Gauss Legendre abscissas and weights

1391

1397

! This routine compute the n Gauss Legendre abscissas and weights

! Adapted from Numerical Recipes routine gauleg

1392

1398

! Adapted from Numerical Recipes routine gauleg

!

1393

1399

!

! n (in) : number of points

1394

1400

! n (in) : number of points

! qx(out) : abscissas

1395

1401

! qx(out) : abscissas

! qw(out) : weights

1396

1402

! qw(out) : weights

!

1397

1403

!

implicit none

1398

1404

implicit none

double precision,parameter :: pi=3.141592653589793d0

1399

1405

double precision,parameter :: pi=3.141592653589793d0

integer, intent(in) :: n

1400

1406

integer, intent(in) :: n

double precision, intent(out) :: qx(n),qw(n)

1401

1407

double precision, intent(out) :: qx(n),qw(n)

1402

1408

integer :: m,i,j

1403

1409

integer :: m,i,j

double precision :: z,z1,p1,p2,p3,pp

1404

1410

double precision :: z,z1,p1,p2,p3,pp

m=(n+1)/2

1405

1411

m=(n+1)/2

do i=1,m

1406

1412

do i=1,m

z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))

1407

1413

z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))

iloop: do

1408

1414

iloop: do

p1=1.d0

1409

1415

p1=1.d0

p2=0.d0

1410

1416

p2=0.d0

do j=1,n

1411

1417

do j=1,n

p3=p2

1412

1418

p3=p2

p2=p1

1413

1419

p2=p1

p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)

1414

1420

p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)

end do

1415

1421

end do

pp=dble(n)*(z*p1-p2)/(z*z-1.d0)

1416

1422

pp=dble(n)*(z*p1-p2)/(z*z-1.d0)

z1=z

1417

1423

z1=z

z=z1-p1/pp

1418

1424

z=z1-p1/pp

if (dabs(z-z1)<=epsilon(z)) then

1419

1425

if (dabs(z-z1)<=epsilon(z)) then

exit iloop

1420

1426

exit iloop

end if

1421

1427

end if

end do iloop

1422

1428

end do iloop

qx(i)=-z

1423

1429

qx(i)=-z

qx(n+1-i)=z

1424

1430

qx(n+1-i)=z

qw(i)=2.d0/((1.d0-z*z)*pp*pp)

1425

1431

qw(i)=2.d0/((1.d0-z*z)*pp*pp)

qw(n+1-i)=qw(i)

1426

1432

qw(n+1-i)=qw(i)

end do

1427

1433

end do

end subroutine

1428

1434

end subroutine

1429

1435

1430

1436

1431

1437

subroutine fvn_d_gl_integ(f,a,b,n,res)

1432

1438

subroutine fvn_d_gl_integ(f,a,b,n,res)

!

1433

1439

!

! This is a simple non adaptative integration routine

1434

1440

! This is a simple non adaptative integration routine

! using n gauss legendre abscissas and weights

1435

1441

! using n gauss legendre abscissas and weights

!

1436

1442

!

! f(in) : the function to integrate

1437

1443

! f(in) : the function to integrate

! a(in) : lower bound

1438

1444

! a(in) : lower bound

! b(in) : higher bound

1439

1445

! b(in) : higher bound

! n(in) : number of gauss legendre pairs

1440

1446

! n(in) : number of gauss legendre pairs

! res(out): the evaluation of the integral

1441

1447

! res(out): the evaluation of the integral

!

1442

1448

!

double precision,external :: f

1443

1449

double precision,external :: f

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

1444

1450

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

integer, intent(in):: n

1445

1451

integer, intent(in):: n

double precision, intent(out) :: res

1446

1452

double precision, intent(out) :: res

1447

1453

double precision, allocatable :: qx(:),qw(:)

1448

1454

double precision, allocatable :: qx(:),qw(:)

double precision :: xm,xr

1449

1455

double precision :: xm,xr

integer :: i

1450

1456

integer :: i

1451

1457

! First compute n gauss legendre abs and weight

1452

1458

! First compute n gauss legendre abs and weight

allocate(qx(n))

1453

1459

allocate(qx(n))

allocate(qw(n))

1454

1460

allocate(qw(n))

call fvn_d_gauss_legendre(n,qx,qw)

1455

1461

call fvn_d_gauss_legendre(n,qx,qw)

1456

1462

xm=0.5d0*(b+a)

1457

1463

xm=0.5d0*(b+a)

xr=0.5d0*(b-a)

1458

1464

xr=0.5d0*(b-a)

1459

1465

res=0.d0

1460

1466

res=0.d0

1461

1467

do i=1,n

1462

1468

do i=1,n

res=res+qw(i)*f(xm+xr*qx(i))

1463

1469

res=res+qw(i)*f(xm+xr*qx(i))

end do

1464

1470

end do

1465

1471

res=xr*res

1466

1472

res=xr*res

1467

1473

deallocate(qw)

1468

1474

deallocate(qw)

deallocate(qx)

1469

1475

deallocate(qx)

1470

1476

end subroutine

1471

1477

end subroutine

1472

1478

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

1473

1479

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

!

1474

1480

!

! Simple and double adaptative Gauss Kronrod integration based on

1475

1481

! Simple and double adaptative Gauss Kronrod integration based on

! a modified version of quadpack ( http://www.netlib.org/quadpack

1476

1482

! a modified version of quadpack ( http://www.netlib.org/quadpack

!

1477

1483

!

! Common parameters :

1478

1484

! Common parameters :

!

1479

1485

!

! key (in)

1480

1486

! key (in)

! epsabs

1481

1487

! epsabs

! epsrel

1482

1488

! epsrel

!

1483

1489

!

1484

1490

!

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

1485

1491

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

1486

1492

subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

1487

1493

subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

!

1488

1494

!

! Evaluate the integral of function f(x) between a and b

1489

1495

! Evaluate the integral of function f(x) between a and b

!

1490

1496

!

! f(in) : the function

1491

1497

! f(in) : the function

! a(in) : lower bound

1492

1498

! a(in) : lower bound

! b(in) : higher bound

1493

1499

! b(in) : higher bound

! epsabs(in) : desired absolute error

1494

1500

! epsabs(in) : desired absolute error

! epsrel(in) : desired relative error

1495

1501

! epsrel(in) : desired relative error

! key(in) : gauss kronrod rule

1496

1502

! key(in) : gauss kronrod rule

! 1: 7 - 15 points

1497

1503

! 1: 7 - 15 points

! 2: 10 - 21 points

1498

1504

! 2: 10 - 21 points

! 3: 15 - 31 points

1499

1505

! 3: 15 - 31 points

! 4: 20 - 41 points

1500

1506

! 4: 20 - 41 points

! 5: 25 - 51 points

1501

1507

! 5: 25 - 51 points

! 6: 30 - 61 points

1502

1508

! 6: 30 - 61 points

!

1503

1509

!

! limit(in) : maximum number of subintervals in the partition of the

1504

1510

! limit(in) : maximum number of subintervals in the partition of the

! given integration interval (a,b). A value of 500 will give the same

1505

1511

! given integration interval (a,b). A value of 500 will give the same

! behaviour as the imsl routine dqdag

1506

1512

! behaviour as the imsl routine dqdag

!

1507

1513

!

! res(out) : estimated integral value

1508

1514

! res(out) : estimated integral value

! abserr(out) : estimated absolute error

1509

1515

! abserr(out) : estimated absolute error

! ier(out) : error flag from quadpack routines

1510

1516

! ier(out) : error flag from quadpack routines

! 0 : no error

1511

1517

! 0 : no error

! 1 : maximum number of subdivisions allowed

1512

1518

! 1 : maximum number of subdivisions allowed

! has been achieved. one can allow more

1513

1519

! has been achieved. one can allow more

! subdivisions by increasing the value of

1514

1520

! subdivisions by increasing the value of

! limit (and taking the according dimension

1515

1521

! limit (and taking the according dimension

! adjustments into account). however, if

1516

1522

! adjustments into account). however, if

! this yield no improvement it is advised

1517

1523

! this yield no improvement it is advised

! to analyze the integrand in order to

1518

1524

! to analyze the integrand in order to

! determine the integration difficulaties.

1519

1525

! determine the integration difficulaties.

! if the position of a local difficulty can

1520

1526

! if the position of a local difficulty can

! be determined (i.e.singularity,

1521

1527

! be determined (i.e.singularity,

! discontinuity within the interval) one

1522

1528

! discontinuity within the interval) one

! will probably gain from splitting up the

1523

1529

! will probably gain from splitting up the

! interval at this point and calling the

1524

1530

! interval at this point and calling the

! integrator on the subranges. if possible,

1525

1531

! integrator on the subranges. if possible,

! an appropriate special-purpose integrator

1526

1532

! an appropriate special-purpose integrator

! should be used which is designed for

1527

1533

! should be used which is designed for

! handling the type of difficulty involved.

1528

1534

! handling the type of difficulty involved.

! 2 : the occurrence of roundoff error is

1529

1535

! 2 : the occurrence of roundoff error is

! detected, which prevents the requested

1530

1536

! detected, which prevents the requested

! tolerance from being achieved.

1531

1537

! tolerance from being achieved.

! 3 : extremely bad integrand behaviour occurs

1532

1538

! 3 : extremely bad integrand behaviour occurs

! at some points of the integration

1533

1539

! at some points of the integration

! interval.

1534

1540

! interval.

! 6 : the input is invalid, because

1535

1541

! 6 : the input is invalid, because

! (epsabs.le.0 and

1536

1542

! (epsabs.le.0 and

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

1537

1543

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

! or limit.lt.1 or lenw.lt.limit*4.

1538

1544

! or limit.lt.1 or lenw.lt.limit*4.

! result, abserr, neval, last are set

1539

1545

! result, abserr, neval, last are set

! to zero.

1540

1546

! to zero.

! except when lenw is invalid, iwork(1),

1541

1547

! except when lenw is invalid, iwork(1),

! work(limit*2+1) and work(limit*3+1) are

1542

1548

! work(limit*2+1) and work(limit*3+1) are

! set to zero, work(1) is set to a and

1543

1549

! set to zero, work(1) is set to a and

! work(limit+1) to b.

1544

1550

! work(limit+1) to b.

1545

1551

implicit none

1546

1552

implicit none

double precision, external :: f

1547

1553

double precision, external :: f

double precision, intent(in) :: a,b,epsabs,epsrel

1548

1554

double precision, intent(in) :: a,b,epsabs,epsrel

integer, intent(in) :: key

1549

1555

integer, intent(in) :: key

integer, intent(in) :: limit

1550

1556

integer, intent(in) :: limit

double precision, intent(out) :: res,abserr

1551

1557

double precision, intent(out) :: res,abserr

integer, intent(out) :: ier

1552

1558

integer, intent(out) :: ier

1553

1559

double precision, allocatable :: work(:)

1554

1560

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

1555

1561

integer, allocatable :: iwork(:)

integer :: lenw,neval,last

1556

1562

integer :: lenw,neval,last

1557

1563

! imsl value for limit is 500

1558

1564

! imsl value for limit is 500

lenw=limit*4

1559

1565

lenw=limit*4

1560

1566

allocate(iwork(limit))

1561

1567

allocate(iwork(limit))

allocate(work(lenw))

1562

1568

allocate(work(lenw))

1563

1569

call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

1564

1570

call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

1565

1571

deallocate(work)

1566

1572

deallocate(work)

deallocate(iwork)

1567

1573

deallocate(iwork)

1568

1574

end subroutine

1569

1575

end subroutine

1570

1576

1571

1577

1572

1578

subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)

1573

1579

subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)

!

1574

1580

!

! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)

1575

1581

! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)

!

1576

1582

!

! f(in) : the function

1577

1583

! f(in) : the function

! a(in) : lower bound

1578

1584

! a(in) : lower bound

! b(in) : higher bound

1579

1585

! b(in) : higher bound

! g(in) : external function describing lower bound for y

1580

1586

! g(in) : external function describing lower bound for y

! h(in) : external function describing higher bound for y

1581

1587

! h(in) : external function describing higher bound for y

! epsabs(in) : desired absolute error

1582

1588

! epsabs(in) : desired absolute error

! epsrel(in) : desired relative error

1583

1589

! epsrel(in) : desired relative error

! key(in) : gauss kronrod rule

1584

1590

! key(in) : gauss kronrod rule

! 1: 7 - 15 points

1585

1591

! 1: 7 - 15 points

! 2: 10 - 21 points

1586

1592

! 2: 10 - 21 points

! 3: 15 - 31 points

1587

1593

! 3: 15 - 31 points

! 4: 20 - 41 points

1588

1594

! 4: 20 - 41 points

! 5: 25 - 51 points

1589

1595

! 5: 25 - 51 points

! 6: 30 - 61 points

1590

1596

! 6: 30 - 61 points

!

1591

1597

!

! limit(in) : maximum number of subintervals in the partition of the

1592

1598

! limit(in) : maximum number of subintervals in the partition of the

! given integration interval (a,b). A value of 500 will give the same

1593

1599

! given integration interval (a,b). A value of 500 will give the same

! behaviour as the imsl routine dqdag

1594

1600

! behaviour as the imsl routine dqdag

!

1595

1601

!

! res(out) : estimated integral value

1596

1602

! res(out) : estimated integral value

! abserr(out) : estimated absolute error

1597

1603

! abserr(out) : estimated absolute error

! ier(out) : error flag from quadpack routines

1598

1604

! ier(out) : error flag from quadpack routines

! 0 : no error

1599

1605

! 0 : no error

! 1 : maximum number of subdivisions allowed

1600

1606

! 1 : maximum number of subdivisions allowed

! has been achieved. one can allow more

1601

1607

! has been achieved. one can allow more

! subdivisions by increasing the value of

1602

1608

! subdivisions by increasing the value of

! limit (and taking the according dimension

1603

1609

! limit (and taking the according dimension

! adjustments into account). however, if

1604

1610

! adjustments into account). however, if

! this yield no improvement it is advised

1605

1611

! this yield no improvement it is advised

! to analyze the integrand in order to

1606

1612

! to analyze the integrand in order to

! determine the integration difficulaties.

1607

1613

! determine the integration difficulaties.

! if the position of a local difficulty can

1608

1614

! if the position of a local difficulty can

! be determined (i.e.singularity,

1609

1615

! be determined (i.e.singularity,

! discontinuity within the interval) one

1610

1616

! discontinuity within the interval) one

! will probably gain from splitting up the

1611

1617

! will probably gain from splitting up the

! interval at this point and calling the

1612

1618

! interval at this point and calling the

! integrator on the subranges. if possible,

1613

1619

! integrator on the subranges. if possible,

! an appropriate special-purpose integrator

1614

1620

! an appropriate special-purpose integrator

! should be used which is designed for

1615

1621

! should be used which is designed for

! handling the type of difficulty involved.

1616

1622

! handling the type of difficulty involved.

! 2 : the occurrence of roundoff error is

1617

1623

! 2 : the occurrence of roundoff error is

! detected, which prevents the requested

1618

1624

! detected, which prevents the requested

! tolerance from being achieved.

1619

1625

! tolerance from being achieved.

! 3 : extremely bad integrand behaviour occurs

1620

1626

! 3 : extremely bad integrand behaviour occurs

! at some points of the integration

1621

1627

! at some points of the integration

! interval.

1622

1628

! interval.

! 6 : the input is invalid, because

1623

1629

! 6 : the input is invalid, because

! (epsabs.le.0 and

1624

1630

! (epsabs.le.0 and

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

1625

1631

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

! or limit.lt.1 or lenw.lt.limit*4.

1626

1632

! or limit.lt.1 or lenw.lt.limit*4.

! result, abserr, neval, last are set

1627

1633

! result, abserr, neval, last are set

! to zero.

1628

1634

! to zero.

! except when lenw is invalid, iwork(1),

1629

1635

! except when lenw is invalid, iwork(1),

! work(limit*2+1) and work(limit*3+1) are

1630

1636

! work(limit*2+1) and work(limit*3+1) are

! set to zero, work(1) is set to a and

1631

1637

! set to zero, work(1) is set to a and

! work(limit+1) to b.

1632

1638

! work(limit+1) to b.

1633

1639

implicit none

1634

1640

implicit none

double precision, external:: f,g,h

1635

1641

double precision, external:: f,g,h

double precision, intent(in) :: a,b,epsabs,epsrel

1636

1642

double precision, intent(in) :: a,b,epsabs,epsrel

integer, intent(in) :: key,limit

1637

1643

integer, intent(in) :: key,limit

integer, intent(out) :: ier

1638

1644

integer, intent(out) :: ier

double precision, intent(out) :: res,abserr

1639

1645

double precision, intent(out) :: res,abserr

1640

1646

1641

1647

double precision, allocatable :: work(:)

1642

1648

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

1643

1649

integer, allocatable :: iwork(:)

integer :: lenw,neval,last

1644

1650

integer :: lenw,neval,last

1645

1651

! imsl value for limit is 500

1646

1652

! imsl value for limit is 500

lenw=limit*4

1647

1653

lenw=limit*4

allocate(work(lenw))

1648

1654

allocate(work(lenw))

allocate(iwork(limit))

1649

1655

allocate(iwork(limit))

1650

1656

call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

1651

1657

call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

1652

1658

deallocate(iwork)

1653

1659

deallocate(iwork)

deallocate(work)

1654

1660

deallocate(work)

end subroutine

1655

1661

end subroutine

1656

1662

1657

1663

1658

1664

subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

1659

1665

subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)

!

1660

1666

!

! Evaluate the single integral of function f(x,y) for y between a and b with a

1661

1667

! Evaluate the single integral of function f(x,y) for y between a and b with a

! given x value

1662

1668

! given x value

!

1663

1669

!

! This function is used for the evaluation of the double integral fvn_d_integ_2_gk

1664

1670

! This function is used for the evaluation of the double integral fvn_d_integ_2_gk

!

1665

1671

!

! f(in) : the function

1666

1672

! f(in) : the function

! x(in) : x

1667

1673

! x(in) : x

! a(in) : lower bound

1668

1674

! a(in) : lower bound

! b(in) : higher bound

1669

1675

! b(in) : higher bound

! epsabs(in) : desired absolute error

1670

1676

! epsabs(in) : desired absolute error

! epsrel(in) : desired relative error

1671

1677

! epsrel(in) : desired relative error

! key(in) : gauss kronrod rule

1672

1678

! key(in) : gauss kronrod rule

! 1: 7 - 15 points

1673

1679

! 1: 7 - 15 points

! 2: 10 - 21 points

1674

1680

! 2: 10 - 21 points

! 3: 15 - 31 points

1675

1681

! 3: 15 - 31 points

! 4: 20 - 41 points

1676

1682

! 4: 20 - 41 points

! 5: 25 - 51 points

1677

1683

! 5: 25 - 51 points

! 6: 30 - 61 points

1678

1684

! 6: 30 - 61 points

!

1679

1685

!

! limit(in) : maximum number of subintervals in the partition of the

1680

1686

! limit(in) : maximum number of subintervals in the partition of the

! given integration interval (a,b). A value of 500 will give the same

1681

1687

! given integration interval (a,b). A value of 500 will give the same

! behaviour as the imsl routine dqdag

1682

1688

! behaviour as the imsl routine dqdag

!

1683

1689

!

! res(out) : estimated integral value

1684

1690

! res(out) : estimated integral value

! abserr(out) : estimated absolute error

1685

1691

! abserr(out) : estimated absolute error

! ier(out) : error flag from quadpack routines

1686

1692

! ier(out) : error flag from quadpack routines

! 0 : no error

1687

1693

! 0 : no error

! 1 : maximum number of subdivisions allowed

1688

1694

! 1 : maximum number of subdivisions allowed

! has been achieved. one can allow more

1689

1695

! has been achieved. one can allow more

! subdivisions by increasing the value of

1690

1696

! subdivisions by increasing the value of

! limit (and taking the according dimension

1691

1697

! limit (and taking the according dimension

! adjustments into account). however, if

1692

1698

! adjustments into account). however, if

! this yield no improvement it is advised

1693

1699

! this yield no improvement it is advised

! to analyze the integrand in order to

1694

1700

! to analyze the integrand in order to

! determine the integration difficulaties.

1695

1701

! determine the integration difficulaties.

! if the position of a local difficulty can

1696

1702

! if the position of a local difficulty can

! be determined (i.e.singularity,

1697

1703

! be determined (i.e.singularity,

! discontinuity within the interval) one

1698

1704

! discontinuity within the interval) one

! will probably gain from splitting up the

1699

1705

! will probably gain from splitting up the

! interval at this point and calling the

1700

1706

! interval at this point and calling the

! integrator on the subranges. if possible,

1701

1707

! integrator on the subranges. if possible,

! an appropriate special-purpose integrator

1702

1708

! an appropriate special-purpose integrator

! should be used which is designed for

1703

1709

! should be used which is designed for

! handling the type of difficulty involved.

1704

1710

! handling the type of difficulty involved.

! 2 : the occurrence of roundoff error is

1705

1711

! 2 : the occurrence of roundoff error is

! detected, which prevents the requested

1706

1712

! detected, which prevents the requested

! tolerance from being achieved.

1707

1713

! tolerance from being achieved.

! 3 : extremely bad integrand behaviour occurs

1708

1714

! 3 : extremely bad integrand behaviour occurs

! at some points of the integration

1709

1715

! at some points of the integration

! interval.

1710

1716

! interval.

! 6 : the input is invalid, because

1711

1717

! 6 : the input is invalid, because

! (epsabs.le.0 and

1712

1718

! (epsabs.le.0 and

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

1713

1719

! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))

! or limit.lt.1 or lenw.lt.limit*4.

1714

1720

! or limit.lt.1 or lenw.lt.limit*4.

! result, abserr, neval, last are set

1715

1721

! result, abserr, neval, last are set

! to zero.

1716

1722

! to zero.

! except when lenw is invalid, iwork(1),

1717

1723

! except when lenw is invalid, iwork(1),

! work(limit*2+1) and work(limit*3+1) are

1718

1724

! work(limit*2+1) and work(limit*3+1) are

! set to zero, work(1) is set to a and

1719

1725

! set to zero, work(1) is set to a and

! work(limit+1) to b.

1720

1726

! work(limit+1) to b.

1721

1727

implicit none

1722

1728

implicit none

double precision, external:: f

1723

1729

double precision, external:: f

double precision, intent(in) :: x,a,b,epsabs,epsrel

1724

1730

double precision, intent(in) :: x,a,b,epsabs,epsrel

integer, intent(in) :: key,limit

1725

1731

integer, intent(in) :: key,limit

integer, intent(out) :: ier

1726

1732

integer, intent(out) :: ier

double precision, intent(out) :: res,abserr

1727

1733

double precision, intent(out) :: res,abserr

1728

1734

1729

1735

double precision, allocatable :: work(:)

1730

1736

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

1731

1737

integer, allocatable :: iwork(:)

integer :: lenw,neval,last

1732

1738

integer :: lenw,neval,last

1733

1739

! imsl value for limit is 500

1734

1740

! imsl value for limit is 500

lenw=limit*4

1735

1741

lenw=limit*4

allocate(work(lenw))

1736

1742

allocate(work(lenw))

allocate(iwork(limit))

1737

1743

allocate(iwork(limit))

1738

1744

call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

1739

1745

call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)

1740

1746

deallocate(iwork)

1741

1747

deallocate(iwork)

deallocate(work)

1742

1748

deallocate(work)

end subroutine

1743

1749

end subroutine

1744

1750

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

1745

1751

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

! Include the modified quadpack files

1746

1752

! Include the modified quadpack files

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

1747

1753

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

include "fvn_quadpack/dqag_2d_inner.f"

1748

1754

include "fvn_quadpack/dqag_2d_inner.f"

include "fvn_quadpack/dqk15_2d_inner.f"

1749

1755

include "fvn_quadpack/dqk15_2d_inner.f"

include "fvn_quadpack/dqk31_2d_outer.f"

1750

1756

include "fvn_quadpack/dqk31_2d_outer.f"

include "fvn_quadpack/d1mach.f"

1751

1757

include "fvn_quadpack/d1mach.f"

include "fvn_quadpack/dqk31_2d_inner.f"

1752

1758

include "fvn_quadpack/dqk31_2d_inner.f"

include "fvn_quadpack/dqage.f"

1753

1759

include "fvn_quadpack/dqage.f"

include "fvn_quadpack/dqk15.f"

1754

1760

include "fvn_quadpack/dqk15.f"

include "fvn_quadpack/dqk21.f"

1755

1761

include "fvn_quadpack/dqk21.f"

include "fvn_quadpack/dqk31.f"

1756

1762

include "fvn_quadpack/dqk31.f"

include "fvn_quadpack/dqk41.f"

1757

1763

include "fvn_quadpack/dqk41.f"

include "fvn_quadpack/dqk51.f"

1758

1764

include "fvn_quadpack/dqk51.f"

include "fvn_quadpack/dqk61.f"

1759

1765

include "fvn_quadpack/dqk61.f"

include "fvn_quadpack/dqk41_2d_outer.f"

1760

1766

include "fvn_quadpack/dqk41_2d_outer.f"

GITLAB

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