lxsd / fvn | GitLab

Browse Code »

Commit d1b0528b40d1b22a5e7f65da95939e13adc95a51

Authored by daniau 2007-06-25 10:00:50 +0200

1 parent 7c27b9e1d2

Exists in master and in 3 other branches

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

Showing 2 changed files with 3558 additions and 0 deletions Inline Diff

stable/fvnlib.f90
testing/fvnlib.f90

stable/fvnlib.f90

Diff comments View file @ d1b0528

File was created	1
	2	module fvn
	3	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	4	!
	5	! fvn : a f95 module replacement for some imsl routines
	6	! it uses lapack for linear algebra
	7	! it uses modified quadpack for integration
	8	!
	9	! William Daniau 2007
	10	! william.daniau@femto-st.fr
	11	!
	12	! Routines naming scheme :
	13	!
	14	! fvn_x_name
	15	! where x can be s : real
	16	! d : real double precision
	17	! c : complex
	18	! z : double complex
	19	!
	20	!
	21	! This piece of code is totally free! Do whatever you want with it. However
	22	! if you find it usefull it would be kind to give credits ;-) Nevertheless, you
	23	! may give credits to quadpack authors.
	24	!
	25	! Version 1.1
	26	!
	27	! TO DO LIST :
	28	! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm
	29	! eigenvalues are given with no particular order.
	30	! + Generic interface for fvn_x_name family -> fvn_name
	31	! + Make some parameters optional, status for example
	32	! + use f95 kinds "double complex" -> complex(kind=8)
	33	! + unify quadpack routines
	34	! + ...
	35	!
	36	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	37
	38	implicit none
	39	! All quadpack routines are private to the module
	40	private :: d1mach,dqag,dqag_2d_inner,dqag_2d_outer,dqage,dqage_2d_inner, &
	41	dqage_2d_outer,dqk15,dqk15_2d_inner,dqk15_2d_outer,dqk21,dqk21_2d_inner,dqk21_2d_outer, &
	42	dqk31,dqk31_2d_inner,dqk31_2d_outer,dqk41,dqk41_2d_inner,dqk41_2d_outer, &
	43	dqk51,dqk51_2d_inner,dqk51_2d_outer,dqk61,dqk61_2d_inner,dqk61_2d_outer,dqpsrt
	44
	45
	46	contains
	47
	48	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	49	!
	50	! Matrix inversion subroutines
	51	!
	52	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	53	subroutine fvn_s_matinv(d,a,inva,status)
	54	!
	55	! Matrix inversion of a real matrix using BLAS and LAPACK
	56	!
	57	! d (in) : matrix rank
	58	! a (in) : input matrix
	59	! inva (out) : inversed matrix
	60	! status (ou) : =0 if something failed
	61	!
	62	integer, intent(in) :: d
	63	real, intent(in) :: a(d,d)
	64	real, intent(out) :: inva(d,d)
	65	integer, intent(out) :: status
	66
	67	integer, allocatable :: ipiv(:)
	68	real, allocatable :: work(:)
	69	real twork(1)
	70	integer :: info
	71	integer :: lwork
	72
	73	status=1
	74
	75	allocate(ipiv(d))
	76	! copy a into inva using BLAS
	77	!call scopy(d*d,a,1,inva,1)
	78	inva(:,:)=a(:,:)
	79	! LU factorization using LAPACK
	80	call sgetrf(d,d,inva,d,ipiv,info)
	81	! if info is not equal to 0, something went wrong we exit setting status to 0
	82	if (info /= 0) then
	83	status=0
	84	deallocate(ipiv)
	85	return
	86	end if
	87	! we use the query fonction of xxxtri to obtain the optimal workspace size
	88	call sgetri(d,inva,d,ipiv,twork,-1,info)
	89	lwork=int(twork(1))
	90	allocate(work(lwork))
	91	! Matrix inversion using LAPACK
	92	call sgetri(d,inva,d,ipiv,work,lwork,info)
	93	! again if info is not equal to 0, we exit setting status to 0
	94	if (info /= 0) then
	95	status=0
	96	end if
	97	deallocate(work)
	98	deallocate(ipiv)
	99	end subroutine
	100
	101	subroutine fvn_d_matinv(d,a,inva,status)
	102	!
	103	! Matrix inversion of a double precision matrix using BLAS and LAPACK
	104	!
	105	! d (in) : matrix rank
	106	! a (in) : input matrix
	107	! inva (out) : inversed matrix
	108	! status (ou) : =0 if something failed
	109	!
	110	integer, intent(in) :: d
	111	double precision, intent(in) :: a(d,d)
	112	double precision, intent(out) :: inva(d,d)
	113	integer, intent(out) :: status
	114
	115	integer, allocatable :: ipiv(:)
	116	double precision, allocatable :: work(:)
	117	double precision :: twork(1)
	118	integer :: info
	119	integer :: lwork
	120
	121	status=1
	122
	123	allocate(ipiv(d))
	124	! copy a into inva using BLAS
	125	!call dcopy(d*d,a,1,inva,1)
	126	inva(:,:)=a(:,:)
	127	! LU factorization using LAPACK
	128	call dgetrf(d,d,inva,d,ipiv,info)
	129	! if info is not equal to 0, something went wrong we exit setting status to 0
	130	if (info /= 0) then
	131	status=0
	132	deallocate(ipiv)
	133	return
	134	end if
	135	! we use the query fonction of xxxtri to obtain the optimal workspace size
	136	call dgetri(d,inva,d,ipiv,twork,-1,info)
	137	lwork=int(twork(1))
	138	allocate(work(lwork))
	139	! Matrix inversion using LAPACK
	140	call dgetri(d,inva,d,ipiv,work,lwork,info)
	141	! again if info is not equal to 0, we exit setting status to 0
	142	if (info /= 0) then
	143	status=0
	144	end if
	145	deallocate(work)
	146	deallocate(ipiv)
	147	end subroutine
	148
	149	subroutine fvn_c_matinv(d,a,inva,status)
	150	!
	151	! Matrix inversion of a complex matrix using BLAS and LAPACK
	152	!
	153	! d (in) : matrix rank
	154	! a (in) : input matrix
	155	! inva (out) : inversed matrix
	156	! status (ou) : =0 if something failed
	157	!
	158	integer, intent(in) :: d
	159	complex, intent(in) :: a(d,d)
	160	complex, intent(out) :: inva(d,d)
	161	integer, intent(out) :: status
	162
	163	integer, allocatable :: ipiv(:)
	164	complex, allocatable :: work(:)
	165	complex :: twork(1)
	166	integer :: info
	167	integer :: lwork
	168
	169	status=1
	170
	171	allocate(ipiv(d))
	172	! copy a into inva using BLAS
	173	!call ccopy(d*d,a,1,inva,1)
	174	inva(:,:)=a(:,:)
	175
	176	! LU factorization using LAPACK
	177	call cgetrf(d,d,inva,d,ipiv,info)
	178	! if info is not equal to 0, something went wrong we exit setting status to 0
	179	if (info /= 0) then
	180	status=0
	181	deallocate(ipiv)
	182	return
	183	end if
	184	! we use the query fonction of xxxtri to obtain the optimal workspace size
	185	call cgetri(d,inva,d,ipiv,twork,-1,info)
	186	lwork=int(twork(1))
	187	allocate(work(lwork))
	188	! Matrix inversion using LAPACK
	189	call cgetri(d,inva,d,ipiv,work,lwork,info)
	190	! again if info is not equal to 0, we exit setting status to 0
	191	if (info /= 0) then
	192	status=0
	193	end if
	194	deallocate(work)
	195	deallocate(ipiv)
	196	end subroutine
	197
	198	subroutine fvn_z_matinv(d,a,inva,status)
	199	!
	200	! Matrix inversion of a double complex matrix using BLAS and LAPACK
	201	!
	202	! d (in) : matrix rank
	203	! a (in) : input matrix
	204	! inva (out) : inversed matrix
	205	! status (ou) : =0 if something failed
	206	!
	207	integer, intent(in) :: d
	208	double complex, intent(in) :: a(d,d)
	209	double complex, intent(out) :: inva(d,d)
	210	integer, intent(out) :: status
	211
	212	integer, allocatable :: ipiv(:)
	213	double complex, allocatable :: work(:)
	214	double complex :: twork(1)
	215	integer :: info
	216	integer :: lwork
	217
	218	status=1
	219
	220	allocate(ipiv(d))
	221	! copy a into inva using BLAS
	222	!call zcopy(d*d,a,1,inva,1)
	223	inva(:,:)=a(:,:)
	224
	225	! LU factorization using LAPACK
	226	call zgetrf(d,d,inva,d,ipiv,info)
	227	! if info is not equal to 0, something went wrong we exit setting status to 0
	228	if (info /= 0) then
	229	status=0
	230	deallocate(ipiv)
	231	return
	232	end if
	233	! we use the query fonction of xxxtri to obtain the optimal workspace size
	234	call zgetri(d,inva,d,ipiv,twork,-1,info)
	235	lwork=int(twork(1))
	236	allocate(work(lwork))
	237	! Matrix inversion using LAPACK
	238	call zgetri(d,inva,d,ipiv,work,lwork,info)
	239	! again if info is not equal to 0, we exit setting status to 0
	240	if (info /= 0) then
	241	status=0
	242	end if
	243	deallocate(work)
	244	deallocate(ipiv)
	245	end subroutine
	246
	247
	248	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	249	!
	250	! Determinants
	251	!
	252	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	253	function fvn_s_det(d,a,status)
	254	!
	255	! Evaluate the determinant of a square matrix using lapack LU factorization
	256	!
	257	! d (in) : matrix rank
	258	! a (in) : The Matrix
	259	! status (out) : =0 if LU factorization failed
	260	!
	261	integer, intent(in) :: d
	262	real, intent(in) :: a(d,d)
	263	integer, intent(out) :: status
	264	real :: fvn_s_det
	265
	266	real, allocatable :: wc_a(:,:)
	267	integer, allocatable :: ipiv(:)
	268	integer :: info,i
	269
	270	status=1
	271	allocate(wc_a(d,d))
	272	allocate(ipiv(d))
	273	wc_a(:,:)=a(:,:)
	274	call sgetrf(d,d,wc_a,d,ipiv,info)
	275	if (info/= 0) then
	276	status=0
	277	fvn_s_det=0.e0
	278	deallocate(ipiv)
	279	deallocate(wc_a)
	280	return
	281	end if
	282	fvn_s_det=1.e0
	283	do i=1,d
	284	if (ipiv(i)==i) then
	285	fvn_s_det=fvn_s_det*wc_a(i,i)
	286	else
	287	fvn_s_det=-fvn_s_det*wc_a(i,i)
	288	end if
	289	end do
	290	deallocate(ipiv)
	291	deallocate(wc_a)
	292
	293	end function
	294
	295	function fvn_d_det(d,a,status)
	296	!
	297	! Evaluate the determinant of a square matrix using lapack LU factorization
	298	!
	299	! d (in) : matrix rank
	300	! a (in) : The Matrix
	301	! status (out) : =0 if LU factorization failed
	302	!
	303	integer, intent(in) :: d
	304	double precision, intent(in) :: a(d,d)
	305	integer, intent(out) :: status
	306	double precision :: fvn_d_det
	307
	308	double precision, allocatable :: wc_a(:,:)
	309	integer, allocatable :: ipiv(:)
	310	integer :: info,i
	311
	312	status=1
	313	allocate(wc_a(d,d))
	314	allocate(ipiv(d))
	315	wc_a(:,:)=a(:,:)
	316	call dgetrf(d,d,wc_a,d,ipiv,info)
	317	if (info/= 0) then
	318	status=0
	319	fvn_d_det=0.d0
	320	deallocate(ipiv)
	321	deallocate(wc_a)
	322	return
	323	end if
	324	fvn_d_det=1.d0
	325	do i=1,d
	326	if (ipiv(i)==i) then
	327	fvn_d_det=fvn_d_det*wc_a(i,i)
	328	else
	329	fvn_d_det=-fvn_d_det*wc_a(i,i)
	330	end if
	331	end do
	332	deallocate(ipiv)
	333	deallocate(wc_a)
	334
	335	end function
	336
	337	function fvn_c_det(d,a,status) !
	338	! Evaluate the determinant of a square matrix using lapack LU factorization
	339	!
	340	! d (in) : matrix rank
	341	! a (in) : The Matrix
	342	! status (out) : =0 if LU factorization failed
	343	!
	344	integer, intent(in) :: d
	345	complex, intent(in) :: a(d,d)
	346	integer, intent(out) :: status
	347	complex :: fvn_c_det
	348
	349	complex, allocatable :: wc_a(:,:)
	350	integer, allocatable :: ipiv(:)
	351	integer :: info,i
	352
	353	status=1
	354	allocate(wc_a(d,d))
	355	allocate(ipiv(d))
	356	wc_a(:,:)=a(:,:)
	357	call cgetrf(d,d,wc_a,d,ipiv,info)
	358	if (info/= 0) then
	359	status=0
	360	fvn_c_det=(0.e0,0.e0)
	361	deallocate(ipiv)
	362	deallocate(wc_a)
	363	return
	364	end if
	365	fvn_c_det=(1.e0,0.e0)
	366	do i=1,d
	367	if (ipiv(i)==i) then
	368	fvn_c_det=fvn_c_det*wc_a(i,i)
	369	else
	370	fvn_c_det=-fvn_c_det*wc_a(i,i)
	371	end if
	372	end do
	373	deallocate(ipiv)
	374	deallocate(wc_a)
	375
	376	end function
	377
	378	function fvn_z_det(d,a,status)
	379	!
	380	! Evaluate the determinant of a square matrix using lapack LU factorization
	381	!
	382	! d (in) : matrix rank
	383	! a (in) : The Matrix
	384	! det (out) : determinant
	385	! status (out) : =0 if LU factorization failed
	386	!
	387	integer, intent(in) :: d
	388	double complex, intent(in) :: a(d,d)
	389	integer, intent(out) :: status
	390	double complex :: fvn_z_det
	391
	392	double complex, allocatable :: wc_a(:,:)
	393	integer, allocatable :: ipiv(:)
	394	integer :: info,i
	395
	396	status=1
	397	allocate(wc_a(d,d))
	398	allocate(ipiv(d))
	399	wc_a(:,:)=a(:,:)
	400	call zgetrf(d,d,wc_a,d,ipiv,info)
	401	if (info/= 0) then
	402	status=0
	403	fvn_z_det=(0.d0,0.d0)
	404	deallocate(ipiv)
	405	deallocate(wc_a)
	406	return
	407	end if
	408	fvn_z_det=(1.d0,0.d0)
	409	do i=1,d
	410	if (ipiv(i)==i) then
	411	fvn_z_det=fvn_z_det*wc_a(i,i)
	412	else
	413	fvn_z_det=-fvn_z_det*wc_a(i,i)
	414	end if
	415	end do
	416	deallocate(ipiv)
	417	deallocate(wc_a)
	418
	419	end function
	420
	421	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	422	!
	423	! Condition test
	424	!
	425	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	426	! 1-norm
	427	! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm
	428	! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond
	429	!
	430	subroutine fvn_s_matcon(d,a,rcond,status)
	431	! Matrix condition (reciprocal of condition number)
	432	!
	433	! d (in) : matrix rank
	434	! a (in) : The Matrix
	435	! rcond (out) : guess what
	436	! status (out) : =0 if something went wrong
	437	!
	438	integer, intent(in) :: d
	439	real, intent(in) :: a(d,d)
	440	real, intent(out) :: rcond
	441	integer, intent(out) :: status
	442
	443	real, allocatable :: work(:)
	444	integer, allocatable :: iwork(:)
	445	real :: anorm
	446	real, allocatable :: wc_a(:,:) ! working copy of a
	447	integer :: info
	448	integer, allocatable :: ipiv(:)
	449
	450	real, external :: slange
	451
	452
	453	status=1
	454
	455	anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
	456
	457	allocate(wc_a(d,d))
	458	!call scopy(d*d,a,1,wc_a,1)
	459	wc_a(:,:)=a(:,:)
	460
	461	allocate(ipiv(d))
	462	call sgetrf(d,d,wc_a,d,ipiv,info)
	463	if (info /= 0) then
	464	status=0
	465	deallocate(ipiv)
	466	deallocate(wc_a)
	467	return
	468	end if
	469	allocate(work(4*d))
	470	allocate(iwork(d))
	471	call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
	472	if (info /= 0) then
	473	status=0
	474	end if
	475	deallocate(iwork)
	476	deallocate(work)
	477	deallocate(ipiv)
	478	deallocate(wc_a)
	479
	480	end subroutine
	481
	482	subroutine fvn_d_matcon(d,a,rcond,status)
	483	! Matrix condition (reciprocal of condition number)
	484	!
	485	! d (in) : matrix rank
	486	! a (in) : The Matrix
	487	! rcond (out) : guess what
	488	! status (out) : =0 if something went wrong
	489	!
	490	integer, intent(in) :: d
	491	double precision, intent(in) :: a(d,d)
	492	double precision, intent(out) :: rcond
	493	integer, intent(out) :: status
	494
	495	double precision, allocatable :: work(:)
	496	integer, allocatable :: iwork(:)
	497	double precision :: anorm
	498	double precision, allocatable :: wc_a(:,:) ! working copy of a
	499	integer :: info
	500	integer, allocatable :: ipiv(:)
	501
	502	double precision, external :: dlange
	503
	504
	505	status=1
	506
	507	anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
	508
	509	allocate(wc_a(d,d))
	510	!call dcopy(d*d,a,1,wc_a,1)
	511	wc_a(:,:)=a(:,:)
	512
	513	allocate(ipiv(d))
	514	call dgetrf(d,d,wc_a,d,ipiv,info)
	515	if (info /= 0) then
	516	status=0
	517	deallocate(ipiv)
	518	deallocate(wc_a)
	519	return
	520	end if
	521
	522	allocate(work(4*d))
	523	allocate(iwork(d))
	524	call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
	525	if (info /= 0) then
	526	status=0
	527	end if
	528	deallocate(iwork)
	529	deallocate(work)
	530	deallocate(ipiv)
	531	deallocate(wc_a)
	532
	533	end subroutine
	534
	535	subroutine fvn_c_matcon(d,a,rcond,status)
	536	! Matrix condition (reciprocal of condition number)
	537	!
	538	! d (in) : matrix rank
	539	! a (in) : The Matrix
	540	! rcond (out) : guess what
	541	! status (out) : =0 if something went wrong
	542	!
	543	integer, intent(in) :: d
	544	complex, intent(in) :: a(d,d)
	545	real, intent(out) :: rcond
	546	integer, intent(out) :: status
	547
	548	real, allocatable :: rwork(:)
	549	complex, allocatable :: work(:)
	550	integer, allocatable :: iwork(:)
	551	real :: anorm
	552	complex, allocatable :: wc_a(:,:) ! working copy of a
	553	integer :: info
	554	integer, allocatable :: ipiv(:)
	555
	556	real, external :: clange
	557
	558
	559	status=1
	560
	561	anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
	562
	563	allocate(wc_a(d,d))
	564	!call ccopy(d*d,a,1,wc_a,1)
	565	wc_a(:,:)=a(:,:)
	566
	567	allocate(ipiv(d))
	568	call cgetrf(d,d,wc_a,d,ipiv,info)
	569	if (info /= 0) then
	570	status=0
	571	deallocate(ipiv)
	572	deallocate(wc_a)
	573	return
	574	end if
	575	allocate(work(2*d))
	576	allocate(rwork(2*d))
	577	call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
	578	if (info /= 0) then
	579	status=0
	580	end if
	581	deallocate(rwork)
	582	deallocate(work)
	583	deallocate(ipiv)
	584	deallocate(wc_a)
	585	end subroutine
	586
	587	subroutine fvn_z_matcon(d,a,rcond,status)
	588	! Matrix condition (reciprocal of condition number)
	589	!
	590	! d (in) : matrix rank
	591	! a (in) : The Matrix
	592	! rcond (out) : guess what
	593	! status (out) : =0 if something went wrong
	594	!
	595	integer, intent(in) :: d
	596	double complex, intent(in) :: a(d,d)
	597	double precision, intent(out) :: rcond
	598	integer, intent(out) :: status
	599
	600	double complex, allocatable :: work(:)
	601	double precision, allocatable :: rwork(:)
	602	double precision :: anorm
	603	double complex, allocatable :: wc_a(:,:) ! working copy of a
	604	integer :: info
	605	integer, allocatable :: ipiv(:)
	606
	607	double precision, external :: zlange
	608
	609
	610	status=1
	611
	612	anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
	613
	614	allocate(wc_a(d,d))
	615	!call zcopy(d*d,a,1,wc_a,1)
	616	wc_a(:,:)=a(:,:)
	617
	618	allocate(ipiv(d))
	619	call zgetrf(d,d,wc_a,d,ipiv,info)
	620	if (info /= 0) then
	621	status=0
	622	deallocate(ipiv)
	623	deallocate(wc_a)
	624	return
	625	end if
	626
	627	allocate(work(2*d))
	628	allocate(rwork(2*d))
	629	call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
	630	if (info /= 0) then
	631	status=0
	632	end if
	633	deallocate(rwork)
	634	deallocate(work)
	635	deallocate(ipiv)
	636	deallocate(wc_a)
	637	end subroutine
	638
	639	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	640	!
	641	! Valeurs propres/ Vecteurs propre
	642	!
	643	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	644
	645	subroutine fvn_s_matev(d,a,evala,eveca,status)
	646	!
	647	! integer d (in) : matrice rank
	648	! real a(d,d) (in) : The Matrix
	649	! complex evala(d) (out) : eigenvalues
	650	! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
	651	! integer (out) : status =0 if something went wrong
	652	!
	653	! interfacing Lapack routine SGEEV
	654
	655	integer, intent(in) :: d
	656	real, intent(in) :: a(d,d)
	657	complex, intent(out) :: evala(d)
	658	complex, intent(out) :: eveca(d,d)
	659	integer, intent(out) :: status
	660
	661	real, allocatable :: wc_a(:,:) ! a working copy of a
	662	integer :: info
	663	integer :: lwork
	664	real, allocatable :: wr(:),wi(:)
	665	real :: vl ! unused but necessary for the call
	666	real, allocatable :: vr(:,:)
	667	real, allocatable :: work(:)
	668	real :: twork(1)
	669	integer i
	670	integer j
	671
	672	! making a working copy of a
	673	allocate(wc_a(d,d))
	674	!call scopy(d*d,a,1,wc_a,1)
	675	wc_a(:,:)=a(:,:)
	676
	677	allocate(wr(d))
	678	allocate(wi(d))
	679	allocate(vr(d,d))
	680	! query optimal work size
	681	call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
	682	lwork=int(twork(1))
	683	allocate(work(lwork))
	684	call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
	685
	686	if (info /= 0) then
	687	status=0
	688	deallocate(work)
	689	deallocate(vr)
	690	deallocate(wi)
	691	deallocate(wr)
	692	deallocate(wc_a)
	693	return
	694	end if
	695
	696	! now fill in the results
	697	i=1
	698	do while(i<=d)
	699	evala(i)=cmplx(wr(i),wi(i))
	700	if (wi(i) == 0.) then ! eigenvalue is real
	701	eveca(:,i)=cmplx(vr(:,i),0.)
	702	else ! eigenvalue is complex
	703	evala(i+1)=cmplx(wr(i+1),wi(i+1))
	704	eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))
	705	eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))
	706	i=i+1
	707	end if
	708	i=i+1
	709	enddo
	710	deallocate(work)
	711	deallocate(vr)
	712	deallocate(wi)
	713	deallocate(wr)
	714	deallocate(wc_a)
	715
	716	end subroutine
	717
	718	subroutine fvn_d_matev(d,a,evala,eveca,status)
	719	!
	720	! integer d (in) : matrice rank
	721	! double precision a(d,d) (in) : The Matrix
	722	! double complex evala(d) (out) : eigenvalues
	723	! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
	724	! integer (out) : status =0 if something went wrong
	725	!
	726	! interfacing Lapack routine DGEEV
	727	integer, intent(in) :: d
	728	double precision, intent(in) :: a(d,d)
	729	double complex, intent(out) :: evala(d)
	730	double complex, intent(out) :: eveca(d,d)
	731	integer, intent(out) :: status
	732
	733	double precision, allocatable :: wc_a(:,:) ! a working copy of a
	734	integer :: info
	735	integer :: lwork
	736	double precision, allocatable :: wr(:),wi(:)
	737	double precision :: vl ! unused but necessary for the call
	738	double precision, allocatable :: vr(:,:)
	739	double precision, allocatable :: work(:)
	740	double precision :: twork(1)
	741	integer i
	742	integer j
	743
	744	! making a working copy of a
	745	allocate(wc_a(d,d))
	746	!call dcopy(d*d,a,1,wc_a,1)
	747	wc_a(:,:)=a(:,:)
	748
	749	allocate(wr(d))
	750	allocate(wi(d))
	751	allocate(vr(d,d))
	752	! query optimal work size
	753	call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
	754	lwork=int(twork(1))
	755	allocate(work(lwork))
	756	call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
	757
	758	if (info /= 0) then
	759	status=0
	760	deallocate(work)
	761	deallocate(vr)
	762	deallocate(wi)
	763	deallocate(wr)
	764	deallocate(wc_a)
	765	return
	766	end if
	767
	768	! now fill in the results
	769	i=1
	770	do while(i<=d)
	771	evala(i)=dcmplx(wr(i),wi(i))
	772	if (wi(i) == 0.) then ! eigenvalue is real
	773	eveca(:,i)=dcmplx(vr(:,i),0.)
	774	else ! eigenvalue is complex
	775	evala(i+1)=dcmplx(wr(i+1),wi(i+1))
	776	eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))
	777	eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))
	778	i=i+1
	779	end if
	780	i=i+1
	781	enddo
	782
	783	deallocate(work)
	784	deallocate(vr)
	785	deallocate(wi)
	786	deallocate(wr)
	787	deallocate(wc_a)
	788
	789	end subroutine
	790
	791	subroutine fvn_c_matev(d,a,evala,eveca,status)
	792	!
	793	! integer d (in) : matrice rank
	794	! complex a(d,d) (in) : The Matrix
	795	! complex evala(d) (out) : eigenvalues
	796	! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
	797	! integer (out) : status =0 if something went wrong
	798	!
	799	! interfacing Lapack routine CGEEV
	800
	801	integer, intent(in) :: d
	802	complex, intent(in) :: a(d,d)
	803	complex, intent(out) :: evala(d)
	804	complex, intent(out) :: eveca(d,d)
	805	integer, intent(out) :: status
	806
	807	complex, allocatable :: wc_a(:,:) ! a working copy of a
	808	integer :: info
	809	integer :: lwork
	810	complex, allocatable :: work(:)
	811	complex :: twork(1)
	812	real, allocatable :: rwork(:)
	813	complex :: vl ! unused but necessary for the call
	814
	815	status=1
	816
	817	! making a working copy of a
	818	allocate(wc_a(d,d))
	819	!call ccopy(d*d,a,1,wc_a,1)
	820	wc_a(:,:)=a(:,:)
	821
	822
	823	! query optimal work size
	824	call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
	825	lwork=int(twork(1))
	826	allocate(work(lwork))
	827	allocate(rwork(2*d))
	828	call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
	829
	830	if (info /= 0) then
	831	status=0
	832	end if
	833	deallocate(rwork)
	834	deallocate(work)
	835	deallocate(wc_a)
	836
	837	end subroutine
	838
	839	subroutine fvn_z_matev(d,a,evala,eveca,status)
	840	!
	841	! integer d (in) : matrice rank
	842	! double complex a(d,d) (in) : The Matrix
	843	! double complex evala(d) (out) : eigenvalues
	844	! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
	845	! integer (out) : status =0 if something went wrong
	846	!
	847	! interfacing Lapack routine ZGEEV
	848
	849	integer, intent(in) :: d
	850	double complex, intent(in) :: a(d,d)
	851	double complex, intent(out) :: evala(d)
	852	double complex, intent(out) :: eveca(d,d)
	853	integer, intent(out) :: status
	854
	855	double complex, allocatable :: wc_a(:,:) ! a working copy of a
	856	integer :: info
	857	integer :: lwork
	858	double complex, allocatable :: work(:)
	859	double complex :: twork(1)
	860	double precision, allocatable :: rwork(:)
	861	double complex :: vl ! unused but necessary for the call
	862
	863	status=1
	864
	865	! making a working copy of a
	866	allocate(wc_a(d,d))
	867	!call zcopy(d*d,a,1,wc_a,1)
	868	wc_a(:,:)=a(:,:)
	869
	870	! query optimal work size
	871	call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
	872	lwork=int(twork(1))
	873	allocate(work(lwork))
	874	allocate(rwork(2*d))
	875	call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
	876
	877	if (info /= 0) then
	878	status=0
	879	end if
	880	deallocate(rwork)
	881	deallocate(work)
	882	deallocate(wc_a)
	883
	884	end subroutine
	885
	886	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	887	!
	888	! Akima spline interpolation and spline evaluation
	889	!
	890	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	891
	892	! Single precision
	893	subroutine fvn_s_akima(n,x,y,br,co)
	894	implicit none
	895	integer, intent(in) :: n
	896	real, intent(in) :: x(n)
	897	real, intent(in) :: y(n)
	898	real, intent(out) :: br(n)
	899	real, intent(out) :: co(4,n)
	900
	901	real, allocatable :: var(:),z(:)
	902	real :: wi_1,wi
	903	integer :: i
	904	real :: dx,a,b
	905
	906	! br is just a copy of x
	907	br(:)=x(:)
	908
	909	allocate(var(n))
	910	allocate(z(n))
	911	! evaluate the variations
	912	do i=1, n-1
	913	var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
	914	end do
	915	var(n+2)=2.e0*var(n+1)-var(n)
	916	var(n+3)=2.e0*var(n+2)-var(n+1)
	917	var(2)=2.e0*var(3)-var(4)
	918	var(1)=2.e0*var(2)-var(3)
	919
	920	do i = 1, n
	921	wi_1=abs(var(i+3)-var(i+2))
	922	wi=abs(var(i+1)-var(i))
	923	if ((wi_1+wi).eq.0.e0) then
	924	z(i)=(var(i+2)+var(i+1))/2.e0
	925	else
	926	z(i)=(wi_1var(i+1)+wivar(i+2))/(wi_1+wi)
	927	end if
	928	end do
	929
	930	do i=1, n-1
	931	dx=x(i+1)-x(i)
	932	a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd
	933	b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
	934	co(1,i)=y(i)
	935	co(2,i)=z(i)
	936	!co(3,i)=-(a-3.b)/dx*2 ! méthode wd
	937	!co(4,i)=(a-2.b)/dx*3 ! méthode wd
	938	co(3,i)=(3.e0var(i+2)-2.e0z(i)-z(i+1))/dx ! méthode JP Moreau
	939	co(4,i)=(z(i)+z(i+1)-2.e0var(i+2))/dx*2 !
	940	! les coefficients donnés par imsl sont co(3,i)2 et co(4,i)6
	941	! etrangement la fonction csval corrige et donne la bonne valeur ...
	942	end do
	943	co(1,n)=y(n)
	944	co(2,n)=z(n)
	945	co(3,n)=0.e0
	946	co(4,n)=0.e0
	947
	948	deallocate(z)
	949	deallocate(var)
	950
	951	end subroutine
	952
	953	! Double precision
	954	subroutine fvn_d_akima(n,x,y,br,co)
	955
	956	implicit none
	957	integer, intent(in) :: n
	958	double precision, intent(in) :: x(n)
	959	double precision, intent(in) :: y(n)
	960	double precision, intent(out) :: br(n)
	961	double precision, intent(out) :: co(4,n)
	962
	963	double precision, allocatable :: var(:),z(:)
	964	double precision :: wi_1,wi
	965	integer :: i
	966	double precision :: dx,a,b
	967
	968	! br is just a copy of x
	969	br(:)=x(:)
	970
	971	allocate(var(n))
	972	allocate(z(n))
	973	! evaluate the variations
	974	do i=1, n-1
	975	var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
	976	end do
	977	var(n+2)=2.d0*var(n+1)-var(n)
	978	var(n+3)=2.d0*var(n+2)-var(n+1)
	979	var(2)=2.d0*var(3)-var(4)
	980	var(1)=2.d0*var(2)-var(3)
	981
	982	do i = 1, n
	983	wi_1=dabs(var(i+3)-var(i+2))
	984	wi=dabs(var(i+1)-var(i))
	985	if ((wi_1+wi).eq.0.d0) then
	986	z(i)=(var(i+2)+var(i+1))/2.d0
	987	else
	988	z(i)=(wi_1var(i+1)+wivar(i+2))/(wi_1+wi)
	989	end if
	990	end do
	991
	992	do i=1, n-1
	993	dx=x(i+1)-x(i)
	994	a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd
	995	b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
	996	co(1,i)=y(i)
	997	co(2,i)=z(i)
	998	!co(3,i)=-(a-3.b)/dx*2 ! méthode wd
	999	!co(4,i)=(a-2.b)/dx*3 ! méthode wd
	1000	co(3,i)=(3.d0var(i+2)-2.d0z(i)-z(i+1))/dx ! méthode JP Moreau
	1001	co(4,i)=(z(i)+z(i+1)-2.d0var(i+2))/dx*2 !
	1002	! les coefficients donnés par imsl sont co(3,i)2 et co(4,i)6
	1003	! etrangement la fonction csval corrige et donne la bonne valeur ...
	1004	end do
	1005	co(1,n)=y(n)
	1006	co(2,n)=z(n)
	1007	co(3,n)=0.d0
	1008	co(4,n)=0.d0
	1009
	1010	deallocate(z)
	1011	deallocate(var)
	1012
	1013	end subroutine
	1014
	1015	!
	1016	! Single precision spline evaluation
	1017	!
	1018	function fvn_s_spline_eval(x,n,br,co)
	1019	implicit none
	1020	real, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
	1021	integer, intent(in) :: n ! number of intervals
	1022	real, intent(in) :: br(n+1) ! breakpoints
	1023	real, intent(in) :: co(4,n+1) ! spline coeeficients
	1024	real :: fvn_s_spline_eval
	1025
	1026	integer :: i
	1027	real :: dx
	1028
	1029	if (x<=br(1)) then
	1030	i=1
	1031	else if (x>=br(n+1)) then
	1032	i=n
	1033	else
	1034	i=1
	1035	do while(x>=br(i))
	1036	i=i+1
	1037	end do
	1038	i=i-1
	1039	end if
	1040	dx=x-br(i)
	1041	fvn_s_spline_eval=co(1,i)+co(2,i)dx+co(3,i)dx*2+co(4,i)dx**3
	1042
	1043	end function
	1044
	1045	! Double precision spline evaluation
	1046	function fvn_d_spline_eval(x,n,br,co)
	1047	implicit none
	1048	double precision, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
	1049	integer, intent(in) :: n ! number of intervals
	1050	double precision, intent(in) :: br(n+1) ! breakpoints
	1051	double precision, intent(in) :: co(4,n+1) ! spline coeeficients
	1052	double precision :: fvn_d_spline_eval
	1053
	1054	integer :: i
	1055	double precision :: dx
	1056
	1057
	1058	if (x<=br(1)) then
	1059	i=1
	1060	else if (x>=br(n+1)) then
	1061	i=n
	1062	else
	1063	i=1
	1064	do while(x>=br(i))
	1065	i=i+1
	1066	end do
	1067	i=i-1
	1068	end if
	1069
	1070	dx=x-br(i)
	1071	fvn_d_spline_eval=co(1,i)+co(2,i)dx+co(3,i)dx*2+co(4,i)dx**3
	1072
	1073	end function
	1074
	1075
	1076	!
	1077	! Muller
	1078	!
	1079	!
	1080	!
	1081	! William Daniau 2007
	1082	!
	1083	! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f
	1084	! http://plato.asu.edu/ftp/other_software/muller.f
	1085	!
	1086	! it can be used as a replacement for imsl routine dzanly with minor changes
	1087	!
	1088	!-----------------------------------------------------------------------
	1089	!
	1090	! purpose - zeros of an analytic complex function
	1091	! using the muller method with deflation
	1092	!
	1093	! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,
	1094	! infer,ier)
	1095	!
	1096	! arguments f - a complex function subprogram, f(z), written
	1097	! by the user specifying the equation whose
	1098	! roots are to be found. f must appear in
	1099	! an external statement in the calling pro-
	1100	! gram.
	1101	! eps - 1st stopping criterion. let fp(z)=f(z)/p
	1102	! where p = (z-z(1))(z-z(2)),,,*(z-z(k-1))
	1103	! and z(1),...,z(k-1) are previously found
	1104	! roots. if ((cdabs(f(z)).le.eps) .and.
	1105	! (cdabs(fp(z)).le.eps)), then z is accepted
	1106	! as a root. (input)
	1107	! eps1 - 2nd stopping criterion. a root is accepted
	1108	! if two successive approximations to a given
	1109	! root agree within eps1. (input)
	1110	! note. if either or both of the stopping
	1111	! criteria are fulfilled, the root is
	1112	! accepted.
	1113	! kn - the number of known roots which must be stored
	1114	! in x(1),...,x(kn), prior to entry to muller
	1115	! nguess - the number of initial guesses provided. these
	1116	! guesses must be stored in x(kn+1),...,
	1117	! x(kn+nguess). nguess must be set equal
	1118	! to zero if no guesses are provided. (input)
	1119	! n - the number of new roots to be found by
	1120	! muller (input)
	1121	! x - a complex vector of length kn+n. x(1),...,
	1122	! x(kn) on input must contain any known
	1123	! roots. x(kn+1),..., x(kn+n) on input may,
	1124	! on user option, contain initial guesses for
	1125	! the n new roots which are to be computed.
	1126	! if the user does not provide an initial
	1127	! guess, zero is used.
	1128	! on output, x(kn+1),...,x(kn+n) contain the
	1129	! approximate roots found by muller.
	1130	! itmax - the maximum allowable number of iterations
	1131	! per root (input)
	1132	! infer - an integer vector of length kn+n. on
	1133	! output infer(j) contains the number of
	1134	! iterations used in finding the j-th root
	1135	! when convergence was achieved. if
	1136	! convergence was not obtained in itmax
	1137	! iterations, infer(j) will be greater than
	1138	! itmax (output).
	1139	! ier - error parameter (output)
	1140	! warning error
	1141	! ier = 33 indicates failure to converge with-
	1142	! in itmax iterations for at least one of
	1143	! the (n) new roots.
	1144	!
	1145	!
	1146	! remarks muller always returns the last approximation for root j
	1147	! in x(j). if the convergence criterion is satisfied,
	1148	! then infer(j) is less than or equal to itmax. if the
	1149	! convergence criterion is not satisified, then infer(j)
	1150	! is set to either itmax+1 or itmax+k, with k greater
	1151	! than 1. infer(j) = itmax+1 indicates that muller did
	1152	! not obtain convergence in the allowed number of iter-
	1153	! ations. in this case, the user may wish to set itmax
	1154	! to a larger value. infer(j) = itmax+k means that con-
	1155	! vergence was obtained (on iteration k) for the defla-
	1156	! ted function
	1157	! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))
	1158	!
	1159	! but failed for f(z). in this case, better initial
	1160	! guesses might help or, it might be necessary to relax
	1161	! the convergence criterion.
	1162	!
	1163	!-----------------------------------------------------------------------
	1164	!
	1165	subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
	1166	implicit none
	1167	double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq
	1168	double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
	1169	tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
	1170	zero,p1,one,four,p5
	1171
	1172	double complex, external :: f
	1173	integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
	1174	knpng,jk,ick,nn,lm1,errcode
	1175	double complex :: x(kn+n)
	1176	integer :: infer(kn+n)
	1177
	1178
	1179	data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &
	1180	one/(1.0,0.0)/,four/(4.0,0.0)/, &
	1181	p5/(0.5,0.0)/, &
	1182	rzero/0.0/,rten/10.0/,rhun/100.0/, &
	1183	ax/0.1/,ickmax/3/,rp01/0.01/
	1184
	1185	ier = 0
	1186	if (n .lt. 1) then ! What the hell are doing here then ...
	1187	return
	1188	end if
	1189	!eps1 = rten **(-nsig)
	1190	eps1 = min(eps1,rp01)
	1191
	1192	knp1 = kn+1
	1193	knpn = kn+n
	1194	knpng = kn+nguess
	1195	do i=1,knpn
	1196	infer(i) = 0
	1197	if (i .gt. knpng) x(i) = zero
	1198	end do
	1199	l= knp1
	1200
	1201	ic=0
	1202	rloop: do while (l<=knpn) ! Main loop over new roots
	1203	jk = 0
	1204	ick = 0
	1205	xl = x(l)
	1206	icloop: do
	1207	ic = 0
	1208	h = ax
	1209	h = p1*h
	1210	if (cdabs(xl) .gt. ax) h = p1*xl
	1211	! first three points are
	1212	! xl+h, xl-h, xl
	1213	rt = xl+h
	1214	call deflated_work(errcode)
	1215	if (errcode == 1) then
	1216	exit icloop
	1217	end if
	1218
	1219	z0 = fprt
	1220	y0 = frt
	1221	x0 = rt
	1222	rt = xl-h
	1223	call deflated_work(errcode)
	1224	if (errcode == 1) then
	1225	exit icloop
	1226	end if
	1227
	1228	z1 = fprt
	1229	y1 = frt
	1230	h = xl-rt
	1231	d = h/(rt-x0)
	1232	rt = xl
	1233
	1234	call deflated_work(errcode)
	1235	if (errcode == 1) then
	1236	exit icloop
	1237	end if
	1238
	1239
	1240	z2 = fprt
	1241	y2 = frt
	1242	! begin main algorithm
	1243	iloop: do
	1244	dd = one + d
	1245	t1 = z0dd
	1246	t2 = z1dddd
	1247	xx = z2*dd
	1248	t3 = z2*d
	1249	bi = t1-t2+xx+t3
	1250	den = bibi-four(xxt1-t3(t2-xx))
	1251	! use denominator of maximum amplitude
	1252	t1 = cdsqrt(den)
	1253	qz = rhun*max(cdabs(bi),cdabs(t1))
	1254	t2 = bi + t1
	1255	tpq = cdabs(t2)+qz
	1256	if (tpq .eq. qz) t2 = zero
	1257	t3 = bi - t1
	1258	tpq = cdabs(t3) + qz
	1259	if (tpq .eq. qz) t3 = zero
	1260	den = t2
	1261	qz = cdabs(t3)-cdabs(t2)
	1262	if (qz .gt. rzero) den = t3
	1263	! test for zero denominator
	1264	if (cdabs(den) .eq. rzero) then
	1265	call trans_rt()
	1266	call deflated_work(errcode)
	1267	if (errcode == 1) then
	1268	exit icloop
	1269	end if
	1270	z2 = fprt
	1271	y2 = frt
	1272	cycle iloop
	1273	end if
	1274
	1275
	1276	d = -xx/den
	1277	d = d+d
	1278	h = d*h
	1279	rt = rt + h
	1280	! check convergence of the first kind
	1281	if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then
	1282	if (ic .ne. 0) then
	1283	exit icloop
	1284	end if
	1285	ic = 1
	1286	z0 = y1
	1287	z1 = y2
	1288	z2 = f(rt)
	1289	xl = rt
	1290	ick = ick+1
	1291	if (ick .le. ickmax) then
	1292	cycle iloop
	1293	end if
	1294	! warning error, itmax = maximum
	1295	jk = itmax + jk
	1296	ier = 33
	1297	end if
	1298	if (ic .ne. 0) then
	1299	cycle icloop
	1300	end if
	1301	call deflated_work(errcode)
	1302	if (errcode == 1) then
	1303	exit icloop
	1304	end if
	1305
	1306	do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)
	1307	! take remedial action to induce
	1308	! convergence
	1309	d = d*p5
	1310	h = h*p5
	1311	rt = rt-h
	1312	call deflated_work(errcode)
	1313	if (errcode == 1) then
	1314	exit icloop
	1315	end if
	1316	end do
	1317	z0 = z1
	1318	z1 = z2
	1319	z2 = fprt
	1320	y0 = y1
	1321	y1 = y2
	1322	y2 = frt
	1323	end do iloop
	1324	end do icloop
	1325	x(l) = rt
	1326	infer(l) = jk
	1327	l = l+1
	1328	end do rloop
	1329
	1330	contains
	1331	subroutine trans_rt()
	1332	tem = rten*eps1
	1333	if (cdabs(rt) .gt. ax) tem = tem*rt
	1334	rt = rt+tem
	1335	d = (h+tem)*d/h
	1336	h = h+tem
	1337	end subroutine trans_rt
	1338
	1339	subroutine deflated_work(errcode)
	1340	! errcode=0 => no errors
	1341	! errcode=1 => jk>itmax or convergence of second kind achieved
	1342	integer :: errcode,flag
	1343
	1344	flag=1
	1345	loop1: do while(flag==1)
	1346	errcode=0
	1347	jk = jk+1
	1348	if (jk .gt. itmax) then
	1349	ier=33
	1350	errcode=1
	1351	return
	1352	end if
	1353	frt = f(rt)
	1354	fprt = frt
	1355	if (l /= 1) then
	1356	lm1 = l-1
	1357	do i=1,lm1
	1358	tem = rt - x(i)
	1359	if (cdabs(tem) .eq. rzero) then
	1360	!if (ic .ne. 0) go to 15 !! ?? possible?
	1361	call trans_rt()
	1362	cycle loop1
	1363	end if
	1364	fprt = fprt/tem
	1365	end do
	1366	end if
	1367	flag=0
	1368	end do loop1
	1369
	1370	if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then
	1371	errcode=1
	1372	return
	1373	end if
	1374
	1375	end subroutine deflated_work
	1376
	1377	end subroutine
	1378
	1379
	1380	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1381	!
	1382	! Integration
	1383	!
	1384	! Only double precision coded atm
	1385	!
	1386	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1387
	1388
	1389	subroutine fvn_d_gauss_legendre(n,qx,qw)
	1390	!
	1391	! This routine compute the n Gauss Legendre abscissas and weights
	1392	! Adapted from Numerical Recipes routine gauleg
	1393	!
	1394	! n (in) : number of points
	1395	! qx(out) : abscissas
	1396	! qw(out) : weights
	1397	!
	1398	implicit none
	1399	double precision,parameter :: pi=3.141592653589793d0
	1400	integer, intent(in) :: n
	1401	double precision, intent(out) :: qx(n),qw(n)
	1402
	1403	integer :: m,i,j
	1404	double precision :: z,z1,p1,p2,p3,pp
	1405	m=(n+1)/2
	1406	do i=1,m
	1407	z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))
	1408	iloop: do
	1409	p1=1.d0
	1410	p2=0.d0
	1411	do j=1,n
	1412	p3=p2
	1413	p2=p1
	1414	p1=((2.d0dble(j)-1.d0)zp2-(dble(j)-1.d0)p3)/dble(j)
	1415	end do
	1416	pp=dble(n)(zp1-p2)/(z*z-1.d0)
	1417	z1=z
	1418	z=z1-p1/pp
	1419	if (dabs(z-z1)<=epsilon(z)) then
	1420	exit iloop
	1421	end if
	1422	end do iloop
	1423	qx(i)=-z
	1424	qx(n+1-i)=z
	1425	qw(i)=2.d0/((1.d0-zz)pp*pp)
	1426	qw(n+1-i)=qw(i)
	1427	end do
	1428	end subroutine
	1429
	1430
	1431
	1432	subroutine fvn_d_gl_integ(f,a,b,n,res)
	1433	!
	1434	! This is a simple non adaptative integration routine
	1435	! using n gauss legendre abscissas and weights
	1436	!
	1437	! f(in) : the function to integrate
	1438	! a(in) : lower bound
	1439	! b(in) : higher bound
	1440	! n(in) : number of gauss legendre pairs
	1441	! res(out): the evaluation of the integral
	1442	!
	1443	double precision,external :: f
	1444	double precision, intent(in) :: a,b
	1445	integer, intent(in):: n
	1446	double precision, intent(out) :: res
	1447
	1448	double precision, allocatable :: qx(:),qw(:)
	1449	double precision :: xm,xr
	1450	integer :: i
	1451
	1452	! First compute n gauss legendre abs and weight
	1453	allocate(qx(n))
	1454	allocate(qw(n))
	1455	call fvn_d_gauss_legendre(n,qx,qw)
	1456
	1457	xm=0.5d0*(b+a)
	1458	xr=0.5d0*(b-a)
	1459
	1460	res=0.d0
	1461
	1462	do i=1,n
	1463	res=res+qw(i)f(xm+xrqx(i))
	1464	end do
	1465
	1466	res=xr*res
	1467
	1468	deallocate(qw)
	1469	deallocate(qx)
	1470
	1471	end subroutine
	1472
	1473	!!!!!!!!!!!!!!!!!!!!!!!!
	1474	!
	1475	! Simple and double adaptative Gauss Kronrod integration based on
	1476	! a modified version of quadpack ( http://www.netlib.org/quadpack
	1477	!
	1478	! Common parameters :
	1479	!
	1480	! key (in)
	1481	! epsabs
	1482	! epsrel
	1483	!
	1484	!
	1485	!!!!!!!!!!!!!!!!!!!!!!!!
	1486
	1487	subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
	1488	!
	1489	! Evaluate the integral of function f(x) between a and b
	1490	!
	1491	! f(in) : the function
	1492	! a(in) : lower bound
	1493	! b(in) : higher bound
	1494	! epsabs(in) : desired absolute error
	1495	! epsrel(in) : desired relative error
	1496	! key(in) : gauss kronrod rule
	1497	! 1: 7 - 15 points
	1498	! 2: 10 - 21 points
	1499	! 3: 15 - 31 points
	1500	! 4: 20 - 41 points
	1501	! 5: 25 - 51 points
	1502	! 6: 30 - 61 points
	1503	!
	1504	! limit(in) : maximum number of subintervals in the partition of the
	1505	! given integration interval (a,b). A value of 500 will give the same
	1506	! behaviour as the imsl routine dqdag
	1507	!
	1508	! res(out) : estimated integral value
	1509	! abserr(out) : estimated absolute error
	1510	! ier(out) : error flag from quadpack routines
	1511	! 0 : no error
	1512	! 1 : maximum number of subdivisions allowed
	1513	! has been achieved. one can allow more
	1514	! subdivisions by increasing the value of
	1515	! limit (and taking the according dimension
	1516	! adjustments into account). however, if
	1517	! this yield no improvement it is advised
	1518	! to analyze the integrand in order to
	1519	! determine the integration difficulaties.
	1520	! if the position of a local difficulty can
	1521	! be determined (i.e.singularity,
	1522	! discontinuity within the interval) one
	1523	! will probably gain from splitting up the
	1524	! interval at this point and calling the
	1525	! integrator on the subranges. if possible,
	1526	! an appropriate special-purpose integrator
	1527	! should be used which is designed for
	1528	! handling the type of difficulty involved.
	1529	! 2 : the occurrence of roundoff error is
	1530	! detected, which prevents the requested
	1531	! tolerance from being achieved.
	1532	! 3 : extremely bad integrand behaviour occurs
	1533	! at some points of the integration
	1534	! interval.
	1535	! 6 : the input is invalid, because
	1536	! (epsabs.le.0 and
	1537	! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	1538	! or limit.lt.1 or lenw.lt.limit*4.
	1539	! result, abserr, neval, last are set
	1540	! to zero.
	1541	! except when lenw is invalid, iwork(1),
	1542	! work(limit2+1) and work(limit3+1) are
	1543	! set to zero, work(1) is set to a and
	1544	! work(limit+1) to b.
	1545
	1546	implicit none
	1547	double precision, external :: f
	1548	double precision, intent(in) :: a,b,epsabs,epsrel
	1549	integer, intent(in) :: key
	1550	integer, intent(in) :: limit
	1551	double precision, intent(out) :: res,abserr
	1552	integer, intent(out) :: ier
	1553
	1554	double precision, allocatable :: work(:)
	1555	integer, allocatable :: iwork(:)
	1556	integer :: lenw,neval,last
	1557
	1558	! imsl value for limit is 500
	1559	lenw=limit*4
	1560
	1561	allocate(iwork(limit))
	1562	allocate(work(lenw))
	1563
	1564	call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
	1565
	1566	deallocate(work)
	1567	deallocate(iwork)
	1568
	1569	end subroutine
	1570
	1571

testing/fvnlib.f90

Diff comments View file @ d1b0528

File was created	1
	2	module fvn
	3	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	4	!
	5	! fvn : a f95 module replacement for some imsl routines
	6	! it uses lapack for linear algebra
	7	! it uses modified quadpack for integration
	8	!
	9	! William Daniau 2007
	10	! william.daniau@femto-st.fr
	11	!
	12	! Routines naming scheme :
	13	!
	14	! fvn_x_name
	15	! where x can be s : real
	16	! d : real double precision
	17	! c : complex
	18	! z : double complex
	19	!
	20	!
	21	! This piece of code is totally free! Do whatever you want with it. However
	22	! if you find it usefull it would be kind to give credits ;-) Nevertheless, you
	23	! may give credits to quadpack authors.
	24	!
	25	! Version 1.1
	26	!
	27	! TO DO LIST :
	28	! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm
	29	! eigenvalues are given with no particular order.
	30	! + Generic interface for fvn_x_name family -> fvn_name
	31	! + Make some parameters optional, status for example
	32	! + use f95 kinds "double complex" -> complex(kind=8)
	33	! + unify quadpack routines
	34	! + ...
	35	!
	36	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	37
	38	implicit none
	39	! All quadpack routines are private to the module
	40	private :: d1mach,dqag,dqag_2d_inner,dqag_2d_outer,dqage,dqage_2d_inner, &
	41	dqage_2d_outer,dqk15,dqk15_2d_inner,dqk15_2d_outer,dqk21,dqk21_2d_inner,dqk21_2d_outer, &
	42	dqk31,dqk31_2d_inner,dqk31_2d_outer,dqk41,dqk41_2d_inner,dqk41_2d_outer, &
	43	dqk51,dqk51_2d_inner,dqk51_2d_outer,dqk61,dqk61_2d_inner,dqk61_2d_outer,dqpsrt
	44
	45
	46	contains
	47
	48	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	49	!
	50	! Matrix inversion subroutines
	51	!
	52	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	53	subroutine fvn_s_matinv(d,a,inva,status)
	54	!
	55	! Matrix inversion of a real matrix using BLAS and LAPACK
	56	!
	57	! d (in) : matrix rank
	58	! a (in) : input matrix
	59	! inva (out) : inversed matrix
	60	! status (ou) : =0 if something failed
	61	!
	62	integer, intent(in) :: d
	63	real, intent(in) :: a(d,d)
	64	real, intent(out) :: inva(d,d)
	65	integer, intent(out) :: status
	66
	67	integer, allocatable :: ipiv(:)
	68	real, allocatable :: work(:)
	69	real twork(1)
	70	integer :: info
	71	integer :: lwork
	72
	73	status=1
	74
	75	allocate(ipiv(d))
	76	! copy a into inva using BLAS
	77	!call scopy(d*d,a,1,inva,1)
	78	inva(:,:)=a(:,:)
	79	! LU factorization using LAPACK
	80	call sgetrf(d,d,inva,d,ipiv,info)
	81	! if info is not equal to 0, something went wrong we exit setting status to 0
	82	if (info /= 0) then
	83	status=0
	84	deallocate(ipiv)
	85	return
	86	end if
	87	! we use the query fonction of xxxtri to obtain the optimal workspace size
	88	call sgetri(d,inva,d,ipiv,twork,-1,info)
	89	lwork=int(twork(1))
	90	allocate(work(lwork))
	91	! Matrix inversion using LAPACK
	92	call sgetri(d,inva,d,ipiv,work,lwork,info)
	93	! again if info is not equal to 0, we exit setting status to 0
	94	if (info /= 0) then
	95	status=0
	96	end if
	97	deallocate(work)
	98	deallocate(ipiv)
	99	end subroutine
	100
	101	subroutine fvn_d_matinv(d,a,inva,status)
	102	!
	103	! Matrix inversion of a double precision matrix using BLAS and LAPACK
	104	!
	105	! d (in) : matrix rank
	106	! a (in) : input matrix
	107	! inva (out) : inversed matrix
	108	! status (ou) : =0 if something failed
	109	!
	110	integer, intent(in) :: d
	111	double precision, intent(in) :: a(d,d)
	112	double precision, intent(out) :: inva(d,d)
	113	integer, intent(out) :: status
	114
	115	integer, allocatable :: ipiv(:)
	116	double precision, allocatable :: work(:)
	117	double precision :: twork(1)
	118	integer :: info
	119	integer :: lwork
	120
	121	status=1
	122
	123	allocate(ipiv(d))
	124	! copy a into inva using BLAS
	125	!call dcopy(d*d,a,1,inva,1)
	126	inva(:,:)=a(:,:)
	127	! LU factorization using LAPACK
	128	call dgetrf(d,d,inva,d,ipiv,info)
	129	! if info is not equal to 0, something went wrong we exit setting status to 0
	130	if (info /= 0) then
	131	status=0
	132	deallocate(ipiv)
	133	return
	134	end if
	135	! we use the query fonction of xxxtri to obtain the optimal workspace size
	136	call dgetri(d,inva,d,ipiv,twork,-1,info)
	137	lwork=int(twork(1))
	138	allocate(work(lwork))
	139	! Matrix inversion using LAPACK
	140	call dgetri(d,inva,d,ipiv,work,lwork,info)
	141	! again if info is not equal to 0, we exit setting status to 0
	142	if (info /= 0) then
	143	status=0
	144	end if
	145	deallocate(work)
	146	deallocate(ipiv)
	147	end subroutine
	148
	149	subroutine fvn_c_matinv(d,a,inva,status)
	150	!
	151	! Matrix inversion of a complex matrix using BLAS and LAPACK
	152	!
	153	! d (in) : matrix rank
	154	! a (in) : input matrix
	155	! inva (out) : inversed matrix
	156	! status (ou) : =0 if something failed
	157	!
	158	integer, intent(in) :: d
	159	complex, intent(in) :: a(d,d)
	160	complex, intent(out) :: inva(d,d)
	161	integer, intent(out) :: status
	162
	163	integer, allocatable :: ipiv(:)
	164	complex, allocatable :: work(:)
	165	complex :: twork(1)
	166	integer :: info
	167	integer :: lwork
	168
	169	status=1
	170
	171	allocate(ipiv(d))
	172	! copy a into inva using BLAS
	173	!call ccopy(d*d,a,1,inva,1)
	174	inva(:,:)=a(:,:)
	175
	176	! LU factorization using LAPACK
	177	call cgetrf(d,d,inva,d,ipiv,info)
	178	! if info is not equal to 0, something went wrong we exit setting status to 0
	179	if (info /= 0) then
	180	status=0
	181	deallocate(ipiv)
	182	return
	183	end if
	184	! we use the query fonction of xxxtri to obtain the optimal workspace size
	185	call cgetri(d,inva,d,ipiv,twork,-1,info)
	186	lwork=int(twork(1))
	187	allocate(work(lwork))
	188	! Matrix inversion using LAPACK
	189	call cgetri(d,inva,d,ipiv,work,lwork,info)
	190	! again if info is not equal to 0, we exit setting status to 0
	191	if (info /= 0) then
	192	status=0
	193	end if
	194	deallocate(work)
	195	deallocate(ipiv)
	196	end subroutine
	197
	198	subroutine fvn_z_matinv(d,a,inva,status)
	199	!
	200	! Matrix inversion of a double complex matrix using BLAS and LAPACK
	201	!
	202	! d (in) : matrix rank
	203	! a (in) : input matrix
	204	! inva (out) : inversed matrix
	205	! status (ou) : =0 if something failed
	206	!
	207	integer, intent(in) :: d
	208	double complex, intent(in) :: a(d,d)
	209	double complex, intent(out) :: inva(d,d)
	210	integer, intent(out) :: status
	211
	212	integer, allocatable :: ipiv(:)
	213	double complex, allocatable :: work(:)
	214	double complex :: twork(1)
	215	integer :: info
	216	integer :: lwork
	217
	218	status=1
	219
	220	allocate(ipiv(d))
	221	! copy a into inva using BLAS
	222	!call zcopy(d*d,a,1,inva,1)
	223	inva(:,:)=a(:,:)
	224
	225	! LU factorization using LAPACK
	226	call zgetrf(d,d,inva,d,ipiv,info)
	227	! if info is not equal to 0, something went wrong we exit setting status to 0
	228	if (info /= 0) then
	229	status=0
	230	deallocate(ipiv)
	231	return
	232	end if
	233	! we use the query fonction of xxxtri to obtain the optimal workspace size
	234	call zgetri(d,inva,d,ipiv,twork,-1,info)
	235	lwork=int(twork(1))
	236	allocate(work(lwork))
	237	! Matrix inversion using LAPACK
	238	call zgetri(d,inva,d,ipiv,work,lwork,info)
	239	! again if info is not equal to 0, we exit setting status to 0
	240	if (info /= 0) then
	241	status=0
	242	end if
	243	deallocate(work)
	244	deallocate(ipiv)
	245	end subroutine
	246
	247
	248	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	249	!
	250	! Determinants
	251	!
	252	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	253	function fvn_s_det(d,a,status)
	254	!
	255	! Evaluate the determinant of a square matrix using lapack LU factorization
	256	!
	257	! d (in) : matrix rank
	258	! a (in) : The Matrix
	259	! status (out) : =0 if LU factorization failed
	260	!
	261	integer, intent(in) :: d
	262	real, intent(in) :: a(d,d)
	263	integer, intent(out) :: status
	264	real :: fvn_s_det
	265
	266	real, allocatable :: wc_a(:,:)
	267	integer, allocatable :: ipiv(:)
	268	integer :: info,i
	269
	270	status=1
	271	allocate(wc_a(d,d))
	272	allocate(ipiv(d))
	273	wc_a(:,:)=a(:,:)
	274	call sgetrf(d,d,wc_a,d,ipiv,info)
	275	if (info/= 0) then
	276	status=0
	277	fvn_s_det=0.e0
	278	deallocate(ipiv)
	279	deallocate(wc_a)
	280	return
	281	end if
	282	fvn_s_det=1.e0
	283	do i=1,d
	284	if (ipiv(i)==i) then
	285	fvn_s_det=fvn_s_det*wc_a(i,i)
	286	else
	287	fvn_s_det=-fvn_s_det*wc_a(i,i)
	288	end if
	289	end do
	290	deallocate(ipiv)
	291	deallocate(wc_a)
	292
	293	end function
	294
	295	function fvn_d_det(d,a,status)
	296	!
	297	! Evaluate the determinant of a square matrix using lapack LU factorization
	298	!
	299	! d (in) : matrix rank
	300	! a (in) : The Matrix
	301	! status (out) : =0 if LU factorization failed
	302	!
	303	integer, intent(in) :: d
	304	double precision, intent(in) :: a(d,d)
	305	integer, intent(out) :: status
	306	double precision :: fvn_d_det
	307
	308	double precision, allocatable :: wc_a(:,:)
	309	integer, allocatable :: ipiv(:)
	310	integer :: info,i
	311
	312	status=1
	313	allocate(wc_a(d,d))
	314	allocate(ipiv(d))
	315	wc_a(:,:)=a(:,:)
	316	call dgetrf(d,d,wc_a,d,ipiv,info)
	317	if (info/= 0) then
	318	status=0
	319	fvn_d_det=0.d0
	320	deallocate(ipiv)
	321	deallocate(wc_a)
	322	return
	323	end if
	324	fvn_d_det=1.d0
	325	do i=1,d
	326	if (ipiv(i)==i) then
	327	fvn_d_det=fvn_d_det*wc_a(i,i)
	328	else
	329	fvn_d_det=-fvn_d_det*wc_a(i,i)
	330	end if
	331	end do
	332	deallocate(ipiv)
	333	deallocate(wc_a)
	334
	335	end function
	336
	337	function fvn_c_det(d,a,status) !
	338	! Evaluate the determinant of a square matrix using lapack LU factorization
	339	!
	340	! d (in) : matrix rank
	341	! a (in) : The Matrix
	342	! status (out) : =0 if LU factorization failed
	343	!
	344	integer, intent(in) :: d
	345	complex, intent(in) :: a(d,d)
	346	integer, intent(out) :: status
	347	complex :: fvn_c_det
	348
	349	complex, allocatable :: wc_a(:,:)
	350	integer, allocatable :: ipiv(:)
	351	integer :: info,i
	352
	353	status=1
	354	allocate(wc_a(d,d))
	355	allocate(ipiv(d))
	356	wc_a(:,:)=a(:,:)
	357	call cgetrf(d,d,wc_a,d,ipiv,info)
	358	if (info/= 0) then
	359	status=0
	360	fvn_c_det=(0.e0,0.e0)
	361	deallocate(ipiv)
	362	deallocate(wc_a)
	363	return
	364	end if
	365	fvn_c_det=(1.e0,0.e0)
	366	do i=1,d
	367	if (ipiv(i)==i) then
	368	fvn_c_det=fvn_c_det*wc_a(i,i)
	369	else
	370	fvn_c_det=-fvn_c_det*wc_a(i,i)
	371	end if
	372	end do
	373	deallocate(ipiv)
	374	deallocate(wc_a)
	375
	376	end function
	377
	378	function fvn_z_det(d,a,status)
	379	!
	380	! Evaluate the determinant of a square matrix using lapack LU factorization
	381	!
	382	! d (in) : matrix rank
	383	! a (in) : The Matrix
	384	! det (out) : determinant
	385	! status (out) : =0 if LU factorization failed
	386	!
	387	integer, intent(in) :: d
	388	double complex, intent(in) :: a(d,d)
	389	integer, intent(out) :: status
	390	double complex :: fvn_z_det
	391
	392	double complex, allocatable :: wc_a(:,:)
	393	integer, allocatable :: ipiv(:)
	394	integer :: info,i
	395
	396	status=1
	397	allocate(wc_a(d,d))
	398	allocate(ipiv(d))
	399	wc_a(:,:)=a(:,:)
	400	call zgetrf(d,d,wc_a,d,ipiv,info)
	401	if (info/= 0) then
	402	status=0
	403	fvn_z_det=(0.d0,0.d0)
	404	deallocate(ipiv)
	405	deallocate(wc_a)
	406	return
	407	end if
	408	fvn_z_det=(1.d0,0.d0)
	409	do i=1,d
	410	if (ipiv(i)==i) then
	411	fvn_z_det=fvn_z_det*wc_a(i,i)
	412	else
	413	fvn_z_det=-fvn_z_det*wc_a(i,i)
	414	end if
	415	end do
	416	deallocate(ipiv)
	417	deallocate(wc_a)
	418
	419	end function
	420
	421	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	422	!
	423	! Condition test
	424	!
	425	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	426	! 1-norm
	427	! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm
	428	! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond
	429	!
	430	subroutine fvn_s_matcon(d,a,rcond,status)
	431	! Matrix condition (reciprocal of condition number)
	432	!
	433	! d (in) : matrix rank
	434	! a (in) : The Matrix
	435	! rcond (out) : guess what
	436	! status (out) : =0 if something went wrong
	437	!
	438	integer, intent(in) :: d
	439	real, intent(in) :: a(d,d)
	440	real, intent(out) :: rcond
	441	integer, intent(out) :: status
	442
	443	real, allocatable :: work(:)
	444	integer, allocatable :: iwork(:)
	445	real :: anorm
	446	real, allocatable :: wc_a(:,:) ! working copy of a
	447	integer :: info
	448	integer, allocatable :: ipiv(:)
	449
	450	real, external :: slange
	451
	452
	453	status=1
	454
	455	anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
	456
	457	allocate(wc_a(d,d))
	458	!call scopy(d*d,a,1,wc_a,1)
	459	wc_a(:,:)=a(:,:)
	460
	461	allocate(ipiv(d))
	462	call sgetrf(d,d,wc_a,d,ipiv,info)
	463	if (info /= 0) then
	464	status=0
	465	deallocate(ipiv)
	466	deallocate(wc_a)
	467	return
	468	end if
	469	allocate(work(4*d))
	470	allocate(iwork(d))
	471	call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
	472	if (info /= 0) then
	473	status=0
	474	end if
	475	deallocate(iwork)
	476	deallocate(work)
	477	deallocate(ipiv)
	478	deallocate(wc_a)
	479
	480	end subroutine
	481
	482	subroutine fvn_d_matcon(d,a,rcond,status)
	483	! Matrix condition (reciprocal of condition number)
	484	!
	485	! d (in) : matrix rank
	486	! a (in) : The Matrix
	487	! rcond (out) : guess what
	488	! status (out) : =0 if something went wrong
	489	!
	490	integer, intent(in) :: d
	491	double precision, intent(in) :: a(d,d)
	492	double precision, intent(out) :: rcond
	493	integer, intent(out) :: status
	494
	495	double precision, allocatable :: work(:)
	496	integer, allocatable :: iwork(:)
	497	double precision :: anorm
	498	double precision, allocatable :: wc_a(:,:) ! working copy of a
	499	integer :: info
	500	integer, allocatable :: ipiv(:)
	501
	502	double precision, external :: dlange
	503
	504
	505	status=1
	506
	507	anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
	508
	509	allocate(wc_a(d,d))
	510	!call dcopy(d*d,a,1,wc_a,1)
	511	wc_a(:,:)=a(:,:)
	512
	513	allocate(ipiv(d))
	514	call dgetrf(d,d,wc_a,d,ipiv,info)
	515	if (info /= 0) then
	516	status=0
	517	deallocate(ipiv)
	518	deallocate(wc_a)
	519	return
	520	end if
	521
	522	allocate(work(4*d))
	523	allocate(iwork(d))
	524	call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
	525	if (info /= 0) then
	526	status=0
	527	end if
	528	deallocate(iwork)
	529	deallocate(work)
	530	deallocate(ipiv)
	531	deallocate(wc_a)
	532
	533	end subroutine
	534
	535	subroutine fvn_c_matcon(d,a,rcond,status)
	536	! Matrix condition (reciprocal of condition number)
	537	!
	538	! d (in) : matrix rank
	539	! a (in) : The Matrix
	540	! rcond (out) : guess what
	541	! status (out) : =0 if something went wrong
	542	!
	543	integer, intent(in) :: d
	544	complex, intent(in) :: a(d,d)
	545	real, intent(out) :: rcond
	546	integer, intent(out) :: status
	547
	548	real, allocatable :: rwork(:)
	549	complex, allocatable :: work(:)
	550	integer, allocatable :: iwork(:)
	551	real :: anorm
	552	complex, allocatable :: wc_a(:,:) ! working copy of a
	553	integer :: info
	554	integer, allocatable :: ipiv(:)
	555
	556	real, external :: clange
	557
	558
	559	status=1
	560
	561	anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
	562
	563	allocate(wc_a(d,d))
	564	!call ccopy(d*d,a,1,wc_a,1)
	565	wc_a(:,:)=a(:,:)
	566
	567	allocate(ipiv(d))
	568	call cgetrf(d,d,wc_a,d,ipiv,info)
	569	if (info /= 0) then
	570	status=0
	571	deallocate(ipiv)
	572	deallocate(wc_a)
	573	return
	574	end if
	575	allocate(work(2*d))
	576	allocate(rwork(2*d))
	577	call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
	578	if (info /= 0) then
	579	status=0
	580	end if
	581	deallocate(rwork)
	582	deallocate(work)
	583	deallocate(ipiv)
	584	deallocate(wc_a)
	585	end subroutine
	586
	587	subroutine fvn_z_matcon(d,a,rcond,status)
	588	! Matrix condition (reciprocal of condition number)
	589	!
	590	! d (in) : matrix rank
	591	! a (in) : The Matrix
	592	! rcond (out) : guess what
	593	! status (out) : =0 if something went wrong
	594	!
	595	integer, intent(in) :: d
	596	double complex, intent(in) :: a(d,d)
	597	double precision, intent(out) :: rcond
	598	integer, intent(out) :: status
	599
	600	double complex, allocatable :: work(:)
	601	double precision, allocatable :: rwork(:)
	602	double precision :: anorm
	603	double complex, allocatable :: wc_a(:,:) ! working copy of a
	604	integer :: info
	605	integer, allocatable :: ipiv(:)
	606
	607	double precision, external :: zlange
	608
	609
	610	status=1
	611
	612	anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
	613
	614	allocate(wc_a(d,d))
	615	!call zcopy(d*d,a,1,wc_a,1)
	616	wc_a(:,:)=a(:,:)
	617
	618	allocate(ipiv(d))
	619	call zgetrf(d,d,wc_a,d,ipiv,info)
	620	if (info /= 0) then
	621	status=0
	622	deallocate(ipiv)
	623	deallocate(wc_a)
	624	return
	625	end if
	626
	627	allocate(work(2*d))
	628	allocate(rwork(2*d))
	629	call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
	630	if (info /= 0) then
	631	status=0
	632	end if
	633	deallocate(rwork)
	634	deallocate(work)
	635	deallocate(ipiv)
	636	deallocate(wc_a)
	637	end subroutine
	638
	639	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	640	!
	641	! Valeurs propres/ Vecteurs propre
	642	!
	643	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	644
	645	subroutine fvn_s_matev(d,a,evala,eveca,status)
	646	!
	647	! integer d (in) : matrice rank
	648	! real a(d,d) (in) : The Matrix
	649	! complex evala(d) (out) : eigenvalues
	650	! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
	651	! integer (out) : status =0 if something went wrong
	652	!
	653	! interfacing Lapack routine SGEEV
	654
	655	integer, intent(in) :: d
	656	real, intent(in) :: a(d,d)
	657	complex, intent(out) :: evala(d)
	658	complex, intent(out) :: eveca(d,d)
	659	integer, intent(out) :: status
	660
	661	real, allocatable :: wc_a(:,:) ! a working copy of a
	662	integer :: info
	663	integer :: lwork
	664	real, allocatable :: wr(:),wi(:)
	665	real :: vl ! unused but necessary for the call
	666	real, allocatable :: vr(:,:)
	667	real, allocatable :: work(:)
	668	real :: twork(1)
	669	integer i
	670	integer j
	671
	672	! making a working copy of a
	673	allocate(wc_a(d,d))
	674	!call scopy(d*d,a,1,wc_a,1)
	675	wc_a(:,:)=a(:,:)
	676
	677	allocate(wr(d))
	678	allocate(wi(d))
	679	allocate(vr(d,d))
	680	! query optimal work size
	681	call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
	682	lwork=int(twork(1))
	683	allocate(work(lwork))
	684	call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
	685
	686	if (info /= 0) then
	687	status=0
	688	deallocate(work)
	689	deallocate(vr)
	690	deallocate(wi)
	691	deallocate(wr)
	692	deallocate(wc_a)
	693	return
	694	end if
	695
	696	! now fill in the results
	697	i=1
	698	do while(i<=d)
	699	evala(i)=cmplx(wr(i),wi(i))
	700	if (wi(i) == 0.) then ! eigenvalue is real
	701	eveca(:,i)=cmplx(vr(:,i),0.)
	702	else ! eigenvalue is complex
	703	evala(i+1)=cmplx(wr(i+1),wi(i+1))
	704	eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))
	705	eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))
	706	i=i+1
	707	end if
	708	i=i+1
	709	enddo
	710	deallocate(work)
	711	deallocate(vr)
	712	deallocate(wi)
	713	deallocate(wr)
	714	deallocate(wc_a)
	715
	716	end subroutine
	717
	718	subroutine fvn_d_matev(d,a,evala,eveca,status)
	719	!
	720	! integer d (in) : matrice rank
	721	! double precision a(d,d) (in) : The Matrix
	722	! double complex evala(d) (out) : eigenvalues
	723	! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
	724	! integer (out) : status =0 if something went wrong
	725	!
	726	! interfacing Lapack routine DGEEV
	727	integer, intent(in) :: d
	728	double precision, intent(in) :: a(d,d)
	729	double complex, intent(out) :: evala(d)
	730	double complex, intent(out) :: eveca(d,d)
	731	integer, intent(out) :: status
	732
	733	double precision, allocatable :: wc_a(:,:) ! a working copy of a
	734	integer :: info
	735	integer :: lwork
	736	double precision, allocatable :: wr(:),wi(:)
	737	double precision :: vl ! unused but necessary for the call
	738	double precision, allocatable :: vr(:,:)
	739	double precision, allocatable :: work(:)
	740	double precision :: twork(1)
	741	integer i
	742	integer j
	743
	744	! making a working copy of a
	745	allocate(wc_a(d,d))
	746	!call dcopy(d*d,a,1,wc_a,1)
	747	wc_a(:,:)=a(:,:)
	748
	749	allocate(wr(d))
	750	allocate(wi(d))
	751	allocate(vr(d,d))
	752	! query optimal work size
	753	call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
	754	lwork=int(twork(1))
	755	allocate(work(lwork))
	756	call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
	757
	758	if (info /= 0) then
	759	status=0
	760	deallocate(work)
	761	deallocate(vr)
	762	deallocate(wi)
	763	deallocate(wr)
	764	deallocate(wc_a)
	765	return
	766	end if
	767
	768	! now fill in the results
	769	i=1
	770	do while(i<=d)
	771	evala(i)=dcmplx(wr(i),wi(i))
	772	if (wi(i) == 0.) then ! eigenvalue is real
	773	eveca(:,i)=dcmplx(vr(:,i),0.)
	774	else ! eigenvalue is complex
	775	evala(i+1)=dcmplx(wr(i+1),wi(i+1))
	776	eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))
	777	eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))
	778	i=i+1
	779	end if
	780	i=i+1
	781	enddo
	782
	783	deallocate(work)
	784	deallocate(vr)
	785	deallocate(wi)
	786	deallocate(wr)
	787	deallocate(wc_a)
	788
	789	end subroutine
	790
	791	subroutine fvn_c_matev(d,a,evala,eveca,status)
	792	!
	793	! integer d (in) : matrice rank
	794	! complex a(d,d) (in) : The Matrix
	795	! complex evala(d) (out) : eigenvalues
	796	! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
	797	! integer (out) : status =0 if something went wrong
	798	!
	799	! interfacing Lapack routine CGEEV
	800
	801	integer, intent(in) :: d
	802	complex, intent(in) :: a(d,d)
	803	complex, intent(out) :: evala(d)
	804	complex, intent(out) :: eveca(d,d)
	805	integer, intent(out) :: status
	806
	807	complex, allocatable :: wc_a(:,:) ! a working copy of a
	808	integer :: info
	809	integer :: lwork
	810	complex, allocatable :: work(:)
	811	complex :: twork(1)
	812	real, allocatable :: rwork(:)
	813	complex :: vl ! unused but necessary for the call
	814
	815	status=1
	816
	817	! making a working copy of a
	818	allocate(wc_a(d,d))
	819	!call ccopy(d*d,a,1,wc_a,1)
	820	wc_a(:,:)=a(:,:)
	821
	822
	823	! query optimal work size
	824	call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
	825	lwork=int(twork(1))
	826	allocate(work(lwork))
	827	allocate(rwork(2*d))
	828	call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
	829
	830	if (info /= 0) then
	831	status=0
	832	end if
	833	deallocate(rwork)
	834	deallocate(work)
	835	deallocate(wc_a)
	836
	837	end subroutine
	838
	839	subroutine fvn_z_matev(d,a,evala,eveca,status)
	840	!
	841	! integer d (in) : matrice rank
	842	! double complex a(d,d) (in) : The Matrix
	843	! double complex evala(d) (out) : eigenvalues
	844	! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
	845	! integer (out) : status =0 if something went wrong
	846	!
	847	! interfacing Lapack routine ZGEEV
	848
	849	integer, intent(in) :: d
	850	double complex, intent(in) :: a(d,d)
	851	double complex, intent(out) :: evala(d)
	852	double complex, intent(out) :: eveca(d,d)
	853	integer, intent(out) :: status
	854
	855	double complex, allocatable :: wc_a(:,:) ! a working copy of a
	856	integer :: info
	857	integer :: lwork
	858	double complex, allocatable :: work(:)
	859	double complex :: twork(1)
	860	double precision, allocatable :: rwork(:)
	861	double complex :: vl ! unused but necessary for the call
	862
	863	status=1
	864
	865	! making a working copy of a
	866	allocate(wc_a(d,d))
	867	!call zcopy(d*d,a,1,wc_a,1)
	868	wc_a(:,:)=a(:,:)
	869
	870	! query optimal work size
	871	call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
	872	lwork=int(twork(1))
	873	allocate(work(lwork))
	874	allocate(rwork(2*d))
	875	call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
	876
	877	if (info /= 0) then
	878	status=0
	879	end if
	880	deallocate(rwork)
	881	deallocate(work)
	882	deallocate(wc_a)
	883
	884	end subroutine
	885
	886	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	887	!
	888	! Akima spline interpolation and spline evaluation
	889	!
	890	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	891
	892	! Single precision
	893	subroutine fvn_s_akima(n,x,y,br,co)
	894	implicit none
	895	integer, intent(in) :: n
	896	real, intent(in) :: x(n)
	897	real, intent(in) :: y(n)
	898	real, intent(out) :: br(n)
	899	real, intent(out) :: co(4,n)
	900
	901	real, allocatable :: var(:),z(:)
	902	real :: wi_1,wi
	903	integer :: i
	904	real :: dx,a,b
	905
	906	! br is just a copy of x
	907	br(:)=x(:)
	908
	909	allocate(var(n))
	910	allocate(z(n))
	911	! evaluate the variations
	912	do i=1, n-1
	913	var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
	914	end do
	915	var(n+2)=2.e0*var(n+1)-var(n)
	916	var(n+3)=2.e0*var(n+2)-var(n+1)
	917	var(2)=2.e0*var(3)-var(4)
	918	var(1)=2.e0*var(2)-var(3)
	919
	920	do i = 1, n
	921	wi_1=abs(var(i+3)-var(i+2))
	922	wi=abs(var(i+1)-var(i))
	923	if ((wi_1+wi).eq.0.e0) then
	924	z(i)=(var(i+2)+var(i+1))/2.e0
	925	else
	926	z(i)=(wi_1var(i+1)+wivar(i+2))/(wi_1+wi)
	927	end if
	928	end do
	929
	930	do i=1, n-1
	931	dx=x(i+1)-x(i)
	932	a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd
	933	b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
	934	co(1,i)=y(i)
	935	co(2,i)=z(i)
	936	!co(3,i)=-(a-3.b)/dx*2 ! méthode wd
	937	!co(4,i)=(a-2.b)/dx*3 ! méthode wd
	938	co(3,i)=(3.e0var(i+2)-2.e0z(i)-z(i+1))/dx ! méthode JP Moreau
	939	co(4,i)=(z(i)+z(i+1)-2.e0var(i+2))/dx*2 !
	940	! les coefficients donnés par imsl sont co(3,i)2 et co(4,i)6
	941	! etrangement la fonction csval corrige et donne la bonne valeur ...
	942	end do
	943	co(1,n)=y(n)
	944	co(2,n)=z(n)
	945	co(3,n)=0.e0
	946	co(4,n)=0.e0
	947
	948	deallocate(z)
	949	deallocate(var)
	950
	951	end subroutine
	952
	953	! Double precision
	954	subroutine fvn_d_akima(n,x,y,br,co)
	955
	956	implicit none
	957	integer, intent(in) :: n
	958	double precision, intent(in) :: x(n)
	959	double precision, intent(in) :: y(n)
	960	double precision, intent(out) :: br(n)
	961	double precision, intent(out) :: co(4,n)
	962
	963	double precision, allocatable :: var(:),z(:)
	964	double precision :: wi_1,wi
	965	integer :: i
	966	double precision :: dx,a,b
	967
	968	! br is just a copy of x
	969	br(:)=x(:)
	970
	971	allocate(var(n))
	972	allocate(z(n))
	973	! evaluate the variations
	974	do i=1, n-1
	975	var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
	976	end do
	977	var(n+2)=2.d0*var(n+1)-var(n)
	978	var(n+3)=2.d0*var(n+2)-var(n+1)
	979	var(2)=2.d0*var(3)-var(4)
	980	var(1)=2.d0*var(2)-var(3)
	981
	982	do i = 1, n
	983	wi_1=dabs(var(i+3)-var(i+2))
	984	wi=dabs(var(i+1)-var(i))
	985	if ((wi_1+wi).eq.0.d0) then
	986	z(i)=(var(i+2)+var(i+1))/2.d0
	987	else
	988	z(i)=(wi_1var(i+1)+wivar(i+2))/(wi_1+wi)
	989	end if
	990	end do
	991
	992	do i=1, n-1
	993	dx=x(i+1)-x(i)
	994	a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd
	995	b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
	996	co(1,i)=y(i)
	997	co(2,i)=z(i)
	998	!co(3,i)=-(a-3.b)/dx*2 ! méthode wd
	999	!co(4,i)=(a-2.b)/dx*3 ! méthode wd
	1000	co(3,i)=(3.d0var(i+2)-2.d0z(i)-z(i+1))/dx ! méthode JP Moreau
	1001	co(4,i)=(z(i)+z(i+1)-2.d0var(i+2))/dx*2 !
	1002	! les coefficients donnés par imsl sont co(3,i)2 et co(4,i)6
	1003	! etrangement la fonction csval corrige et donne la bonne valeur ...
	1004	end do
	1005	co(1,n)=y(n)
	1006	co(2,n)=z(n)
	1007	co(3,n)=0.d0
	1008	co(4,n)=0.d0
	1009
	1010	deallocate(z)
	1011	deallocate(var)
	1012
	1013	end subroutine
	1014
	1015	!
	1016	! Single precision spline evaluation
	1017	!
	1018	function fvn_s_spline_eval(x,n,br,co)
	1019	implicit none
	1020	real, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
	1021	integer, intent(in) :: n ! number of intervals
	1022	real, intent(in) :: br(n+1) ! breakpoints
	1023	real, intent(in) :: co(4,n+1) ! spline coeeficients
	1024	real :: fvn_s_spline_eval
	1025
	1026	integer :: i
	1027	real :: dx
	1028
	1029	if (x<=br(1)) then
	1030	i=1
	1031	else if (x>=br(n+1)) then
	1032	i=n
	1033	else
	1034	i=1
	1035	do while(x>=br(i))
	1036	i=i+1
	1037	end do
	1038	i=i-1
	1039	end if
	1040	dx=x-br(i)
	1041	fvn_s_spline_eval=co(1,i)+co(2,i)dx+co(3,i)dx*2+co(4,i)dx**3
	1042
	1043	end function
	1044
	1045	! Double precision spline evaluation
	1046	function fvn_d_spline_eval(x,n,br,co)
	1047	implicit none
	1048	double precision, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
	1049	integer, intent(in) :: n ! number of intervals
	1050	double precision, intent(in) :: br(n+1) ! breakpoints
	1051	double precision, intent(in) :: co(4,n+1) ! spline coeeficients
	1052	double precision :: fvn_d_spline_eval
	1053
	1054	integer :: i
	1055	double precision :: dx
	1056
	1057
	1058	if (x<=br(1)) then
	1059	i=1
	1060	else if (x>=br(n+1)) then
	1061	i=n
	1062	else
	1063	i=1
	1064	do while(x>=br(i))
	1065	i=i+1
	1066	end do
	1067	i=i-1
	1068	end if
	1069
	1070	dx=x-br(i)
	1071	fvn_d_spline_eval=co(1,i)+co(2,i)dx+co(3,i)dx*2+co(4,i)dx**3
	1072
	1073	end function
	1074
	1075
	1076	!
	1077	! Muller
	1078	!
	1079	!
	1080	!
	1081	! William Daniau 2007
	1082	!
	1083	! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f
	1084	! http://plato.asu.edu/ftp/other_software/muller.f
	1085	!
	1086	! it can be used as a replacement for imsl routine dzanly with minor changes
	1087	!
	1088	!-----------------------------------------------------------------------
	1089	!
	1090	! purpose - zeros of an analytic complex function
	1091	! using the muller method with deflation
	1092	!
	1093	! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,
	1094	! infer,ier)
	1095	!
	1096	! arguments f - a complex function subprogram, f(z), written
	1097	! by the user specifying the equation whose
	1098	! roots are to be found. f must appear in
	1099	! an external statement in the calling pro-
	1100	! gram.
	1101	! eps - 1st stopping criterion. let fp(z)=f(z)/p
	1102	! where p = (z-z(1))(z-z(2)),,,*(z-z(k-1))
	1103	! and z(1),...,z(k-1) are previously found
	1104	! roots. if ((cdabs(f(z)).le.eps) .and.
	1105	! (cdabs(fp(z)).le.eps)), then z is accepted
	1106	! as a root. (input)
	1107	! eps1 - 2nd stopping criterion. a root is accepted
	1108	! if two successive approximations to a given
	1109	! root agree within eps1. (input)
	1110	! note. if either or both of the stopping
	1111	! criteria are fulfilled, the root is
	1112	! accepted.
	1113	! kn - the number of known roots which must be stored
	1114	! in x(1),...,x(kn), prior to entry to muller
	1115	! nguess - the number of initial guesses provided. these
	1116	! guesses must be stored in x(kn+1),...,
	1117	! x(kn+nguess). nguess must be set equal
	1118	! to zero if no guesses are provided. (input)
	1119	! n - the number of new roots to be found by
	1120	! muller (input)
	1121	! x - a complex vector of length kn+n. x(1),...,
	1122	! x(kn) on input must contain any known
	1123	! roots. x(kn+1),..., x(kn+n) on input may,
	1124	! on user option, contain initial guesses for
	1125	! the n new roots which are to be computed.
	1126	! if the user does not provide an initial
	1127	! guess, zero is used.
	1128	! on output, x(kn+1),...,x(kn+n) contain the
	1129	! approximate roots found by muller.
	1130	! itmax - the maximum allowable number of iterations
	1131	! per root (input)
	1132	! infer - an integer vector of length kn+n. on
	1133	! output infer(j) contains the number of
	1134	! iterations used in finding the j-th root
	1135	! when convergence was achieved. if
	1136	! convergence was not obtained in itmax
	1137	! iterations, infer(j) will be greater than
	1138	! itmax (output).
	1139	! ier - error parameter (output)
	1140	! warning error
	1141	! ier = 33 indicates failure to converge with-
	1142	! in itmax iterations for at least one of
	1143	! the (n) new roots.
	1144	!
	1145	!
	1146	! remarks muller always returns the last approximation for root j
	1147	! in x(j). if the convergence criterion is satisfied,
	1148	! then infer(j) is less than or equal to itmax. if the
	1149	! convergence criterion is not satisified, then infer(j)
	1150	! is set to either itmax+1 or itmax+k, with k greater
	1151	! than 1. infer(j) = itmax+1 indicates that muller did
	1152	! not obtain convergence in the allowed number of iter-
	1153	! ations. in this case, the user may wish to set itmax
	1154	! to a larger value. infer(j) = itmax+k means that con-
	1155	! vergence was obtained (on iteration k) for the defla-
	1156	! ted function
	1157	! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))
	1158	!
	1159	! but failed for f(z). in this case, better initial
	1160	! guesses might help or, it might be necessary to relax
	1161	! the convergence criterion.
	1162	!
	1163	!-----------------------------------------------------------------------
	1164	!
	1165	subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
	1166	implicit none
	1167	double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq
	1168	double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
	1169	tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
	1170	zero,p1,one,four,p5
	1171
	1172	double complex, external :: f
	1173	integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
	1174	knpng,jk,ick,nn,lm1,errcode
	1175	double complex :: x(kn+n)
	1176	integer :: infer(kn+n)
	1177
	1178
	1179	data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &
	1180	one/(1.0,0.0)/,four/(4.0,0.0)/, &
	1181	p5/(0.5,0.0)/, &
	1182	rzero/0.0/,rten/10.0/,rhun/100.0/, &
	1183	ax/0.1/,ickmax/3/,rp01/0.01/
	1184
	1185	ier = 0
	1186	if (n .lt. 1) then ! What the hell are doing here then ...
	1187	return
	1188	end if
	1189	!eps1 = rten **(-nsig)
	1190	eps1 = min(eps1,rp01)
	1191
	1192	knp1 = kn+1
	1193	knpn = kn+n
	1194	knpng = kn+nguess
	1195	do i=1,knpn
	1196	infer(i) = 0
	1197	if (i .gt. knpng) x(i) = zero
	1198	end do
	1199	l= knp1
	1200
	1201	ic=0
	1202	rloop: do while (l<=knpn) ! Main loop over new roots
	1203	jk = 0
	1204	ick = 0
	1205	xl = x(l)
	1206	icloop: do
	1207	ic = 0
	1208	h = ax
	1209	h = p1*h
	1210	if (cdabs(xl) .gt. ax) h = p1*xl
	1211	! first three points are
	1212	! xl+h, xl-h, xl
	1213	rt = xl+h
	1214	call deflated_work(errcode)
	1215	if (errcode == 1) then
	1216	exit icloop
	1217	end if
	1218
	1219	z0 = fprt
	1220	y0 = frt
	1221	x0 = rt
	1222	rt = xl-h
	1223	call deflated_work(errcode)
	1224	if (errcode == 1) then
	1225	exit icloop
	1226	end if
	1227
	1228	z1 = fprt
	1229	y1 = frt
	1230	h = xl-rt
	1231	d = h/(rt-x0)
	1232	rt = xl
	1233
	1234	call deflated_work(errcode)
	1235	if (errcode == 1) then
	1236	exit icloop
	1237	end if
	1238
	1239
	1240	z2 = fprt
	1241	y2 = frt
	1242	! begin main algorithm
	1243	iloop: do
	1244	dd = one + d
	1245	t1 = z0dd
	1246	t2 = z1dddd
	1247	xx = z2*dd
	1248	t3 = z2*d
	1249	bi = t1-t2+xx+t3
	1250	den = bibi-four(xxt1-t3(t2-xx))
	1251	! use denominator of maximum amplitude
	1252	t1 = cdsqrt(den)
	1253	qz = rhun*max(cdabs(bi),cdabs(t1))
	1254	t2 = bi + t1
	1255	tpq = cdabs(t2)+qz
	1256	if (tpq .eq. qz) t2 = zero
	1257	t3 = bi - t1
	1258	tpq = cdabs(t3) + qz
	1259	if (tpq .eq. qz) t3 = zero
	1260	den = t2
	1261	qz = cdabs(t3)-cdabs(t2)
	1262	if (qz .gt. rzero) den = t3
	1263	! test for zero denominator
	1264	if (cdabs(den) .eq. rzero) then
	1265	call trans_rt()
	1266	call deflated_work(errcode)
	1267	if (errcode == 1) then
	1268	exit icloop
	1269	end if
	1270	z2 = fprt
	1271	y2 = frt
	1272	cycle iloop
	1273	end if
	1274
	1275
	1276	d = -xx/den
	1277	d = d+d
	1278	h = d*h
	1279	rt = rt + h
	1280	! check convergence of the first kind
	1281	if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then
	1282	if (ic .ne. 0) then
	1283	exit icloop
	1284	end if
	1285	ic = 1
	1286	z0 = y1
	1287	z1 = y2
	1288	z2 = f(rt)
	1289	xl = rt
	1290	ick = ick+1
	1291	if (ick .le. ickmax) then
	1292	cycle iloop
	1293	end if
	1294	! warning error, itmax = maximum
	1295	jk = itmax + jk
	1296	ier = 33
	1297	end if
	1298	if (ic .ne. 0) then
	1299	cycle icloop
	1300	end if
	1301	call deflated_work(errcode)
	1302	if (errcode == 1) then
	1303	exit icloop
	1304	end if
	1305
	1306	do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)
	1307	! take remedial action to induce
	1308	! convergence
	1309	d = d*p5
	1310	h = h*p5
	1311	rt = rt-h
	1312	call deflated_work(errcode)
	1313	if (errcode == 1) then
	1314	exit icloop
	1315	end if
	1316	end do
	1317	z0 = z1
	1318	z1 = z2
	1319	z2 = fprt
	1320	y0 = y1
	1321	y1 = y2
	1322	y2 = frt
	1323	end do iloop
	1324	end do icloop
	1325	x(l) = rt
	1326	infer(l) = jk
	1327	l = l+1
	1328	end do rloop
	1329
	1330	contains
	1331	subroutine trans_rt()
	1332	tem = rten*eps1
	1333	if (cdabs(rt) .gt. ax) tem = tem*rt
	1334	rt = rt+tem
	1335	d = (h+tem)*d/h
	1336	h = h+tem
	1337	end subroutine trans_rt
	1338
	1339	subroutine deflated_work(errcode)
	1340	! errcode=0 => no errors
	1341	! errcode=1 => jk>itmax or convergence of second kind achieved
	1342	integer :: errcode,flag
	1343
	1344	flag=1
	1345	loop1: do while(flag==1)
	1346	errcode=0
	1347	jk = jk+1
	1348	if (jk .gt. itmax) then
	1349	ier=33
	1350	errcode=1
	1351	return
	1352	end if
	1353	frt = f(rt)
	1354	fprt = frt
	1355	if (l /= 1) then
	1356	lm1 = l-1
	1357	do i=1,lm1
	1358	tem = rt - x(i)
	1359	if (cdabs(tem) .eq. rzero) then
	1360	!if (ic .ne. 0) go to 15 !! ?? possible?
	1361	call trans_rt()
	1362	cycle loop1
	1363	end if
	1364	fprt = fprt/tem
	1365	end do
	1366	end if
	1367	flag=0
	1368	end do loop1
	1369
	1370	if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then
	1371	errcode=1
	1372	return
	1373	end if
	1374
	1375	end subroutine deflated_work
	1376
	1377	end subroutine
	1378
	1379
	1380	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1381	!
	1382	! Integration
	1383	!
	1384	! Only double precision coded atm
	1385	!
	1386	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1387
	1388
	1389	subroutine fvn_d_gauss_legendre(n,qx,qw)
	1390	!
	1391	! This routine compute the n Gauss Legendre abscissas and weights
	1392	! Adapted from Numerical Recipes routine gauleg
	1393	!
	1394	! n (in) : number of points
	1395	! qx(out) : abscissas
	1396	! qw(out) : weights
	1397	!
	1398	implicit none
	1399	double precision,parameter :: pi=3.141592653589793d0
	1400	integer, intent(in) :: n
	1401	double precision, intent(out) :: qx(n),qw(n)
	1402
	1403	integer :: m,i,j
	1404	double precision :: z,z1,p1,p2,p3,pp
	1405	m=(n+1)/2
	1406	do i=1,m
	1407	z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))
	1408	iloop: do
	1409	p1=1.d0
	1410	p2=0.d0
	1411	do j=1,n
	1412	p3=p2
	1413	p2=p1
	1414	p1=((2.d0dble(j)-1.d0)zp2-(dble(j)-1.d0)p3)/dble(j)
	1415	end do
	1416	pp=dble(n)(zp1-p2)/(z*z-1.d0)
	1417	z1=z
	1418	z=z1-p1/pp
	1419	if (dabs(z-z1)<=epsilon(z)) then
	1420	exit iloop
	1421	end if
	1422	end do iloop
	1423	qx(i)=-z
	1424	qx(n+1-i)=z
	1425	qw(i)=2.d0/((1.d0-zz)pp*pp)
	1426	qw(n+1-i)=qw(i)
	1427	end do
	1428	end subroutine
	1429
	1430
	1431
	1432	subroutine fvn_d_gl_integ(f,a,b,n,res)
	1433	!
	1434	! This is a simple non adaptative integration routine
	1435	! using n gauss legendre abscissas and weights
	1436	!
	1437	! f(in) : the function to integrate
	1438	! a(in) : lower bound
	1439	! b(in) : higher bound
	1440	! n(in) : number of gauss legendre pairs
	1441	! res(out): the evaluation of the integral
	1442	!
	1443	double precision,external :: f
	1444	double precision, intent(in) :: a,b
	1445	integer, intent(in):: n
	1446	double precision, intent(out) :: res
	1447
	1448	double precision, allocatable :: qx(:),qw(:)
	1449	double precision :: xm,xr
	1450	integer :: i
	1451
	1452	! First compute n gauss legendre abs and weight
	1453	allocate(qx(n))
	1454	allocate(qw(n))
	1455	call fvn_d_gauss_legendre(n,qx,qw)
	1456
	1457	xm=0.5d0*(b+a)
	1458	xr=0.5d0*(b-a)
	1459
	1460	res=0.d0
	1461
	1462	do i=1,n
	1463	res=res+qw(i)f(xm+xrqx(i))
	1464	end do
	1465
	1466	res=xr*res
	1467
	1468	deallocate(qw)
	1469	deallocate(qx)
	1470
	1471	end subroutine
	1472
	1473	!!!!!!!!!!!!!!!!!!!!!!!!
	1474	!
	1475	! Simple and double adaptative Gauss Kronrod integration based on
	1476	! a modified version of quadpack ( http://www.netlib.org/quadpack
	1477	!
	1478	! Common parameters :
	1479	!
	1480	! key (in)
	1481	! epsabs
	1482	! epsrel
	1483	!
	1484	!
	1485	!!!!!!!!!!!!!!!!!!!!!!!!
	1486
	1487	subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
	1488	!
	1489	! Evaluate the integral of function f(x) between a and b
	1490	!
	1491	! f(in) : the function
	1492	! a(in) : lower bound
	1493	! b(in) : higher bound
	1494	! epsabs(in) : desired absolute error
	1495	! epsrel(in) : desired relative error
	1496	! key(in) : gauss kronrod rule
	1497	! 1: 7 - 15 points
	1498	! 2: 10 - 21 points
	1499	! 3: 15 - 31 points
	1500	! 4: 20 - 41 points
	1501	! 5: 25 - 51 points
	1502	! 6: 30 - 61 points
	1503	!
	1504	! limit(in) : maximum number of subintervals in the partition of the
	1505	! given integration interval (a,b). A value of 500 will give the same
	1506	! behaviour as the imsl routine dqdag
	1507	!
	1508	! res(out) : estimated integral value
	1509	! abserr(out) : estimated absolute error
	1510	! ier(out) : error flag from quadpack routines
	1511	! 0 : no error
	1512	! 1 : maximum number of subdivisions allowed
	1513	! has been achieved. one can allow more
	1514	! subdivisions by increasing the value of
	1515	! limit (and taking the according dimension
	1516	! adjustments into account). however, if
	1517	! this yield no improvement it is advised
	1518	! to analyze the integrand in order to
	1519	! determine the integration difficulaties.
	1520	! if the position of a local difficulty can
	1521	! be determined (i.e.singularity,
	1522	! discontinuity within the interval) one
	1523	! will probably gain from splitting up the
	1524	! interval at this point and calling the
	1525	! integrator on the subranges. if possible,
	1526	! an appropriate special-purpose integrator
	1527	! should be used which is designed for
	1528	! handling the type of difficulty involved.
	1529	! 2 : the occurrence of roundoff error is
	1530	! detected, which prevents the requested
	1531	! tolerance from being achieved.
	1532	! 3 : extremely bad integrand behaviour occurs
	1533	! at some points of the integration
	1534	! interval.
	1535	! 6 : the input is invalid, because
	1536	! (epsabs.le.0 and
	1537	! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	1538	! or limit.lt.1 or lenw.lt.limit*4.
	1539	! result, abserr, neval, last are set
	1540	! to zero.
	1541	! except when lenw is invalid, iwork(1),
	1542	! work(limit2+1) and work(limit3+1) are
	1543	! set to zero, work(1) is set to a and
	1544	! work(limit+1) to b.
	1545
	1546	implicit none
	1547	double precision, external :: f
	1548	double precision, intent(in) :: a,b,epsabs,epsrel
	1549	integer, intent(in) :: key
	1550	integer, intent(in) :: limit
	1551	double precision, intent(out) :: res,abserr
	1552	integer, intent(out) :: ier
	1553
	1554	double precision, allocatable :: work(:)
	1555	integer, allocatable :: iwork(:)
	1556	integer :: lenw,neval,last
	1557
	1558	! imsl value for limit is 500
	1559	lenw=limit*4
	1560
	1561	allocate(iwork(limit))
	1562	allocate(work(lenw))
	1563
	1564	call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
	1565
	1566	deallocate(work)
	1567	deallocate(iwork)
	1568
	1569	end subroutine
	1570
	1571