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

!

! svn version

25

! svn version

! June 2007 : added some complex trigonometric functions

26

! June 2007 : added some complex trigonometric functions

!

27

!

! TO DO LIST :

28

! TO DO LIST :

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

29

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

! eigenvalues are given with no particular order.

30

! eigenvalues are given with no particular order.

! + Generic interface for fvn_x_name family -> fvn_name

31

! + Generic interface for fvn_x_name family -> fvn_name

! + Make some parameters optional, status for example

32

! + Make some parameters optional, status for example

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

33

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

! + unify quadpack routines

34

! + unify quadpack routines

! + ...

35

! + ...

!

36

!

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

37

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

38

implicit none

39

implicit none

! We define pi and i for the module

40

! We define pi and i for the module

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

41

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

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

42

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

43

44

! All quadpack routines are private to the module

45

! All quadpack routines are private to the module

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

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

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

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

49

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

50

51

contains

52

contains

53

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

54

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

!

55

!

! Matrix inversion subroutines

56

! Matrix inversion subroutines

!

57

!

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

58

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

subroutine fvn_s_matinv(d,a,inva,status)

59

subroutine fvn_s_matinv(d,a,inva,status)

!

60

!

! Matrix inversion of a real matrix using BLAS and LAPACK

61

! Matrix inversion of a real matrix using BLAS and LAPACK

!

62

!

! d (in) : matrix rank

63

! d (in) : matrix rank

! a (in) : input matrix

64

! a (in) : input matrix

! inva (out) : inversed matrix

65

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

66

! status (ou) : =0 if something failed

!

67

!

68

implicit none

integer, intent(in) :: d

68

69

integer, intent(in) :: d

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

69

70

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

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

70

71

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

integer, intent(out) :: status

71

72

integer, intent(out) :: status

72

73

integer, allocatable :: ipiv(:)

73

74

integer, allocatable :: ipiv(:)

real, allocatable :: work(:)

74

75

real, allocatable :: work(:)

real twork(1)

75

76

real twork(1)

integer :: info

76

77

integer :: info

integer :: lwork

77

78

integer :: lwork

78

79

status=1

79

80

status=1

80

81

allocate(ipiv(d))

81

82

allocate(ipiv(d))

! copy a into inva using BLAS

82

83

! copy a into inva using BLAS

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

83

84

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

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

84

85

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

! LU factorization using LAPACK

85

86

! LU factorization using LAPACK

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

86

87

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

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

87

88

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

if (info /= 0) then

88

89

if (info /= 0) then

status=0

89

90

status=0

deallocate(ipiv)

90

91

deallocate(ipiv)

return

91

92

return

end if

92

93

end if

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

93

94

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

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

94

95

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

lwork=int(twork(1))

95

96

lwork=int(twork(1))

allocate(work(lwork))

96

97

allocate(work(lwork))

! Matrix inversion using LAPACK

97

98

! Matrix inversion using LAPACK

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

98

99

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

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

99

100

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

if (info /= 0) then

100

101

if (info /= 0) then

status=0

101

102

status=0

end if

102

103

end if

deallocate(work)

103

104

deallocate(work)

deallocate(ipiv)

104

105

deallocate(ipiv)

end subroutine

105

106

end subroutine

106

107

subroutine fvn_d_matinv(d,a,inva,status)

107

108

subroutine fvn_d_matinv(d,a,inva,status)

!

108

109

!

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

109

110

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

!

110

111

!

! d (in) : matrix rank

111

112

! d (in) : matrix rank

! a (in) : input matrix

112

113

! a (in) : input matrix

! inva (out) : inversed matrix

113

114

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

114

115

! status (ou) : =0 if something failed

!

115

116

!

117

implicit none

integer, intent(in) :: d

116

118

integer, intent(in) :: d

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

117

119

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

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

118

120

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

integer, intent(out) :: status

119

121

integer, intent(out) :: status

120

122

integer, allocatable :: ipiv(:)

121

123

integer, allocatable :: ipiv(:)

double precision, allocatable :: work(:)

122

124

double precision, allocatable :: work(:)

double precision :: twork(1)

123

125

double precision :: twork(1)

integer :: info

124

126

integer :: info

integer :: lwork

125

127

integer :: lwork

126

128

status=1

127

129

status=1

128

130

allocate(ipiv(d))

129

131

allocate(ipiv(d))

! copy a into inva using BLAS

130

132

! copy a into inva using BLAS

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

131

133

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

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

132

134

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

! LU factorization using LAPACK

133

135

! LU factorization using LAPACK

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

134

136

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

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

135

137

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

if (info /= 0) then

136

138

if (info /= 0) then

status=0

137

139

status=0

deallocate(ipiv)

138

140

deallocate(ipiv)

return

139

141

return

end if

140

142

end if

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

141

143

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

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

142

144

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

lwork=int(twork(1))

143

145

lwork=int(twork(1))

allocate(work(lwork))

144

146

allocate(work(lwork))

! Matrix inversion using LAPACK

145

147

! Matrix inversion using LAPACK

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

146

148

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

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

147

149

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

if (info /= 0) then

148

150

if (info /= 0) then

status=0

149

151

status=0

end if

150

152

end if

deallocate(work)

151

153

deallocate(work)

deallocate(ipiv)

152

154

deallocate(ipiv)

end subroutine

153

155

end subroutine

154

156

subroutine fvn_c_matinv(d,a,inva,status)

155

157

subroutine fvn_c_matinv(d,a,inva,status)

!

156

158

!

! Matrix inversion of a complex matrix using BLAS and LAPACK

157

159

! Matrix inversion of a complex matrix using BLAS and LAPACK

!

158

160

!

! d (in) : matrix rank

159

161

! d (in) : matrix rank

! a (in) : input matrix

160

162

! a (in) : input matrix

! inva (out) : inversed matrix

161

163

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

162

164

! status (ou) : =0 if something failed

!

163

165

!

166

implicit none

integer, intent(in) :: d

164

167

integer, intent(in) :: d

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

165

168

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

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

166

169

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

integer, intent(out) :: status

167

170

integer, intent(out) :: status

168

171

integer, allocatable :: ipiv(:)

169

172

integer, allocatable :: ipiv(:)

complex, allocatable :: work(:)

170

173

complex, allocatable :: work(:)

complex :: twork(1)

171

174

complex :: twork(1)

integer :: info

172

175

integer :: info

integer :: lwork

173

176

integer :: lwork

174

177

status=1

175

178

status=1

176

179

allocate(ipiv(d))

177

180

allocate(ipiv(d))

! copy a into inva using BLAS

178

181

! copy a into inva using BLAS

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

179

182

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

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

180

183

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

181

184

! LU factorization using LAPACK

182

185

! LU factorization using LAPACK

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

183

186

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

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

184

187

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

if (info /= 0) then

185

188

if (info /= 0) then

status=0

186

189

status=0

deallocate(ipiv)

187

190

deallocate(ipiv)

return

188

191

return

end if

189

192

end if

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

190

193

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

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

191

194

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

lwork=int(twork(1))

192

195

lwork=int(twork(1))

allocate(work(lwork))

193

196

allocate(work(lwork))

! Matrix inversion using LAPACK

194

197

! Matrix inversion using LAPACK

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

195

198

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

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

196

199

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

if (info /= 0) then

197

200

if (info /= 0) then

status=0

198

201

status=0

end if

199

202

end if

deallocate(work)

200

203

deallocate(work)

deallocate(ipiv)

201

204

deallocate(ipiv)

end subroutine

202

205

end subroutine

203

206

subroutine fvn_z_matinv(d,a,inva,status)

204

207

subroutine fvn_z_matinv(d,a,inva,status)

!

205

208

!

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

206

209

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

!

207

210

!

! d (in) : matrix rank

208

211

! d (in) : matrix rank

! a (in) : input matrix

209

212

! a (in) : input matrix

! inva (out) : inversed matrix

210

213

! inva (out) : inversed matrix

! status (ou) : =0 if something failed

211

214

! status (ou) : =0 if something failed

!

212

215

!

216

implicit none

integer, intent(in) :: d

213

217

integer, intent(in) :: d

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

214

218

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

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

215

219

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

integer, intent(out) :: status

216

220

integer, intent(out) :: status

217

221

integer, allocatable :: ipiv(:)

218

222

integer, allocatable :: ipiv(:)

double complex, allocatable :: work(:)

219

223

double complex, allocatable :: work(:)

double complex :: twork(1)

220

224

double complex :: twork(1)

integer :: info

221

225

integer :: info

integer :: lwork

222

226

integer :: lwork

223

227

status=1

224

228

status=1

225

229

allocate(ipiv(d))

226

230

allocate(ipiv(d))

! copy a into inva using BLAS

227

231

! copy a into inva using BLAS

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

228

232

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

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

229

233

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

230

234

! LU factorization using LAPACK

231

235

! LU factorization using LAPACK

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

232

236

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

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

233

237

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

if (info /= 0) then

234

238

if (info /= 0) then

status=0

235

239

status=0

deallocate(ipiv)

236

240

deallocate(ipiv)

return

237

241

return

end if

238

242

end if

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

239

243

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

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

240

244

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

lwork=int(twork(1))

241

245

lwork=int(twork(1))

allocate(work(lwork))

242

246

allocate(work(lwork))

! Matrix inversion using LAPACK

243

247

! Matrix inversion using LAPACK

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

244

248

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

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

245

249

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

if (info /= 0) then

246

250

if (info /= 0) then

status=0

247

251

status=0

end if

248

252

end if

deallocate(work)

249

253

deallocate(work)

deallocate(ipiv)

250

254

deallocate(ipiv)

end subroutine

251

255

end subroutine

252

256

253

257

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

254

258

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

!

255

259

!

! Determinants

256

260

! Determinants

!

257

261

!

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

258

262

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

function fvn_s_det(d,a,status)

259

263

function fvn_s_det(d,a,status)

!

260

264

!

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

261

265

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

!

262

266

!

! d (in) : matrix rank

263

267

! d (in) : matrix rank

! a (in) : The Matrix

264

268

! a (in) : The Matrix

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

265

269

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

!

266

270

!

271

implicit none

integer, intent(in) :: d

267

272

integer, intent(in) :: d

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

268

273

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

integer, intent(out) :: status

269

274

integer, intent(out) :: status

real :: fvn_s_det

270

275

real :: fvn_s_det

271

276

real, allocatable :: wc_a(:,:)

272

277

real, allocatable :: wc_a(:,:)

integer, allocatable :: ipiv(:)

273

278

integer, allocatable :: ipiv(:)

integer :: info,i

274

279

integer :: info,i

275

280

status=1

276

281

status=1

allocate(wc_a(d,d))

277

282

allocate(wc_a(d,d))

allocate(ipiv(d))

278

283

allocate(ipiv(d))

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

279

284

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

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

280

285

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

if (info/= 0) then

281

286

if (info/= 0) then

status=0

282

287

status=0

fvn_s_det=0.e0

283

288

fvn_s_det=0.e0

deallocate(ipiv)

284

289

deallocate(ipiv)

deallocate(wc_a)

285

290

deallocate(wc_a)

return

286

291

return

end if

287

292

end if

fvn_s_det=1.e0

288

293

fvn_s_det=1.e0

do i=1,d

289

294

do i=1,d

if (ipiv(i)==i) then

290

295

if (ipiv(i)==i) then

fvn_s_det=fvn_s_det*wc_a(i,i)

291

296

fvn_s_det=fvn_s_det*wc_a(i,i)

else

292

297

else

fvn_s_det=-fvn_s_det*wc_a(i,i)

293

298

fvn_s_det=-fvn_s_det*wc_a(i,i)

end if

294

299

end if

end do

295

300

end do

deallocate(ipiv)

296

301

deallocate(ipiv)

deallocate(wc_a)

297

302

deallocate(wc_a)

298

303

end function

299

304

end function

300

305

function fvn_d_det(d,a,status)

301

306

function fvn_d_det(d,a,status)

!

302

307

!

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

303

308

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

!

304

309

!

! d (in) : matrix rank

305

310

! d (in) : matrix rank

! a (in) : The Matrix

306

311

! a (in) : The Matrix

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

307

312

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

!

308

313

!

314

implicit none

integer, intent(in) :: d

309

315

integer, intent(in) :: d

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

310

316

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

integer, intent(out) :: status

311

317

integer, intent(out) :: status

double precision :: fvn_d_det

312

318

double precision :: fvn_d_det

313

319

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

314

320

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

integer, allocatable :: ipiv(:)

315

321

integer, allocatable :: ipiv(:)

integer :: info,i

316

322

integer :: info,i

317

323

status=1

318

324

status=1

allocate(wc_a(d,d))

319

325

allocate(wc_a(d,d))

allocate(ipiv(d))

320

326

allocate(ipiv(d))

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

321

327

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

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

322

328

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

if (info/= 0) then

323

329

if (info/= 0) then

status=0

324

330

status=0

fvn_d_det=0.d0

325

331

fvn_d_det=0.d0

deallocate(ipiv)

326

332

deallocate(ipiv)

deallocate(wc_a)

327

333

deallocate(wc_a)

return

328

334

return

end if

329

335

end if

fvn_d_det=1.d0

330

336

fvn_d_det=1.d0

do i=1,d

331

337

do i=1,d

if (ipiv(i)==i) then

332

338

if (ipiv(i)==i) then

fvn_d_det=fvn_d_det*wc_a(i,i)

333

339

fvn_d_det=fvn_d_det*wc_a(i,i)

else

334

340

else

fvn_d_det=-fvn_d_det*wc_a(i,i)

335

341

fvn_d_det=-fvn_d_det*wc_a(i,i)

end if

336

342

end if

end do

337

343

end do

deallocate(ipiv)

338

344

deallocate(ipiv)

deallocate(wc_a)

339

345

deallocate(wc_a)

340

346

end function

341

347

end function

342

348

function fvn_c_det(d,a,status) !

343

349

function fvn_c_det(d,a,status) !

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

344

350

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

!

345

351

!

! d (in) : matrix rank

346

352

! d (in) : matrix rank

! a (in) : The Matrix

347

353

! a (in) : The Matrix

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

348

354

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

!

349

355

!

356

implicit none

integer, intent(in) :: d

350

357

integer, intent(in) :: d

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

351

358

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

integer, intent(out) :: status

352

359

integer, intent(out) :: status

complex :: fvn_c_det

353

360

complex :: fvn_c_det

354

361

complex, allocatable :: wc_a(:,:)

355

362

complex, allocatable :: wc_a(:,:)

integer, allocatable :: ipiv(:)

356

363

integer, allocatable :: ipiv(:)

integer :: info,i

357

364

integer :: info,i

358

365

status=1

359

366

status=1

allocate(wc_a(d,d))

360

367

allocate(wc_a(d,d))

allocate(ipiv(d))

361

368

allocate(ipiv(d))

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

362

369

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

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

363

370

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

if (info/= 0) then

364

371

if (info/= 0) then

status=0

365

372

status=0

fvn_c_det=(0.e0,0.e0)

366

373

fvn_c_det=(0.e0,0.e0)

deallocate(ipiv)

367

374

deallocate(ipiv)

deallocate(wc_a)

368

375

deallocate(wc_a)

return

369

376

return

end if

370

377

end if

fvn_c_det=(1.e0,0.e0)

371

378

fvn_c_det=(1.e0,0.e0)

do i=1,d

372

379

do i=1,d

if (ipiv(i)==i) then

373

380

if (ipiv(i)==i) then

fvn_c_det=fvn_c_det*wc_a(i,i)

374

381

fvn_c_det=fvn_c_det*wc_a(i,i)

else

375

382

else

fvn_c_det=-fvn_c_det*wc_a(i,i)

376

383

fvn_c_det=-fvn_c_det*wc_a(i,i)

end if

377

384

end if

end do

378

385

end do

deallocate(ipiv)

379

386

deallocate(ipiv)

deallocate(wc_a)

380

387

deallocate(wc_a)

381

388

end function

382

389

end function

383

390

function fvn_z_det(d,a,status)

384

391

function fvn_z_det(d,a,status)

!

385

392

!

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

386

393

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

!

387

394

!

! d (in) : matrix rank

388

395

! d (in) : matrix rank

! a (in) : The Matrix

389

396

! a (in) : The Matrix

! det (out) : determinant

390

397

! det (out) : determinant

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

391

398

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

!

392

399

!

400

implicit none

integer, intent(in) :: d

393

401

integer, intent(in) :: d

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

394

402

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

integer, intent(out) :: status

395

403

integer, intent(out) :: status

double complex :: fvn_z_det

396

404

double complex :: fvn_z_det

397

405

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

398

406

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

integer, allocatable :: ipiv(:)

399

407

integer, allocatable :: ipiv(:)

integer :: info,i

400

408

integer :: info,i

401

409

status=1

402

410

status=1

allocate(wc_a(d,d))

403

411

allocate(wc_a(d,d))

allocate(ipiv(d))

404

412

allocate(ipiv(d))

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

405

413

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

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

406

414

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

if (info/= 0) then

407

415

if (info/= 0) then

status=0

408

416

status=0

fvn_z_det=(0.d0,0.d0)

409

417

fvn_z_det=(0.d0,0.d0)

deallocate(ipiv)

410

418

deallocate(ipiv)

deallocate(wc_a)

411

419

deallocate(wc_a)

return

412

420

return

end if

413

421

end if

fvn_z_det=(1.d0,0.d0)

414

422

fvn_z_det=(1.d0,0.d0)

do i=1,d

415

423

do i=1,d

if (ipiv(i)==i) then

416

424

if (ipiv(i)==i) then

fvn_z_det=fvn_z_det*wc_a(i,i)

417

425

fvn_z_det=fvn_z_det*wc_a(i,i)

else

418

426

else

fvn_z_det=-fvn_z_det*wc_a(i,i)

419

427

fvn_z_det=-fvn_z_det*wc_a(i,i)

end if

420

428

end if

end do

421

429

end do

deallocate(ipiv)

422

430

deallocate(ipiv)

deallocate(wc_a)

423

431

deallocate(wc_a)

424

432

end function

425

433

end function

426

434

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

427

435

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

!

428

436

!

! Condition test

429

437

! Condition test

!

430

438

!

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

431

439

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

! 1-norm

432

440

! 1-norm

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

433

441

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

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

434

442

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

!

435

443

!

subroutine fvn_s_matcon(d,a,rcond,status)

436

444

subroutine fvn_s_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

437

445

! Matrix condition (reciprocal of condition number)

!

438

446

!

! d (in) : matrix rank

439

447

! d (in) : matrix rank

! a (in) : The Matrix

440

448

! a (in) : The Matrix

! rcond (out) : guess what

441

449

! rcond (out) : guess what

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

442

450

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

!

443

451

!

452

implicit none

integer, intent(in) :: d

444

453

integer, intent(in) :: d

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

445

454

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

real, intent(out) :: rcond

446

455

real, intent(out) :: rcond

integer, intent(out) :: status

447

456

integer, intent(out) :: status

448

457

real, allocatable :: work(:)

449

458

real, allocatable :: work(:)

integer, allocatable :: iwork(:)

450

459

integer, allocatable :: iwork(:)

real :: anorm

451

460

real :: anorm

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

452

461

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

integer :: info

453

462

integer :: info

integer, allocatable :: ipiv(:)

454

463

integer, allocatable :: ipiv(:)

455

464

real, external :: slange

456

465

real, external :: slange

457

466

458

467

status=1

459

468

status=1

460

469

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

461

470

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

462

471

allocate(wc_a(d,d))

463

472

allocate(wc_a(d,d))

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

464

473

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

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

465

474

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

466

475

allocate(ipiv(d))

467

476

allocate(ipiv(d))

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

468

477

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

if (info /= 0) then

469

478

if (info /= 0) then

status=0

470

479

status=0

deallocate(ipiv)

471

480

deallocate(ipiv)

deallocate(wc_a)

472

481

deallocate(wc_a)

return

473

482

return

end if

474

483

end if

allocate(work(4*d))

475

484

allocate(work(4*d))

allocate(iwork(d))

476

485

allocate(iwork(d))

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

477

486

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

if (info /= 0) then

478

487

if (info /= 0) then

status=0

479

488

status=0

end if

480

489

end if

deallocate(iwork)

481

490

deallocate(iwork)

deallocate(work)

482

491

deallocate(work)

deallocate(ipiv)

483

492

deallocate(ipiv)

deallocate(wc_a)

484

493

deallocate(wc_a)

485

494

end subroutine

486

495

end subroutine

487

496

subroutine fvn_d_matcon(d,a,rcond,status)

488

497

subroutine fvn_d_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

489

498

! Matrix condition (reciprocal of condition number)

!

490

499

!

! d (in) : matrix rank

491

500

! d (in) : matrix rank

! a (in) : The Matrix

492

501

! a (in) : The Matrix

! rcond (out) : guess what

493

502

! rcond (out) : guess what

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

494

503

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

!

495

504

!

505

implicit none

integer, intent(in) :: d

496

506

integer, intent(in) :: d

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

497

507

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

double precision, intent(out) :: rcond

498

508

double precision, intent(out) :: rcond

integer, intent(out) :: status

499

509

integer, intent(out) :: status

500

510

double precision, allocatable :: work(:)

501

511

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

502

512

integer, allocatable :: iwork(:)

double precision :: anorm

503

513

double precision :: anorm

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

504

514

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

integer :: info

505

515

integer :: info

integer, allocatable :: ipiv(:)

506

516

integer, allocatable :: ipiv(:)

507

517

double precision, external :: dlange

508

518

double precision, external :: dlange

509

519

510

520

status=1

511

521

status=1

512

522

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

513

523

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

514

524

allocate(wc_a(d,d))

515

525

allocate(wc_a(d,d))

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

516

526

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

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

517

527

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

518

528

allocate(ipiv(d))

519

529

allocate(ipiv(d))

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

520

530

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

if (info /= 0) then

521

531

if (info /= 0) then

status=0

522

532

status=0

deallocate(ipiv)

523

533

deallocate(ipiv)

deallocate(wc_a)

524

534

deallocate(wc_a)

return

525

535

return

end if

526

536

end if

527

537

allocate(work(4*d))

528

538

allocate(work(4*d))

allocate(iwork(d))

529

539

allocate(iwork(d))

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

530

540

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

if (info /= 0) then

531

541

if (info /= 0) then

status=0

532

542

status=0

end if

533

543

end if

deallocate(iwork)

534

544

deallocate(iwork)

deallocate(work)

535

545

deallocate(work)

deallocate(ipiv)

536

546

deallocate(ipiv)

deallocate(wc_a)

537

547

deallocate(wc_a)

538

548

end subroutine

539

549

end subroutine

540

550

subroutine fvn_c_matcon(d,a,rcond,status)

541

551

subroutine fvn_c_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

542

552

! Matrix condition (reciprocal of condition number)

!

543

553

!

! d (in) : matrix rank

544

554

! d (in) : matrix rank

! a (in) : The Matrix

545

555

! a (in) : The Matrix

! rcond (out) : guess what

546

556

! rcond (out) : guess what

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

547

557

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

!

548

558

!

559

implicit none

integer, intent(in) :: d

549

560

integer, intent(in) :: d

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

550

561

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

real, intent(out) :: rcond

551

562

real, intent(out) :: rcond

integer, intent(out) :: status

552

563

integer, intent(out) :: status

553

564

real, allocatable :: rwork(:)

554

565

real, allocatable :: rwork(:)

complex, allocatable :: work(:)

555

566

complex, allocatable :: work(:)

integer, allocatable :: iwork(:)

556

567

integer, allocatable :: iwork(:)

real :: anorm

557

568

real :: anorm

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

558

569

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

integer :: info

559

570

integer :: info

integer, allocatable :: ipiv(:)

560

571

integer, allocatable :: ipiv(:)

561

572

real, external :: clange

562

573

real, external :: clange

563

574

564

575

status=1

565

576

status=1

566

577

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

567

578

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

568

579

allocate(wc_a(d,d))

569

580

allocate(wc_a(d,d))

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

570

581

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

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

571

582

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

572

583

allocate(ipiv(d))

573

584

allocate(ipiv(d))

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

574

585

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

if (info /= 0) then

575

586

if (info /= 0) then

status=0

576

587

status=0

deallocate(ipiv)

577

588

deallocate(ipiv)

deallocate(wc_a)

578

589

deallocate(wc_a)

return

579

590

return

end if

580

591

end if

allocate(work(2*d))

581

592

allocate(work(2*d))

allocate(rwork(2*d))

582

593

allocate(rwork(2*d))

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

583

594

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

if (info /= 0) then

584

595

if (info /= 0) then

status=0

585

596

status=0

end if

586

597

end if

deallocate(rwork)

587

598

deallocate(rwork)

deallocate(work)

588

599

deallocate(work)

deallocate(ipiv)

589

600

deallocate(ipiv)

deallocate(wc_a)

590

601

deallocate(wc_a)

end subroutine

591

602

end subroutine

592

603

subroutine fvn_z_matcon(d,a,rcond,status)

593

604

subroutine fvn_z_matcon(d,a,rcond,status)

! Matrix condition (reciprocal of condition number)

594

605

! Matrix condition (reciprocal of condition number)

!

595

606

!

! d (in) : matrix rank

596

607

! d (in) : matrix rank

! a (in) : The Matrix

597

608

! a (in) : The Matrix

! rcond (out) : guess what

598

609

! rcond (out) : guess what

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

599

610

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

!

600

611

!

612

implicit none

integer, intent(in) :: d

601

613

integer, intent(in) :: d

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

602

614

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

double precision, intent(out) :: rcond

603

615

double precision, intent(out) :: rcond

integer, intent(out) :: status

604

616

integer, intent(out) :: status

605

617

double complex, allocatable :: work(:)

606

618

double complex, allocatable :: work(:)

double precision, allocatable :: rwork(:)

607

619

double precision, allocatable :: rwork(:)

double precision :: anorm

608

620

double precision :: anorm

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

609

621

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

integer :: info

610

622

integer :: info

integer, allocatable :: ipiv(:)

611

623

integer, allocatable :: ipiv(:)

612

624

double precision, external :: zlange

613

625

double precision, external :: zlange

614

626

615

627

status=1

616

628

status=1

617

629

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

618

630

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

619

631

allocate(wc_a(d,d))

620

632

allocate(wc_a(d,d))

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

621

633

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

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

622

634

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

623

635

allocate(ipiv(d))

624

636

allocate(ipiv(d))

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

625

637

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

if (info /= 0) then

626

638

if (info /= 0) then

status=0

627

639

status=0

deallocate(ipiv)

628

640

deallocate(ipiv)

deallocate(wc_a)

629

641

deallocate(wc_a)

return

630

642

return

end if

631

643

end if

632

644

allocate(work(2*d))

633

645

allocate(work(2*d))

allocate(rwork(2*d))

634

646

allocate(rwork(2*d))

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

635

647

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

if (info /= 0) then

636

648

if (info /= 0) then

status=0

637

649

status=0

end if

638

650

end if

deallocate(rwork)

639

651

deallocate(rwork)

deallocate(work)

640

652

deallocate(work)

deallocate(ipiv)

641

653

deallocate(ipiv)

deallocate(wc_a)

642

654

deallocate(wc_a)

end subroutine

643

655

end subroutine

644

656

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

645

657

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

!

646

658

!

! Valeurs propres/ Vecteurs propre

647

659

! Valeurs propres/ Vecteurs propre

!

648

660

!

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

649

661

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

650

662

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

651

663

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

!

652

664

!

! integer d (in) : matrice rank

653

665

! integer d (in) : matrice rank

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

654

666

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

! complex evala(d) (out) : eigenvalues

655

667

! complex evala(d) (out) : eigenvalues

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

656

668

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

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

657

669

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

!

658

670

!

! interfacing Lapack routine SGEEV

659

671

! interfacing Lapack routine SGEEV

660

672

implicit none

integer, intent(in) :: d

661

673

integer, intent(in) :: d

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

662

674

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

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

663

675

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

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

664

676

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

integer, intent(out) :: status

665

677

integer, intent(out) :: status

666

678

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

667

679

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

integer :: info

668

680

integer :: info

integer :: lwork

669

681

integer :: lwork

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

670

682

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

real :: vl ! unused but necessary for the call

671

683

real :: vl ! unused but necessary for the call

real, allocatable :: vr(:,:)

672

684

real, allocatable :: vr(:,:)

real, allocatable :: work(:)

673

685

real, allocatable :: work(:)

real :: twork(1)

674

686

real :: twork(1)

integer i

675

687

integer i

integer j

676

688

integer j

677

689

! making a working copy of a

678

690

! making a working copy of a

allocate(wc_a(d,d))

679

691

allocate(wc_a(d,d))

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

680

692

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

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

681

693

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

682

694

allocate(wr(d))

683

695

allocate(wr(d))

allocate(wi(d))

684

696

allocate(wi(d))

allocate(vr(d,d))

685

697

allocate(vr(d,d))

! query optimal work size

686

698

! query optimal work size

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

687

699

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

lwork=int(twork(1))

688

700

lwork=int(twork(1))

allocate(work(lwork))

689

701

allocate(work(lwork))

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

690

702

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

691

703

if (info /= 0) then

692

704

if (info /= 0) then

status=0

693

705

status=0

deallocate(work)

694

706

deallocate(work)

deallocate(vr)

695

707

deallocate(vr)

deallocate(wi)

696

708

deallocate(wi)

deallocate(wr)

697

709

deallocate(wr)

deallocate(wc_a)

698

710

deallocate(wc_a)

return

699

711

return

end if

700

712

end if

701

713

! now fill in the results

702

714

! now fill in the results

i=1

703

715

i=1

do while(i<=d)

704

716

do while(i<=d)

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

705

717

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

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

706

718

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

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

707

719

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

else ! eigenvalue is complex

708

720

else ! eigenvalue is complex

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

709

721

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

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

710

722

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

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

711

723

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

i=i+1

712

724

i=i+1

end if

713

725

end if

i=i+1

714

726

i=i+1

enddo

715

727

enddo

deallocate(work)

716

728

deallocate(work)

deallocate(vr)

717

729

deallocate(vr)

deallocate(wi)

718

730

deallocate(wi)

deallocate(wr)

719

731

deallocate(wr)

deallocate(wc_a)

720

732

deallocate(wc_a)

721

733

end subroutine

722

734

end subroutine

723

735

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

724

736

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

!

725

737

!

! integer d (in) : matrice rank

726

738

! integer d (in) : matrice rank

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

727

739

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

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

728

740

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

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

729

741

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

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

730

742

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

!

731

743

!

! interfacing Lapack routine DGEEV

732

744

! interfacing Lapack routine DGEEV

745

implicit none

integer, intent(in) :: d

733

746

integer, intent(in) :: d

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

734

747

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

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

735

748

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

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

736

749

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

integer, intent(out) :: status

737

750

integer, intent(out) :: status

738

751

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

739

752

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

integer :: info

740

753

integer :: info

integer :: lwork

741

754

integer :: lwork

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

742

755

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

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

743

756

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

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

744

757

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

double precision, allocatable :: work(:)

745

758

double precision, allocatable :: work(:)

double precision :: twork(1)

746

759

double precision :: twork(1)

integer i

747

760

integer i

integer j

748

761

integer j

749

762

! making a working copy of a

750

763

! making a working copy of a

allocate(wc_a(d,d))

751

764

allocate(wc_a(d,d))

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

752

765

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

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

753

766

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

754

767

allocate(wr(d))

755

768

allocate(wr(d))

allocate(wi(d))

756

769

allocate(wi(d))

allocate(vr(d,d))

757

770

allocate(vr(d,d))

! query optimal work size

758

771

! query optimal work size

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

759

772

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

lwork=int(twork(1))

760

773

lwork=int(twork(1))

allocate(work(lwork))

761

774

allocate(work(lwork))

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

762

775

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

763

776

if (info /= 0) then

764

777

if (info /= 0) then

status=0

765

778

status=0

deallocate(work)

766

779

deallocate(work)

deallocate(vr)

767

780

deallocate(vr)

deallocate(wi)

768

781

deallocate(wi)

deallocate(wr)

769

782

deallocate(wr)

deallocate(wc_a)

770

783

deallocate(wc_a)

return

771

784

return

end if

772

785

end if

773

786

! now fill in the results

774

787

! now fill in the results

i=1

775

788

i=1

do while(i<=d)

776

789

do while(i<=d)

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

777

790

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

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

778

791

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

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

779

792

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

else ! eigenvalue is complex

780

793

else ! eigenvalue is complex

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

781

794

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

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

782

795

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

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

783

796

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

i=i+1

784

797

i=i+1

end if

785

798

end if

i=i+1

786

799

i=i+1

enddo

787

800

enddo

788

801

deallocate(work)

789

802

deallocate(work)

deallocate(vr)

790

803

deallocate(vr)

deallocate(wi)

791

804

deallocate(wi)

deallocate(wr)

792

805

deallocate(wr)

deallocate(wc_a)

793

806

deallocate(wc_a)

794

807

end subroutine

795

808

end subroutine

796

809

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

797

810

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

!

798

811

!

! integer d (in) : matrice rank

799

812

! integer d (in) : matrice rank

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

800

813

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

! complex evala(d) (out) : eigenvalues

801

814

! complex evala(d) (out) : eigenvalues

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

802

815

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

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

803

816

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

!

804

817

!

! interfacing Lapack routine CGEEV

805

818

! interfacing Lapack routine CGEEV

806

819

implicit none

integer, intent(in) :: d

807

820

integer, intent(in) :: d

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

808

821

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

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

809

822

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

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

810

823

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

integer, intent(out) :: status

811

824

integer, intent(out) :: status

812

825

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

813

826

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

integer :: info

814

827

integer :: info

integer :: lwork

815

828

integer :: lwork

complex, allocatable :: work(:)

816

829

complex, allocatable :: work(:)

complex :: twork(1)

817

830

complex :: twork(1)

real, allocatable :: rwork(:)

818

831

real, allocatable :: rwork(:)

complex :: vl ! unused but necessary for the call

819

832

complex :: vl ! unused but necessary for the call

820

833

status=1

821

834

status=1

822

835

! making a working copy of a

823

836

! making a working copy of a

allocate(wc_a(d,d))

824

837

allocate(wc_a(d,d))

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

825

838

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

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

826

839

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

827

840

828

841

! query optimal work size

829

842

! query optimal work size

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

830

843

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

lwork=int(twork(1))

831

844

lwork=int(twork(1))

allocate(work(lwork))

832

845

allocate(work(lwork))

allocate(rwork(2*d))

833

846

allocate(rwork(2*d))

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

834

847

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

835

848

if (info /= 0) then

836

849

if (info /= 0) then

status=0

837

850

status=0

end if

838

851

end if

deallocate(rwork)

839

852

deallocate(rwork)

deallocate(work)

840

853

deallocate(work)

deallocate(wc_a)

841

854

deallocate(wc_a)

842

855

end subroutine

843

856

end subroutine

844

857

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

845

858

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

!

846

859

!

! integer d (in) : matrice rank

847

860

! integer d (in) : matrice rank

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

848

861

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

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

849

862

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

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

850

863

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

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

851

864

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

!

852

865

!

! interfacing Lapack routine ZGEEV

853

866

! interfacing Lapack routine ZGEEV

854

867

implicit none

integer, intent(in) :: d

855

868

integer, intent(in) :: d

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

856

869

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

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

857

870

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

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

858

871

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

integer, intent(out) :: status

859

872

integer, intent(out) :: status

860

873

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

861

874

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

integer :: info

862

875

integer :: info

integer :: lwork

863

876

integer :: lwork

double complex, allocatable :: work(:)

864

877

double complex, allocatable :: work(:)

double complex :: twork(1)

865

878

double complex :: twork(1)

double precision, allocatable :: rwork(:)

866

879

double precision, allocatable :: rwork(:)

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

867

880

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

868

881

status=1

869

882

status=1

870

883

! making a working copy of a

871

884

! making a working copy of a

allocate(wc_a(d,d))

872

885

allocate(wc_a(d,d))

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

873

886

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

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

874

887

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

875

888

! query optimal work size

876

889

! query optimal work size

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

877

890

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

lwork=int(twork(1))

878

891

lwork=int(twork(1))

allocate(work(lwork))

879

892

allocate(work(lwork))

allocate(rwork(2*d))

880

893

allocate(rwork(2*d))

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

881

894

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

882

895

if (info /= 0) then

883

896

if (info /= 0) then

status=0

884

897

status=0

end if

885

898

end if

deallocate(rwork)

886

899

deallocate(rwork)

deallocate(work)

887

900

deallocate(work)

deallocate(wc_a)

888

901

deallocate(wc_a)

889

902

end subroutine

890

903

end subroutine

891

904

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

892

905

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

!

893

906

!

! Akima spline interpolation and spline evaluation

894

907

! Akima spline interpolation and spline evaluation

!

895

908

!

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

896

909

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

897

910

! Single precision

898

911

! Single precision

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

899

912

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

implicit none

900

913

implicit none

integer, intent(in) :: n

901

914

integer, intent(in) :: n

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

902

915

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

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

903

916

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

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

904

917

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

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

905

918

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

906

919

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

907

920

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

real :: wi_1,wi

908

921

real :: wi_1,wi

integer :: i

909

922

integer :: i

real :: dx,a,b

910

923

real :: dx,a,b

911

924

! br is just a copy of x

912

925

! br is just a copy of x

br(:)=x(:)

913

926

br(:)=x(:)

914

927

allocate(var(n))

915

928

allocate(var(n))

allocate(z(n))

916

929

allocate(z(n))

! evaluate the variations

917

930

! evaluate the variations

do i=1, n-1

918

931

do i=1, n-1

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

919

932

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

end do

920

933

end do

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

921

934

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

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

922

935

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

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

923

936

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

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

924

937

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

925

938

do i = 1, n

926

939

do i = 1, n

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

927

940

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

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

928

941

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

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

929

942

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

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

930

943

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

else

931

944

else

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

932

945

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

end if

933

946

end if

end do

934

947

end do

935

948

do i=1, n-1

936

949

do i=1, n-1

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

937

950

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

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

938

951

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

939

952

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

co(1,i)=y(i)

940

953

co(1,i)=y(i)

co(2,i)=z(i)

941

954

co(2,i)=z(i)

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

942

955

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

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

943

956

!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

944

957

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 !

945

958

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

946

959

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

947

960

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

end do

948

961

end do

co(1,n)=y(n)

949

962

co(1,n)=y(n)

co(2,n)=z(n)

950

963

co(2,n)=z(n)

co(3,n)=0.e0

951

964

co(3,n)=0.e0

co(4,n)=0.e0

952

965

co(4,n)=0.e0

953

966

deallocate(z)

954

967

deallocate(z)

deallocate(var)

955

968

deallocate(var)

956

969

end subroutine

957

970

end subroutine

958

971

! Double precision

959

972

! Double precision

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

960

973

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

961

974

implicit none

962

975

implicit none

integer, intent(in) :: n

963

976

integer, intent(in) :: n

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

964

977

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

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

965

978

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

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

966

979

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

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

967

980

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

968

981

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

969

982

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

double precision :: wi_1,wi

970

983

double precision :: wi_1,wi

integer :: i

971

984

integer :: i

double precision :: dx,a,b

972

985

double precision :: dx,a,b

973

986

! br is just a copy of x

974

987

! br is just a copy of x

br(:)=x(:)

975

988

br(:)=x(:)

976

989

allocate(var(n))

977

990

allocate(var(n))

allocate(z(n))

978

991

allocate(z(n))

! evaluate the variations

979

992

! evaluate the variations

do i=1, n-1

980

993

do i=1, n-1

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

981

994

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

end do

982

995

end do

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

983

996

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

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

984

997

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

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

985

998

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

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

986

999

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

987

1000

do i = 1, n

988

1001

do i = 1, n

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

989

1002

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

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

990

1003

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

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

991

1004

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

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

992

1005

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

else

993

1006

else

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

994

1007

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

end if

995

1008

end if

end do

996

1009

end do

997

1010

do i=1, n-1

998

1011

do i=1, n-1

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

999

1012

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

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

1000

1013

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

1001

1014

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

co(1,i)=y(i)

1002

1015

co(1,i)=y(i)

co(2,i)=z(i)

1003

1016

co(2,i)=z(i)

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

1004

1017

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

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

1005

1018

!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

1006

1019

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 !

1007

1020

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

1008

1021

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

1009

1022

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

end do

1010

1023

end do

co(1,n)=y(n)

1011

1024

co(1,n)=y(n)

co(2,n)=z(n)

1012

1025

co(2,n)=z(n)

co(3,n)=0.d0

1013

1026

co(3,n)=0.d0

co(4,n)=0.d0

1014

1027

co(4,n)=0.d0

1015

1028

deallocate(z)

1016

1029

deallocate(z)

deallocate(var)

1017

1030

deallocate(var)

1018

1031

end subroutine

1019

1032

end subroutine

1020

1033

!

1021

1034

!

! Single precision spline evaluation

1022

1035

! Single precision spline evaluation

!

1023

1036

!

function fvn_s_spline_eval(x,n,br,co)

1024

1037

function fvn_s_spline_eval(x,n,br,co)

implicit none

1025

1038

implicit none

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

1026

1039

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

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

1027

1040

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

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

1028

1041

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

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

1029

1042

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

real :: fvn_s_spline_eval

1030

1043

real :: fvn_s_spline_eval

1031

1044

integer :: i

1032

1045

integer :: i

real :: dx

1033

1046

real :: dx

1034

1047

if (x<=br(1)) then

1035

1048

if (x<=br(1)) then

i=1

1036

1049

i=1

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

1037

1050

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

i=n

1038

1051

i=n

else

1039

1052

else

i=1

1040

1053

i=1

do while(x>=br(i))

1041

1054

do while(x>=br(i))

i=i+1

1042

1055

i=i+1

end do

1043

1056

end do

i=i-1

1044

1057

i=i-1

end if

1045

1058

end if

dx=x-br(i)

1046

1059

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

1047

1060

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

1048

1061

end function

1049

1062

end function

1050

1063

! Double precision spline evaluation

1051

1064

! Double precision spline evaluation

function fvn_d_spline_eval(x,n,br,co)

1052

1065

function fvn_d_spline_eval(x,n,br,co)

implicit none

1053

1066

implicit none

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

1054

1067

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

1055

1068

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

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

1056

1069

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

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

1057

1070

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

double precision :: fvn_d_spline_eval

1058

1071

double precision :: fvn_d_spline_eval

1059

1072

integer :: i

1060

1073

integer :: i

double precision :: dx

1061

1074

double precision :: dx

1062

1075

1063

1076

if (x<=br(1)) then

1064

1077

if (x<=br(1)) then

i=1

1065

1078

i=1

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

1066

1079

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

i=n

1067

1080

i=n

else

1068

1081

else

i=1

1069

1082

i=1

do while(x>=br(i))

1070

1083

do while(x>=br(i))

i=i+1

1071

1084

i=i+1

end do

1072

1085

end do

i=i-1

1073

1086

i=i-1

end if

1074

1087

end if

1075

1088

dx=x-br(i)

1076

1089

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

1077

1090

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

1078

1091

end function

1079

1092

end function

1080

1093

1081

1094

1095

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

1096

!

1097

! Least square problem

!

1082

1098

!

1099

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

1100

!

1101

!

1102

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

1103

!

1104

! Least square polynomial fitting

1105

!

1106

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

1107

! for a set of coordinates.

1108

!

1109

! The degree must be lower than the number of points

1110

!

1111

! np (in) : number of points

1112

! x(np) (in) : x data

1113

! y(np) (in) : y data

1114

! deg (in) : polynomial's degree

1115

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

1116

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

1117

!

1118

implicit none

1119

1120

integer, intent(in) :: np,deg

1121

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

1122

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

1123

integer, intent(out) :: status

1124

1125

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

1126

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

1127

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

1128

integer :: lwork,info,rank

1129

1130

integer :: i,j

1131

1132

status=1

1133

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

1134

1135

! Design matrix valorisation

1136

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

1137

1138

! second member valorisation

1139

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

1140

1141

! query workspace size

1142

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

1143

lwork=twork(1)

1144

allocate(work(int(lwork)))

1145

! real call

1146

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

1147

1148

if (info /= 0) then

1149

status=0

1150

end if

1151

1152

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

1153

1154

deallocate(work)

1155

deallocate(mat,bmat,singval)

1156

end subroutine

1157

1158

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

1159

!

1160

! Least square polynomial fitting

1161

!

1162

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

1163

! for a set of coordinates.

1164

!

1165

! The degree must be lower than the number of points

1166

!

1167

! np (in) : number of points

1168

! x(np) (in) : x data

1169

! y(np) (in) : y data

1170

! deg (in) : polynomial's degree

1171

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

1172

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

1173

!

1174

implicit none

1175

1176

integer, intent(in) :: np,deg

1177

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

1178

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

1179

integer, intent(out) :: status

1180

1181

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

1182

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

1183

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

1184

integer :: lwork,info,rank

1185

1186

integer :: i,j

1187

1188

status=1

1189

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

1190

1191

! Design matrix valorisation

1192

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

1193

1194

! second member valorisation

1195

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

1196

1197

! query workspace size

1198

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

1199

lwork=twork(1)

1200

allocate(work(int(lwork)))

1201

! real call

1202

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

1203

1204

if (info /= 0) then

1205

status=0

1206

end if

1207

1208

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

1209

1210

deallocate(work)

1211

deallocate(mat,bmat,singval)

1212

end subroutine

1213

1214

1215

!

! Muller

1083

1216

! Muller

!

1084

1217

!

1085

1218

!

1086

1219

!

! William Daniau 2007

1087

1220

! William Daniau 2007

!

1088

1221

!

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

1089

1222

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

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

1090

1223

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

!

1091

1224

!

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

1092

1225

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

!

1093

1226

!

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

1094

1227

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

!

1095

1228

!

! purpose - zeros of an analytic complex function

1096

1229

! purpose - zeros of an analytic complex function

! using the muller method with deflation

1097

1230

! using the muller method with deflation

!

1098

1231

!

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

1099

1232

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

! infer,ier)

1100

1233

! infer,ier)

!

1101

1234

!

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

1102

1235

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

! by the user specifying the equation whose

1103

1236

! by the user specifying the equation whose

! roots are to be found. f must appear in

1104

1237

! roots are to be found. f must appear in

! an external statement in the calling pro-

1105

1238

! an external statement in the calling pro-

! gram.

1106

1239

! gram.

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

1107

1240

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

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

1108

1241

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

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

1109

1242

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

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

1110

1243

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

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

1111

1244

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

! as a root. (input)

1112

1245

! as a root. (input)

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

1113

1246

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

! if two successive approximations to a given

1114

1247

! if two successive approximations to a given

! root agree within eps1. (input)

1115

1248

! root agree within eps1. (input)

! note. if either or both of the stopping

1116

1249

! note. if either or both of the stopping

! criteria are fulfilled, the root is

1117

1250

! criteria are fulfilled, the root is

! accepted.

1118

1251

! accepted.

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

1119

1252

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

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

1120

1253

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

! nguess - the number of initial guesses provided. these

1121

1254

! nguess - the number of initial guesses provided. these

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

1122

1255

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

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

1123

1256

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

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

1124

1257

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

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

1125

1258

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

! muller (input)

1126

1259

! muller (input)

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

1127

1260

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

! x(kn) on input must contain any known

1128

1261

! x(kn) on input must contain any known

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

1129

1262

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

! on user option, contain initial guesses for

1130

1263

! on user option, contain initial guesses for

! the n new roots which are to be computed.

1131

1264

! the n new roots which are to be computed.

! if the user does not provide an initial

1132

1265

! if the user does not provide an initial

! guess, zero is used.

1133

1266

! guess, zero is used.

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

1134

1267

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

! approximate roots found by muller.

1135

1268

! approximate roots found by muller.

! itmax - the maximum allowable number of iterations

1136

1269

! itmax - the maximum allowable number of iterations

! per root (input)

1137

1270

! per root (input)

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

1138

1271

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

! output infer(j) contains the number of

1139

1272

! output infer(j) contains the number of

! iterations used in finding the j-th root

1140

1273

! iterations used in finding the j-th root

! when convergence was achieved. if

1141

1274

! when convergence was achieved. if

! convergence was not obtained in itmax

1142

1275

! convergence was not obtained in itmax

! iterations, infer(j) will be greater than

1143

1276

! iterations, infer(j) will be greater than

! itmax (output).

1144

1277

! itmax (output).

! ier - error parameter (output)

1145

1278

! ier - error parameter (output)

! warning error

1146

1279

! warning error

! ier = 33 indicates failure to converge with-

1147

1280

! ier = 33 indicates failure to converge with-

! in itmax iterations for at least one of

1148

1281

! in itmax iterations for at least one of

! the (n) new roots.

1149

1282

! the (n) new roots.

!

1150

1283

!

1151

1284

!

! remarks muller always returns the last approximation for root j

1152

1285

! remarks muller always returns the last approximation for root j

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

1153

1286

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

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

1154

1287

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

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

1155

1288

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

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

1156

1289

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

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

1157

1290

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

! not obtain convergence in the allowed number of iter-

1158

1291

! not obtain convergence in the allowed number of iter-

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

1159

1292

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

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

1160

1293

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

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

1161

1294

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

! ted function

1162

1295

! ted function

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

1163

1296

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

!

1164

1297

!

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

1165

1298

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

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

1166

1299

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

! the convergence criterion.

1167

1300

! the convergence criterion.

!

1168

1301

!

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

1169

1302

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

!

1170

1303

!

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

1171

1304

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

implicit none

1172

1305

implicit none

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

1173

1306

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

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

1174

1307

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

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

1175

1308

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

zero,p1,one,four,p5

1176

1309

zero,p1,one,four,p5

1177

1310

double complex, external :: f

1178

1311

double complex, external :: f

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

1179

1312

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

knpng,jk,ick,nn,lm1,errcode

1180

1313

knpng,jk,ick,nn,lm1,errcode

double complex :: x(kn+n)

1181

1314

double complex :: x(kn+n)

integer :: infer(kn+n)

1182

1315

integer :: infer(kn+n)

1183

1316

1184

1317

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

1185

1318

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

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

1186

1319

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

p5/(0.5,0.0)/, &

1187

1320

p5/(0.5,0.0)/, &

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

1188

1321

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

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

1189

1322

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

1190

1323

ier = 0

1191

1324

ier = 0

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

1192

1325

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

return

1193

1326

return

end if

1194

1327

end if

!eps1 = rten **(-nsig)

1195

1328

!eps1 = rten **(-nsig)

eps1 = min(eps1,rp01)

1196

1329

eps1 = min(eps1,rp01)

1197

1330

knp1 = kn+1

1198

1331

knp1 = kn+1

knpn = kn+n

1199

1332

knpn = kn+n

knpng = kn+nguess

1200

1333

knpng = kn+nguess

do i=1,knpn

1201

1334

do i=1,knpn

infer(i) = 0

1202

1335

infer(i) = 0

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

1203

1336

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

end do

1204

1337

end do

l= knp1

1205

1338

l= knp1

1206

1339

ic=0

1207

1340

ic=0

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

1208

1341

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

jk = 0

1209

1342

jk = 0

ick = 0

1210

1343

ick = 0

xl = x(l)

1211

1344

xl = x(l)

icloop: do

1212

1345

icloop: do

ic = 0

1213

1346

ic = 0

h = ax

1214

1347

h = ax

h = p1*h

1215

1348

h = p1*h

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

1216

1349

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

! first three points are

1217

1350

! first three points are

! xl+h, xl-h, xl

1218

1351

! xl+h, xl-h, xl

rt = xl+h

1219

1352

rt = xl+h

call deflated_work(errcode)

1220

1353

call deflated_work(errcode)

if (errcode == 1) then

1221

1354

if (errcode == 1) then

exit icloop

1222

1355

exit icloop

end if

1223

1356

end if

1224

1357

z0 = fprt

1225

1358

z0 = fprt

y0 = frt

1226

1359

y0 = frt

x0 = rt

1227

1360

x0 = rt

rt = xl-h

1228

1361

rt = xl-h

call deflated_work(errcode)

1229

1362

call deflated_work(errcode)

if (errcode == 1) then

1230

1363

if (errcode == 1) then

exit icloop

1231

1364

exit icloop

end if

1232

1365

end if

1233

1366

z1 = fprt

1234

1367

z1 = fprt

y1 = frt

1235

1368

y1 = frt

h = xl-rt

1236

1369

h = xl-rt

d = h/(rt-x0)

1237

1370

d = h/(rt-x0)

rt = xl

1238

1371

rt = xl

1239

1372

call deflated_work(errcode)

1240

1373

call deflated_work(errcode)

if (errcode == 1) then

1241

1374

if (errcode == 1) then

exit icloop

1242

1375

exit icloop

end if

1243

1376

end if

1244

1377

1245

1378

z2 = fprt

1246

1379

z2 = fprt

y2 = frt

1247

1380

y2 = frt

! begin main algorithm

1248

1381

! begin main algorithm

iloop: do

1249

1382

iloop: do

dd = one + d

1250

1383

dd = one + d

t1 = z0*d*d

1251

1384

t1 = z0*d*d

t2 = z1*dd*dd

1252

1385

t2 = z1*dd*dd

xx = z2*dd

1253

1386

xx = z2*dd

t3 = z2*d

1254

1387

t3 = z2*d

bi = t1-t2+xx+t3

1255

1388

bi = t1-t2+xx+t3

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

1256

1389

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

! use denominator of maximum amplitude

1257

1390

! use denominator of maximum amplitude

t1 = cdsqrt(den)

1258

1391

t1 = cdsqrt(den)

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

1259

1392

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

t2 = bi + t1

1260

1393

t2 = bi + t1

tpq = cdabs(t2)+qz

1261

1394

tpq = cdabs(t2)+qz

if (tpq .eq. qz) t2 = zero

1262

1395

if (tpq .eq. qz) t2 = zero

t3 = bi - t1

1263

1396

t3 = bi - t1

tpq = cdabs(t3) + qz

1264

1397

tpq = cdabs(t3) + qz

if (tpq .eq. qz) t3 = zero

1265

1398

if (tpq .eq. qz) t3 = zero

den = t2

1266

1399

den = t2

qz = cdabs(t3)-cdabs(t2)

1267

1400

qz = cdabs(t3)-cdabs(t2)

if (qz .gt. rzero) den = t3

1268

1401

if (qz .gt. rzero) den = t3

! test for zero denominator

1269

1402

! test for zero denominator

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

1270

1403

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

call trans_rt()

1271

1404

call trans_rt()

call deflated_work(errcode)

1272

1405

call deflated_work(errcode)

if (errcode == 1) then

1273

1406

if (errcode == 1) then

exit icloop

1274

1407

exit icloop

end if

1275

1408

end if

z2 = fprt

1276

1409

z2 = fprt

y2 = frt

1277

1410

y2 = frt

cycle iloop

1278

1411

cycle iloop

end if

1279

1412

end if

1280

1413

1281

1414

d = -xx/den

1282

1415

d = -xx/den

d = d+d

1283

1416

d = d+d

h = d*h

1284

1417

h = d*h

rt = rt + h

1285

1418

rt = rt + h

! check convergence of the first kind

1286

1419

! check convergence of the first kind

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

1287

1420

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

if (ic .ne. 0) then

1288

1421

if (ic .ne. 0) then

exit icloop

1289

1422

exit icloop

end if

1290

1423

end if

ic = 1

1291

1424

ic = 1

z0 = y1

1292

1425

z0 = y1

z1 = y2

1293

1426

z1 = y2

z2 = f(rt)

1294

1427

z2 = f(rt)

xl = rt

1295

1428

xl = rt

ick = ick+1

1296

1429

ick = ick+1

if (ick .le. ickmax) then

1297

1430

if (ick .le. ickmax) then

cycle iloop

1298

1431

cycle iloop

end if

1299

1432

end if

! warning error, itmax = maximum

1300

1433

! warning error, itmax = maximum

jk = itmax + jk

1301

1434

jk = itmax + jk

ier = 33

1302

1435

ier = 33

end if

1303

1436

end if

if (ic .ne. 0) then

1304

1437

if (ic .ne. 0) then

cycle icloop

1305

1438

cycle icloop

end if

1306

1439

end if

call deflated_work(errcode)

1307

1440

call deflated_work(errcode)

if (errcode == 1) then

1308

1441

if (errcode == 1) then

exit icloop

1309

1442

exit icloop

end if

1310

1443

end if

1311

1444

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

1312

1445

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

! take remedial action to induce

1313

1446

! take remedial action to induce

! convergence

1314

1447

! convergence

d = d*p5

1315

1448

d = d*p5

h = h*p5

1316

1449

h = h*p5

rt = rt-h

1317

1450

rt = rt-h

call deflated_work(errcode)

1318

1451

call deflated_work(errcode)

if (errcode == 1) then

1319

1452

if (errcode == 1) then

exit icloop

1320

1453

exit icloop

end if

1321

1454

end if

end do

1322

1455

end do

z0 = z1

1323

1456

z0 = z1

z1 = z2

1324

1457

z1 = z2

z2 = fprt

1325

1458

z2 = fprt

y0 = y1

1326

1459

y0 = y1

y1 = y2

1327

1460

y1 = y2

y2 = frt

1328

1461

y2 = frt

end do iloop

1329

1462

end do iloop

end do icloop

1330

1463

end do icloop

x(l) = rt

1331

1464

x(l) = rt

infer(l) = jk

1332

1465

infer(l) = jk

l = l+1

1333

1466

l = l+1

end do rloop

1334

1467

end do rloop

1335

1468

contains

1336

1469

contains

subroutine trans_rt()

1337

1470

subroutine trans_rt()

tem = rten*eps1

1338

1471

tem = rten*eps1

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

1339

1472

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

rt = rt+tem

1340

1473

rt = rt+tem

d = (h+tem)*d/h

1341

1474

d = (h+tem)*d/h

h = h+tem

1342

1475

h = h+tem

end subroutine trans_rt

1343

1476

end subroutine trans_rt

1344

1477

subroutine deflated_work(errcode)

1345

1478

subroutine deflated_work(errcode)

! errcode=0 => no errors

1346

1479

! errcode=0 => no errors

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

1347

1480

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

integer :: errcode,flag

1348

1481

integer :: errcode,flag

1349

1482

flag=1

1350

1483

flag=1

loop1: do while(flag==1)

1351

1484

loop1: do while(flag==1)

errcode=0

1352

1485

errcode=0

jk = jk+1

1353

1486

jk = jk+1

if (jk .gt. itmax) then

1354

1487

if (jk .gt. itmax) then

ier=33

1355

1488

ier=33

errcode=1

1356

1489

errcode=1

return

1357

1490

return

end if

1358

1491

end if

frt = f(rt)

1359

1492

frt = f(rt)

fprt = frt

1360

1493

fprt = frt

if (l /= 1) then

1361

1494

if (l /= 1) then

lm1 = l-1

1362

1495

lm1 = l-1

do i=1,lm1

1363

1496

do i=1,lm1

tem = rt - x(i)

1364

1497

tem = rt - x(i)

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

1365

1498

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

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

1366

1499

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

call trans_rt()

1367

1500

call trans_rt()

cycle loop1

1368

1501

cycle loop1

end if

1369

1502

end if

fprt = fprt/tem

1370

1503

fprt = fprt/tem

end do

1371

1504

end do

end if

1372

1505

end if

flag=0

1373

1506

flag=0

end do loop1

1374

1507

end do loop1

1375

1508

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

1376

1509

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

errcode=1

1377

1510

errcode=1

return

1378

1511

return

end if

1379

1512

end if

1380

1513

end subroutine deflated_work

1381

1514

end subroutine deflated_work

1382

1515

end subroutine

1383

1516

end subroutine

1384

1517

1385

1518

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

1386

1519

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

!

1387

1520

!

! Integration

1388

1521

! Integration

!

1389

1522

!

! Only double precision coded atm

1390

1523

! Only double precision coded atm

!

1391

1524

!

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

1392

1525

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

1393

1526

1394

1527

subroutine fvn_d_gauss_legendre(n,qx,qw)

1395

1528

subroutine fvn_d_gauss_legendre(n,qx,qw)

!

1396

1529

!

! This routine compute the n Gauss Legendre abscissas and weights

1397

1530

! This routine compute the n Gauss Legendre abscissas and weights

! Adapted from Numerical Recipes routine gauleg

1398

1531

! Adapted from Numerical Recipes routine gauleg

!

1399

1532

!

! n (in) : number of points

1400

1533

! n (in) : number of points

! qx(out) : abscissas

1401

1534

! qx(out) : abscissas

! qw(out) : weights

1402

1535

! qw(out) : weights

!

1403

1536

!

implicit none

1404

1537

implicit none

double precision,parameter :: pi=3.141592653589793d0

1405

1538

double precision,parameter :: pi=3.141592653589793d0

integer, intent(in) :: n

1406

1539

integer, intent(in) :: n

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

1407

1540

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

1408

1541

integer :: m,i,j

1409

1542

integer :: m,i,j

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

1410

1543

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

m=(n+1)/2

1411

1544

m=(n+1)/2

do i=1,m

1412

1545

do i=1,m

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

1413

1546

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

iloop: do

1414

1547

iloop: do

p1=1.d0

1415

1548

p1=1.d0

p2=0.d0

1416

1549

p2=0.d0

do j=1,n

1417

1550

do j=1,n

p3=p2

1418

1551

p3=p2

p2=p1

1419

1552

p2=p1

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

1420

1553

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

end do

1421

1554

end do

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

1422

1555

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

z1=z

1423

1556

z1=z

z=z1-p1/pp

1424

1557

z=z1-p1/pp

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

1425

1558

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

exit iloop

1426

1559

exit iloop

end if

1427

1560

end if

end do iloop

1428

1561

end do iloop

qx(i)=-z

1429

1562

qx(i)=-z

qx(n+1-i)=z

1430

1563

qx(n+1-i)=z

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

1431

1564

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

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

1432

1565

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

end do

1433

1566

end do

end subroutine

1434

1567

end subroutine

1435

1568

1436

1569

1437

1570

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

1438

1571

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

!

1439

1572

!

! This is a simple non adaptative integration routine

1440

1573

! This is a simple non adaptative integration routine

! using n gauss legendre abscissas and weights

1441

1574

! using n gauss legendre abscissas and weights

!

1442

1575

!

! f(in) : the function to integrate

1443

1576

! f(in) : the function to integrate

! a(in) : lower bound

1444

1577

! a(in) : lower bound

! b(in) : higher bound

1445

1578

! b(in) : higher bound

! n(in) : number of gauss legendre pairs

1446

1579

! n(in) : number of gauss legendre pairs

! res(out): the evaluation of the integral

1447

1580

! res(out): the evaluation of the integral

!

1448

1581

!

double precision,external :: f

1449

1582

double precision,external :: f

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

1450

1583

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

integer, intent(in):: n

1451

1584

integer, intent(in):: n

double precision, intent(out) :: res

1452

1585

double precision, intent(out) :: res

1453

1586

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

1454

1587

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

double precision :: xm,xr

1455

1588

double precision :: xm,xr

integer :: i

1456

1589

integer :: i

1457

1590

! First compute n gauss legendre abs and weight

1458

1591

! First compute n gauss legendre abs and weight

allocate(qx(n))

1459

1592

allocate(qx(n))

allocate(qw(n))

1460

1593

allocate(qw(n))

call fvn_d_gauss_legendre(n,qx,qw)

1461

1594

call fvn_d_gauss_legendre(n,qx,qw)

1462

1595

xm=0.5d0*(b+a)

1463

1596

xm=0.5d0*(b+a)

xr=0.5d0*(b-a)

1464

1597

xr=0.5d0*(b-a)

1465

1598

res=0.d0

1466

1599

res=0.d0

1467

1600

do i=1,n

1468

1601

do i=1,n

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

1469

1602

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

end do

1470

1603

end do

1471

1604

res=xr*res

1472

1605

res=xr*res

1473

1606

deallocate(qw)

1474

1607

deallocate(qw)

deallocate(qx)

1475

1608

deallocate(qx)

1476

1609

end subroutine

1477

1610

end subroutine

1478

1611

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

1479

1612

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

!

1480

1613

!

! Simple and double adaptative Gauss Kronrod integration based on

1481

1614

! Simple and double adaptative Gauss Kronrod integration based on

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

1482

1615

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

!

1483

1616

!

! Common parameters :

1484

1617

! Common parameters :

!

1485

1618

!

! key (in)

1486

1619

! key (in)

! epsabs

1487

1620

! epsabs

! epsrel

1488

1621

! epsrel

!

1489

1622

!

1490

1623

!

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

1491

1624

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

1492

1625

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

1493

1626

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

!

1494

1627

!

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

1495

1628

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

!

1496

1629

!

! f(in) : the function

1497

1630

! f(in) : the function

! a(in) : lower bound

1498

1631

! a(in) : lower bound

! b(in) : higher bound

1499

1632

! b(in) : higher bound

! epsabs(in) : desired absolute error

1500

1633

! epsabs(in) : desired absolute error

! epsrel(in) : desired relative error

1501

1634

! epsrel(in) : desired relative error

! key(in) : gauss kronrod rule

1502

1635

! key(in) : gauss kronrod rule

! 1: 7 - 15 points

1503

1636

! 1: 7 - 15 points

! 2: 10 - 21 points

1504

1637

! 2: 10 - 21 points

! 3: 15 - 31 points

1505

1638

! 3: 15 - 31 points

! 4: 20 - 41 points

1506

1639

! 4: 20 - 41 points

! 5: 25 - 51 points

1507

1640

! 5: 25 - 51 points

! 6: 30 - 61 points

1508

1641

! 6: 30 - 61 points

!

1509

1642

!

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

1510

1643

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

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

1511

1644

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

! behaviour as the imsl routine dqdag

1512

1645

! behaviour as the imsl routine dqdag

!

1513

1646

!

! res(out) : estimated integral value

1514

1647

! res(out) : estimated integral value

! abserr(out) : estimated absolute error

1515

1648

! abserr(out) : estimated absolute error

! ier(out) : error flag from quadpack routines

1516

1649

! ier(out) : error flag from quadpack routines

! 0 : no error

1517

1650

! 0 : no error

! 1 : maximum number of subdivisions allowed

1518

1651

! 1 : maximum number of subdivisions allowed

! has been achieved. one can allow more

1519

1652

! has been achieved. one can allow more

! subdivisions by increasing the value of

1520

1653

! subdivisions by increasing the value of

! limit (and taking the according dimension

1521

1654

! limit (and taking the according dimension

! adjustments into account). however, if

1522

1655

! adjustments into account). however, if

! this yield no improvement it is advised

1523

1656

! this yield no improvement it is advised

! to analyze the integrand in order to

1524

1657

! to analyze the integrand in order to

! determine the integration difficulaties.

1525

1658

! determine the integration difficulaties.

! if the position of a local difficulty can

1526

1659

! if the position of a local difficulty can

! be determined (i.e.singularity,

1527

1660

! be determined (i.e.singularity,

! discontinuity within the interval) one

1528

1661

! discontinuity within the interval) one

! will probably gain from splitting up the

1529

1662

! will probably gain from splitting up the

! interval at this point and calling the

1530

1663

! interval at this point and calling the

! integrator on the subranges. if possible,

1531

1664

! integrator on the subranges. if possible,

! an appropriate special-purpose integrator

1532

1665

! an appropriate special-purpose integrator

! should be used which is designed for

1533

1666

! should be used which is designed for

! handling the type of difficulty involved.

1534

1667

! handling the type of difficulty involved.

! 2 : the occurrence of roundoff error is

1535

1668

! 2 : the occurrence of roundoff error is

! detected, which prevents the requested

1536

1669

! detected, which prevents the requested

! tolerance from being achieved.

1537

1670

! tolerance from being achieved.

! 3 : extremely bad integrand behaviour occurs

1538

1671

! 3 : extremely bad integrand behaviour occurs

! at some points of the integration

1539

1672

! at some points of the integration

! interval.

1540

1673

! interval.

! 6 : the input is invalid, because

1541

1674

! 6 : the input is invalid, because

! (epsabs.le.0 and

1542

1675

! (epsabs.le.0 and

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

1543

1676

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

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

1544

1677

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

! result, abserr, neval, last are set

1545

1678

! result, abserr, neval, last are set

! to zero.

1546

1679

! to zero.

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

1547

1680

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

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

1548

1681

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

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

1549

1682

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

! work(limit+1) to b.

1550

1683

! work(limit+1) to b.

1551

1684

implicit none

1552

1685

implicit none

double precision, external :: f

1553

1686

double precision, external :: f

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

1554

1687

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

integer, intent(in) :: key

1555

1688

integer, intent(in) :: key

integer, intent(in) :: limit

1556

1689

integer, intent(in) :: limit

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

1557

1690

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

integer, intent(out) :: ier

1558

1691

integer, intent(out) :: ier

1559

1692

double precision, allocatable :: work(:)

1560

1693

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

1561

1694

integer, allocatable :: iwork(:)

integer :: lenw,neval,last

1562

1695

integer :: lenw,neval,last

1563

1696

! imsl value for limit is 500

1564

1697

! imsl value for limit is 500

lenw=limit*4

1565

1698

lenw=limit*4

1566

1699

allocate(iwork(limit))

1567

1700

allocate(iwork(limit))

allocate(work(lenw))

1568

1701

allocate(work(lenw))

1569

1702

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

1570

1703

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

1571

1704

deallocate(work)

1572

1705

deallocate(work)

deallocate(iwork)

1573

1706

deallocate(iwork)

1574

1707

end subroutine

1575

1708

end subroutine

1576

1709

1577

1710

1578

1711

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

1579

1712

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

!

1580

1713

!

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

1581

1714

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

!

1582

1715

!

! f(in) : the function

1583

1716

! f(in) : the function

! a(in) : lower bound

1584

1717

! a(in) : lower bound

! b(in) : higher bound

1585

1718

! b(in) : higher bound

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

1586

1719

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

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

1587

1720

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

! epsabs(in) : desired absolute error

1588

1721

! epsabs(in) : desired absolute error

! epsrel(in) : desired relative error

1589

1722

! epsrel(in) : desired relative error

! key(in) : gauss kronrod rule

1590

1723

! key(in) : gauss kronrod rule

! 1: 7 - 15 points

1591

1724

! 1: 7 - 15 points

! 2: 10 - 21 points

1592

1725

! 2: 10 - 21 points

! 3: 15 - 31 points

1593

1726

! 3: 15 - 31 points

! 4: 20 - 41 points

1594

1727

! 4: 20 - 41 points

! 5: 25 - 51 points

1595

1728

! 5: 25 - 51 points

! 6: 30 - 61 points

1596

1729

! 6: 30 - 61 points

!

1597

1730

!

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

1598

1731

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

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

1599

1732

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

! behaviour as the imsl routine dqdag

1600

1733

! behaviour as the imsl routine dqdag

!

1601

1734

!

! res(out) : estimated integral value

1602

1735

! res(out) : estimated integral value

! abserr(out) : estimated absolute error

1603

1736

! abserr(out) : estimated absolute error

! ier(out) : error flag from quadpack routines

1604

1737

! ier(out) : error flag from quadpack routines

! 0 : no error

1605

1738

! 0 : no error

! 1 : maximum number of subdivisions allowed

1606

1739

! 1 : maximum number of subdivisions allowed

! has been achieved. one can allow more

1607

1740

! has been achieved. one can allow more

! subdivisions by increasing the value of

1608

1741

! subdivisions by increasing the value of

! limit (and taking the according dimension

1609

1742

! limit (and taking the according dimension

! adjustments into account). however, if

1610

1743

! adjustments into account). however, if

! this yield no improvement it is advised

1611

1744

! this yield no improvement it is advised

! to analyze the integrand in order to

1612

1745

! to analyze the integrand in order to

! determine the integration difficulaties.

1613

1746

! determine the integration difficulaties.

! if the position of a local difficulty can

1614

1747

! if the position of a local difficulty can

! be determined (i.e.singularity,

1615

1748

! be determined (i.e.singularity,

! discontinuity within the interval) one

1616

1749

! discontinuity within the interval) one

! will probably gain from splitting up the

1617

1750

! will probably gain from splitting up the

! interval at this point and calling the

1618

1751

! interval at this point and calling the

! integrator on the subranges. if possible,

1619

1752

! integrator on the subranges. if possible,

! an appropriate special-purpose integrator

1620

1753

! an appropriate special-purpose integrator

! should be used which is designed for

1621

1754

! should be used which is designed for

! handling the type of difficulty involved.

1622

1755

! handling the type of difficulty involved.

! 2 : the occurrence of roundoff error is

1623

1756

! 2 : the occurrence of roundoff error is

! detected, which prevents the requested

1624

1757

! detected, which prevents the requested

! tolerance from being achieved.

1625

1758

! tolerance from being achieved.

! 3 : extremely bad integrand behaviour occurs

1626

1759

! 3 : extremely bad integrand behaviour occurs

! at some points of the integration

1627

1760

! at some points of the integration

! interval.

1628

1761

! interval.

! 6 : the input is invalid, because

1629

1762

! 6 : the input is invalid, because

! (epsabs.le.0 and

1630

1763

! (epsabs.le.0 and

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

1631

1764

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

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

1632

1765

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

! result, abserr, neval, last are set

1633

1766

! result, abserr, neval, last are set

! to zero.

1634

1767

! to zero.

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

1635

1768

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

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

1636

1769

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

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

1637

1770

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

! work(limit+1) to b.

1638

1771

! work(limit+1) to b.

1639

1772

implicit none

1640

1773

implicit none

double precision, external:: f,g,h

1641

1774

double precision, external:: f,g,h

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

1642

1775

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

integer, intent(in) :: key,limit

1643

1776

integer, intent(in) :: key,limit

integer, intent(out) :: ier

1644

1777

integer, intent(out) :: ier

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

1645

1778

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

1646

1779

1647

1780

double precision, allocatable :: work(:)

1648

1781

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

1649

1782

integer, allocatable :: iwork(:)

integer :: lenw,neval,last

1650

1783

integer :: lenw,neval,last

1651

1784

! imsl value for limit is 500

1652

1785

! imsl value for limit is 500

lenw=limit*4

1653

1786

lenw=limit*4

allocate(work(lenw))

1654

1787

allocate(work(lenw))

allocate(iwork(limit))

1655

1788

allocate(iwork(limit))

1656

1789

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

1657

1790

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

1658

1791

deallocate(iwork)

1659

1792

deallocate(iwork)

deallocate(work)

1660

1793

deallocate(work)

end subroutine

1661

1794

end subroutine

1662

1795

1663

1796

1664

1797

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

1665

1798

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

!

1666

1799

!

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

1667

1800

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

! given x value

1668

1801

! given x value

!

1669

1802

!

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

1670

1803

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

!

1671

1804

!

! f(in) : the function

1672

1805

! f(in) : the function

! x(in) : x

1673

1806

! x(in) : x

! a(in) : lower bound

1674

1807

! a(in) : lower bound

! b(in) : higher bound

1675

1808

! b(in) : higher bound

! epsabs(in) : desired absolute error

1676

1809

! epsabs(in) : desired absolute error

! epsrel(in) : desired relative error

1677

1810

! epsrel(in) : desired relative error

! key(in) : gauss kronrod rule

1678

1811

! key(in) : gauss kronrod rule

! 1: 7 - 15 points

1679

1812

! 1: 7 - 15 points

! 2: 10 - 21 points

1680

1813

! 2: 10 - 21 points

! 3: 15 - 31 points

1681

1814

! 3: 15 - 31 points

! 4: 20 - 41 points

1682

1815

! 4: 20 - 41 points

! 5: 25 - 51 points

1683

1816

! 5: 25 - 51 points

! 6: 30 - 61 points

1684

1817

! 6: 30 - 61 points

!

1685

1818

!

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

1686

1819

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

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

1687

1820

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

! behaviour as the imsl routine dqdag

1688

1821

! behaviour as the imsl routine dqdag

!

1689

1822

!

! res(out) : estimated integral value

1690

1823

! res(out) : estimated integral value

! abserr(out) : estimated absolute error

1691

1824

! abserr(out) : estimated absolute error

! ier(out) : error flag from quadpack routines

1692

1825

! ier(out) : error flag from quadpack routines

! 0 : no error

1693

1826

! 0 : no error

! 1 : maximum number of subdivisions allowed

1694

1827

! 1 : maximum number of subdivisions allowed

! has been achieved. one can allow more

1695

1828

! has been achieved. one can allow more

! subdivisions by increasing the value of

1696

1829

! subdivisions by increasing the value of

! limit (and taking the according dimension

1697

1830

! limit (and taking the according dimension

! adjustments into account). however, if

1698

1831

! adjustments into account). however, if

! this yield no improvement it is advised

1699

1832

! this yield no improvement it is advised

! to analyze the integrand in order to

1700

1833

! to analyze the integrand in order to

! determine the integration difficulaties.

1701

1834

! determine the integration difficulaties.

! if the position of a local difficulty can

1702

1835

! if the position of a local difficulty can

! be determined (i.e.singularity,

1703

1836

! be determined (i.e.singularity,

! discontinuity within the interval) one

1704

1837

! discontinuity within the interval) one

! will probably gain from splitting up the

1705

1838

! will probably gain from splitting up the

! interval at this point and calling the

1706

1839

! interval at this point and calling the

! integrator on the subranges. if possible,

1707

1840

! integrator on the subranges. if possible,

! an appropriate special-purpose integrator

1708

1841

! an appropriate special-purpose integrator

! should be used which is designed for

1709

1842

! should be used which is designed for

! handling the type of difficulty involved.

1710

1843

! handling the type of difficulty involved.

! 2 : the occurrence of roundoff error is

1711

1844

! 2 : the occurrence of roundoff error is

! detected, which prevents the requested

1712

1845

! detected, which prevents the requested

! tolerance from being achieved.

1713

1846

! tolerance from being achieved.

! 3 : extremely bad integrand behaviour occurs

1714

1847

! 3 : extremely bad integrand behaviour occurs

! at some points of the integration

1715

1848

! at some points of the integration

! interval.

1716

1849

! interval.

! 6 : the input is invalid, because

1717

1850

! 6 : the input is invalid, because

! (epsabs.le.0 and

1718

1851

! (epsabs.le.0 and

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

1719

1852

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

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

1720

1853

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

! result, abserr, neval, last are set

1721

1854

! result, abserr, neval, last are set

! to zero.

1722

1855

! to zero.

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

1723

1856

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

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

1724

1857

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

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

1725

1858

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

! work(limit+1) to b.

1726

1859

! work(limit+1) to b.

1727

1860

implicit none

1728

1861

implicit none

double precision, external:: f

1729

1862

double precision, external:: f

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

1730

1863

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

integer, intent(in) :: key,limit

1731

1864

integer, intent(in) :: key,limit

integer, intent(out) :: ier

1732

1865

integer, intent(out) :: ier

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

1733

1866

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

1734

1867

1735

1868

double precision, allocatable :: work(:)

1736

1869

double precision, allocatable :: work(:)

integer, allocatable :: iwork(:)

1737

1870

integer, allocatable :: iwork(:)

integer :: lenw,neval,last

1738

1871

integer :: lenw,neval,last

1739

1872

! imsl value for limit is 500

1740

1873

! imsl value for limit is 500

lenw=limit*4

1741

1874

lenw=limit*4

allocate(work(lenw))

1742

1875

allocate(work(lenw))

allocate(iwork(limit))

1743

1876

allocate(iwork(limit))

1744

1877

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

1745

1878

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

1746

1879

deallocate(iwork)

1747

1880

deallocate(iwork)

deallocate(work)

1748

1881

deallocate(work)

end subroutine

1749

1882

end subroutine

1750

1883

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

1751

1884

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

GITLAB

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