lxsd / fvn | GitLab

Browse Code »

Commit f85d3b317a5f308dca1d53a0eb67db8d54d41c36

Authored by daniau 2007-09-14 11:33:16 +0200

1 parent 856145b45a

Exists in master and in 3 other branches

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

Showing 3 changed files with 257 additions and 0 deletions Inline Diff

doc/fvn.pdf
fvn_sparse/umfpack_wrapper.c
fvnlib.f90

doc/fvn.pdf

Diff comments View file @ f85d3b3

No preview for this file type

fvn_sparse/umfpack_wrapper.c

Diff comments View file @ f85d3b3

#include "umfpack.h"	1	1	#include "umfpack.h"
#include <ctype.h>	2	2	#include <ctype.h>
#include <stdio.h>	3	3	#include <stdio.h>
#ifdef NULL	4	4	#ifdef NULL
#undef NULL	5	5	#undef NULL
#endif	6	6	#endif
#define NULL 0	7	7	#define NULL 0
#define LEN 200	8	8	#define LEN 200
	9	9
/* for umfpack 4.1 */	10	10	/* for umfpack 4.1 */
/* #define UF_long long */	11	11	/* #define UF_long long */
/*	12	12	/*
	13	13
set of routines used by fvn_sparse	14	14	set of routines used by fvn_sparse
	15	15
complex(8) and integer(8)	16	16	complex(8) and integer(8)
	17	17
*/	18	18	*/
	19	19
/* defaults */	20	20	/* defaults */
	21	21
void umfpack_zl_defaults_ (double Control [UMFPACK_CONTROL])	22	22	void umfpack_zl_defaults_ (double Control [UMFPACK_CONTROL])
{	23	23	{
umfpack_zl_defaults (Control) ;	24	24	umfpack_zl_defaults (Control) ;
}	25	25	}
	26	26
/* free the Numeric object */	27	27	/* free the Numeric object */
void umfpack_zl_free_numeric_(void **Numeric)	28	28	void umfpack_zl_free_numeric_(void **Numeric)
{	29	29	{
umfpack_zl_free_numeric (Numeric) ;	30	30	umfpack_zl_free_numeric (Numeric) ;
}	31	31	}
	32	32
/* free th Symbolic object */	33	33	/* free th Symbolic object */
void umfpack_zl_free_symbolic_(void **Symbolic)	34	34	void umfpack_zl_free_symbolic_(void **Symbolic)
{	35	35	{
umfpack_zl_free_symbolic (Symbolic) ;	36	36	umfpack_zl_free_symbolic (Symbolic) ;
}	37	37	}
	38	38
/* numeric factorization */	39	39	/* numeric factorization */
void umfpack_zl_numeric_ (UF_long Ap [ ], UF_long Ai [ ], double Ax [ ], double Az [ ],	40	40	void umfpack_zl_numeric_ (UF_long Ap [ ], UF_long Ai [ ], double Ax [ ], double Az [ ],
void Symbolic, void Numeric,	41	41	void Symbolic, void Numeric,
double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])	42	42	double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
{	43	43	{
(void) umfpack_zl_numeric (Ap, Ai, Ax, Az, *Symbolic, Numeric, Control, Info);	44	44	(void) umfpack_zl_numeric (Ap, Ai, Ax, Az, *Symbolic, Numeric, Control, Info);
}	45	45	}
	46	46
/* pre-ordering and symbolic factorization */	47	47	/* pre-ordering and symbolic factorization */
void umfpack_zl_symbolic_ (UF_long m, UF_long n, UF_long Ap [ ], UF_long Ai [ ],	48	48	void umfpack_zl_symbolic_ (UF_long m, UF_long n, UF_long Ap [ ], UF_long Ai [ ],
double Ax [ ], double Az [ ], void **Symbolic,	49	49	double Ax [ ], double Az [ ], void **Symbolic,
double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])	50	50	double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
{	51	51	{
(void) umfpack_zl_symbolic (m, n, Ap, Ai, Ax, Az, Symbolic, Control, Info) ;	52	52	(void) umfpack_zl_symbolic (m, n, Ap, Ai, Ax, Az, Symbolic, Control, Info) ;
}	53	53	}
	54	54
/* solve a linear system */	55	55	/* solve a linear system */
	56	56
void umfpack_zl_solve_ (UF_long *sys, UF_long Ap [ ], UF_long Ai [ ], double Ax [ ], double Az [ ],	57	57	void umfpack_zl_solve_ (UF_long *sys, UF_long Ap [ ], UF_long Ai [ ], double Ax [ ], double Az [ ],
double x [ ], double xz [ ], double b [ ], double bz [ ], void **Numeric,	58	58	double x [ ], double xz [ ], double b [ ], double bz [ ], void **Numeric,
double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])	59	59	double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
{	60	60	{
(void) umfpack_zl_solve (*sys, Ap, Ai, Ax, Az, x, xz, b, bz,	61	61	(void) umfpack_zl_solve (*sys, Ap, Ai, Ax, Az, x, xz, b, bz,
*Numeric, Control, Info) ;	62	62	*Numeric, Control, Info) ;
}	63	63	}
	64	64
/* triplet 2 col */	65	65	/* triplet 2 col */
	66	66
void umfpack_zl_triplet_to_col_ (UF_long m, UF_long n, UF_long *nz, UF_long Ti [ ], UF_long Tj [ ],	67	67	void umfpack_zl_triplet_to_col_ (UF_long m, UF_long n, UF_long *nz, UF_long Ti [ ], UF_long Tj [ ],
double Tx [ ], double Tz [ ], UF_long Ap [ ], UF_long Ai [ ],	68	68	double Tx [ ], double Tz [ ], UF_long Ap [ ], UF_long Ai [ ],
double Ax [ ], double Az [ ], UF_long *status)	69	69	double Ax [ ], double Az [ ], UF_long *status)
{	70	70	{
status = umfpack_zl_triplet_to_col (m, n, nz, Ti, Tj, Tx, Tz, Ap, Ai, Ax, Az, (UF_long *) NULL);	71	71	status = umfpack_zl_triplet_to_col (m, n, nz, Ti, Tj, Tx, Tz, Ap, Ai, Ax, Az, (UF_long *) NULL);
}	72	72	}
	73	73
	74	74
	75	75
	76	76
/* complex(8) and integer(4) */	77	77	/* complex(8) and integer(4) */
	78	78
/* defaults */	79	79	/* defaults */
	80	80
void umfpack_zi_defaults_ (double Control [UMFPACK_CONTROL])	81	81	void umfpack_zi_defaults_ (double Control [UMFPACK_CONTROL])
{	82	82	{
umfpack_zi_defaults (Control) ;	83	83	umfpack_zi_defaults (Control) ;
}	84	84	}
	85	85
	86	86
/* free the Numeric object */	87	87	/* free the Numeric object */
void umfpack_zi_free_numeric_(void **Numeric)	88	88	void umfpack_zi_free_numeric_(void **Numeric)
{	89	89	{
umfpack_zi_free_numeric (Numeric) ;	90	90	umfpack_zi_free_numeric (Numeric) ;
}	91	91	}
	92	92
/* free the Symbolic object */	93	93	/* free the Symbolic object */
void umfpack_zi_free_symbolic_(void **Symbolic)	94	94	void umfpack_zi_free_symbolic_(void **Symbolic)
{	95	95	{
umfpack_zi_free_symbolic (Symbolic) ;	96	96	umfpack_zi_free_symbolic (Symbolic) ;
}	97	97	}
	98	98
/* numeric factorization */	99	99	/* numeric factorization */
void umfpack_zi_numeric_ (int Ap [ ], int Ai [ ], double Ax [ ], double Az [ ],	100	100	void umfpack_zi_numeric_ (int Ap [ ], int Ai [ ], double Ax [ ], double Az [ ],
void Symbolic, void Numeric,	101	101	void Symbolic, void Numeric,
double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])	102	102	double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
{	103	103	{
(void) umfpack_zi_numeric (Ap, Ai, Ax, Az, *Symbolic, Numeric, Control, Info);	104	104	(void) umfpack_zi_numeric (Ap, Ai, Ax, Az, *Symbolic, Numeric, Control, Info);
}	105	105	}
	106	106
/* pre-ordering and symbolic factorization */	107	107	/* pre-ordering and symbolic factorization */
void umfpack_zi_symbolic_ (int m, int n, int Ap [ ], int Ai [ ],	108	108	void umfpack_zi_symbolic_ (int m, int n, int Ap [ ], int Ai [ ],
double Ax [ ], double Az [ ], void **Symbolic,	109	109	double Ax [ ], double Az [ ], void **Symbolic,
double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])	110	110	double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
{	111	111	{
(void) umfpack_zi_symbolic (m, n, Ap, Ai, Ax, Az, Symbolic, Control, Info) ;	112	112	(void) umfpack_zi_symbolic (m, n, Ap, Ai, Ax, Az, Symbolic, Control, Info) ;
}	113	113	}
	114	114
/* solve a linear system */	115	115	/* solve a linear system */
	116	116
void umfpack_zi_solve_ (int *sys, int Ap [ ], int Ai [ ], double Ax [ ], double Az [ ],	117	117	void umfpack_zi_solve_ (int *sys, int Ap [ ], int Ai [ ], double Ax [ ], double Az [ ],
double x [ ], double xz [ ], double b [ ], double bz [ ], void **Numeric,	118	118	double x [ ], double xz [ ], double b [ ], double bz [ ], void **Numeric,
double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])	119	119	double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
{	120	120	{
(void) umfpack_zi_solve (*sys, Ap, Ai, Ax, Az, x, xz, b, bz,	121	121	(void) umfpack_zi_solve (*sys, Ap, Ai, Ax, Az, x, xz, b, bz,
*Numeric, Control, Info) ;	122	122	*Numeric, Control, Info) ;
}	123	123	}
	124	124
/* triplet 2 col */	125	125	/* triplet 2 col */
	126	126
void umfpack_zi_triplet_to_col_ (int m, int n, int *nz, int Ti [ ], int Tj [ ],	127	127	void umfpack_zi_triplet_to_col_ (int m, int n, int *nz, int Ti [ ], int Tj [ ],
double Tx [ ], double Tz [ ], int Ap [ ], int Ai [ ],	128	128	double Tx [ ], double Tz [ ], int Ap [ ], int Ai [ ],
double Ax [ ], double Az [ ], int *status)	129	129	double Ax [ ], double Az [ ], int *status)
{	130	130	{
status = umfpack_zi_triplet_to_col (m, n, nz, Ti, Tj, Tx, Tz, Ap, Ai, Ax, Az, (int *) NULL);	131	131	status = umfpack_zi_triplet_to_col (m, n, nz, Ti, Tj, Tx, Tz, Ap, Ai, Ax, Az, (int *) NULL);
}	132	132	}
		133
		134	/* real(8) and integer(8) */
		135
		136	/* defaults */
		137
		138	void umfpack_dl_defaults_ (double Control [UMFPACK_CONTROL])
		139	{
		140	umfpack_dl_defaults (Control) ;
		141	}
		142
		143	void umfpack_dl_free_numeric_ (void **Numeric)
		144	{
		145	umfpack_dl_free_numeric (Numeric) ;
		146	}
		147
		148	void umfpack_dl_free_symbolic_ (void **Symbolic)
		149	{
		150	umfpack_dl_free_symbolic (Symbolic) ;
		151	}
		152
		153	void umfpack_dl_numeric_ (UF_long Ap [ ], UF_long Ai [ ], double Ax [ ],
		154	void Symbolic, void Numeric,
		155	double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
		156	{
		157	(void) umfpack_dl_numeric (Ap, Ai, Ax, *Symbolic, Numeric, Control, Info);
		158	}
		159
		160	void umfpack_dl_symbolic_ (UF_long m, UF_long n, UF_long Ap [ ], UF_long Ai [ ],
		161	double Ax [ ], void **Symbolic,
		162	double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
		163	{
		164	(void) umfpack_dl_symbolic (m, n, Ap, Ai, Ax, Symbolic, Control, Info) ;
		165	}
		166
		167	void umfpack_dl_solve_ (UF_long *sys, UF_long Ap [ ], UF_long Ai [ ], double Ax [ ],
		168	double x [ ], double b [ ], void **Numeric,
		169	double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
		170	{
		171	(void) umfpack_dl_solve (sys, Ap, Ai, Ax, x, b, Numeric, Control, Info) ;
		172	}
		173
		174	void umfpack_dl_triplet_to_col_ (UF_long m, UF_long n, UF_long *nz, UF_long Ti [ ], UF_long Tj [ ],
		175	double T [ ], UF_long Ap [ ], UF_long Ai [ ],
		176	double A [ ], UF_long *status)
		177	{
		178	status = umfpack_dl_triplet_to_col (m, n, nz, Ti, Tj, T, Ap, Ai, A, (UF_long *) NULL);
		179	}
		180
		181	/* real(8) and integer(4) */
		182
		183	/* defaults */
		184
		185	void umfpack_di_defaults_ (double Control [UMFPACK_CONTROL])
		186	{
		187	umfpack_di_defaults (Control) ;
		188	}
		189
		190	void umfpack_di_free_numeric_ (void **Numeric)
		191	{
		192	umfpack_di_free_numeric (Numeric) ;
		193	}
		194
		195	void umfpack_di_free_symbolic_ (void **Symbolic)
		196	{
		197	umfpack_di_free_symbolic (Symbolic) ;

fvnlib.f90

Diff comments View file @ f85d3b3

	1	1
module fvn	2	2	module fvn
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	3	3	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	4	4	!
! fvn : a f95 module replacement for some imsl routines	5	5	! fvn : a f95 module replacement for some imsl routines
! it uses lapack for linear algebra	6	6	! it uses lapack for linear algebra
! it uses modified quadpack for integration	7	7	! it uses modified quadpack for integration
!	8	8	!
! William Daniau 2007	9	9	! William Daniau 2007
! william.daniau@femto-st.fr	10	10	! william.daniau@femto-st.fr
!	11	11	!
! Routines naming scheme :	12	12	! Routines naming scheme :
!	13	13	!
! fvn_x_name	14	14	! fvn_x_name
! where x can be s : real	15	15	! where x can be s : real
! d : real double precision	16	16	! d : real double precision
! c : complex	17	17	! c : complex
! z : double complex	18	18	! z : double complex
!	19	19	!
!	20	20	!
! This piece of code is totally free! Do whatever you want with it. However	21	21	! This piece of code is totally free! Do whatever you want with it. However
! if you find it usefull it would be kind to give credits ;-)	22	22	! if you find it usefull it would be kind to give credits ;-)
!	23	23	!
! svn version	24	24	! svn version
! September 2007 : added sparse system solving by interfacing umfpack	25	25	! September 2007 : added sparse system solving by interfacing umfpack
! June 2007 : added some complex trigonometric functions	26	26	! June 2007 : added some complex trigonometric functions
!	27	27	!
! TO DO LIST :	28	28	! TO DO LIST :
! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm	29	29	! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm
! eigenvalues are given with no particular order.	30	30	! eigenvalues are given with no particular order.
! + Generic interface for fvn_x_name family -> fvn_name	31	31	! + Generic interface for fvn_x_name family -> fvn_name
! + Make some parameters optional, status for example	32	32	! + Make some parameters optional, status for example
! + use f95 kinds "double complex" -> complex(kind=8)	33	33	! + use f95 kinds "double complex" -> complex(kind=8)
! + unify quadpack routines	34	34	! + unify quadpack routines
! + ...	35	35	! + ...
!	36	36	!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	37	37	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	38	38
implicit none	39	39	implicit none
! We define pi and i for the module	40	40	! We define pi and i for the module
real(kind=8),parameter :: fvn_pi = 3.141592653589793_8	41	41	real(kind=8),parameter :: fvn_pi = 3.141592653589793_8
complex(kind=8),parameter :: fvn_i = (0._8,1._8)	42	42	complex(kind=8),parameter :: fvn_i = (0._8,1._8)
	43	43
	44	44
! All quadpack routines are private to the module	45	45	! All quadpack routines are private to the module
private :: d1mach,dqag,dqag_2d_inner,dqag_2d_outer,dqage,dqage_2d_inner, &	46	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	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	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	49	dqk51,dqk51_2d_inner,dqk51_2d_outer,dqk61,dqk61_2d_inner,dqk61_2d_outer,dqpsrt
	50	50
	51	51
contains	52	52	contains
	53	53
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	54	54	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	55	55	!
! Matrix inversion subroutines	56	56	! Matrix inversion subroutines
!	57	57	!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	58	58	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine fvn_s_matinv(d,a,inva,status)	59	59	subroutine fvn_s_matinv(d,a,inva,status)
!	60	60	!
! Matrix inversion of a real matrix using BLAS and LAPACK	61	61	! Matrix inversion of a real matrix using BLAS and LAPACK
!	62	62	!
! d (in) : matrix rank	63	63	! d (in) : matrix rank
! a (in) : input matrix	64	64	! a (in) : input matrix
! inva (out) : inversed matrix	65	65	! inva (out) : inversed matrix
! status (ou) : =0 if something failed	66	66	! status (ou) : =0 if something failed
!	67	67	!
implicit none	68	68	implicit none
integer, intent(in) :: d	69	69	integer, intent(in) :: d
real, intent(in) :: a(d,d)	70	70	real, intent(in) :: a(d,d)
real, intent(out) :: inva(d,d)	71	71	real, intent(out) :: inva(d,d)
integer, intent(out) :: status	72	72	integer, intent(out) :: status
	73	73
integer, allocatable :: ipiv(:)	74	74	integer, allocatable :: ipiv(:)
real, allocatable :: work(:)	75	75	real, allocatable :: work(:)
real twork(1)	76	76	real twork(1)
integer :: info	77	77	integer :: info
integer :: lwork	78	78	integer :: lwork
	79	79
status=1	80	80	status=1
	81	81
allocate(ipiv(d))	82	82	allocate(ipiv(d))
! copy a into inva using BLAS	83	83	! copy a into inva using BLAS
!call scopy(d*d,a,1,inva,1)	84	84	!call scopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:)	85	85	inva(:,:)=a(:,:)
! LU factorization using LAPACK	86	86	! LU factorization using LAPACK
call sgetrf(d,d,inva,d,ipiv,info)	87	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	88	88	! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then	89	89	if (info /= 0) then
status=0	90	90	status=0
deallocate(ipiv)	91	91	deallocate(ipiv)
return	92	92	return
end if	93	93	end if
! we use the query fonction of xxxtri to obtain the optimal workspace size	94	94	! we use the query fonction of xxxtri to obtain the optimal workspace size
call sgetri(d,inva,d,ipiv,twork,-1,info)	95	95	call sgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1))	96	96	lwork=int(twork(1))
allocate(work(lwork))	97	97	allocate(work(lwork))
! Matrix inversion using LAPACK	98	98	! Matrix inversion using LAPACK
call sgetri(d,inva,d,ipiv,work,lwork,info)	99	99	call sgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0	100	100	! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then	101	101	if (info /= 0) then
status=0	102	102	status=0
end if	103	103	end if
deallocate(work)	104	104	deallocate(work)
deallocate(ipiv)	105	105	deallocate(ipiv)
end subroutine	106	106	end subroutine
	107	107
subroutine fvn_d_matinv(d,a,inva,status)	108	108	subroutine fvn_d_matinv(d,a,inva,status)
!	109	109	!
! Matrix inversion of a double precision matrix using BLAS and LAPACK	110	110	! Matrix inversion of a double precision matrix using BLAS and LAPACK
!	111	111	!
! d (in) : matrix rank	112	112	! d (in) : matrix rank
! a (in) : input matrix	113	113	! a (in) : input matrix
! inva (out) : inversed matrix	114	114	! inva (out) : inversed matrix
! status (ou) : =0 if something failed	115	115	! status (ou) : =0 if something failed
!	116	116	!
implicit none	117	117	implicit none
integer, intent(in) :: d	118	118	integer, intent(in) :: d
double precision, intent(in) :: a(d,d)	119	119	double precision, intent(in) :: a(d,d)
double precision, intent(out) :: inva(d,d)	120	120	double precision, intent(out) :: inva(d,d)
integer, intent(out) :: status	121	121	integer, intent(out) :: status
	122	122
integer, allocatable :: ipiv(:)	123	123	integer, allocatable :: ipiv(:)
double precision, allocatable :: work(:)	124	124	double precision, allocatable :: work(:)
double precision :: twork(1)	125	125	double precision :: twork(1)
integer :: info	126	126	integer :: info
integer :: lwork	127	127	integer :: lwork
	128	128
status=1	129	129	status=1
	130	130
allocate(ipiv(d))	131	131	allocate(ipiv(d))
! copy a into inva using BLAS	132	132	! copy a into inva using BLAS
!call dcopy(d*d,a,1,inva,1)	133	133	!call dcopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:)	134	134	inva(:,:)=a(:,:)
! LU factorization using LAPACK	135	135	! LU factorization using LAPACK
call dgetrf(d,d,inva,d,ipiv,info)	136	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	137	137	! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then	138	138	if (info /= 0) then
status=0	139	139	status=0
deallocate(ipiv)	140	140	deallocate(ipiv)
return	141	141	return
end if	142	142	end if
! we use the query fonction of xxxtri to obtain the optimal workspace size	143	143	! we use the query fonction of xxxtri to obtain the optimal workspace size
call dgetri(d,inva,d,ipiv,twork,-1,info)	144	144	call dgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1))	145	145	lwork=int(twork(1))
allocate(work(lwork))	146	146	allocate(work(lwork))
! Matrix inversion using LAPACK	147	147	! Matrix inversion using LAPACK
call dgetri(d,inva,d,ipiv,work,lwork,info)	148	148	call dgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0	149	149	! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then	150	150	if (info /= 0) then
status=0	151	151	status=0
end if	152	152	end if
deallocate(work)	153	153	deallocate(work)
deallocate(ipiv)	154	154	deallocate(ipiv)
end subroutine	155	155	end subroutine
	156	156
subroutine fvn_c_matinv(d,a,inva,status)	157	157	subroutine fvn_c_matinv(d,a,inva,status)
!	158	158	!
! Matrix inversion of a complex matrix using BLAS and LAPACK	159	159	! Matrix inversion of a complex matrix using BLAS and LAPACK
!	160	160	!
! d (in) : matrix rank	161	161	! d (in) : matrix rank
! a (in) : input matrix	162	162	! a (in) : input matrix
! inva (out) : inversed matrix	163	163	! inva (out) : inversed matrix
! status (ou) : =0 if something failed	164	164	! status (ou) : =0 if something failed
!	165	165	!
implicit none	166	166	implicit none
integer, intent(in) :: d	167	167	integer, intent(in) :: d
complex, intent(in) :: a(d,d)	168	168	complex, intent(in) :: a(d,d)
complex, intent(out) :: inva(d,d)	169	169	complex, intent(out) :: inva(d,d)
integer, intent(out) :: status	170	170	integer, intent(out) :: status
	171	171
integer, allocatable :: ipiv(:)	172	172	integer, allocatable :: ipiv(:)
complex, allocatable :: work(:)	173	173	complex, allocatable :: work(:)
complex :: twork(1)	174	174	complex :: twork(1)
integer :: info	175	175	integer :: info
integer :: lwork	176	176	integer :: lwork
	177	177
status=1	178	178	status=1
	179	179
allocate(ipiv(d))	180	180	allocate(ipiv(d))
! copy a into inva using BLAS	181	181	! copy a into inva using BLAS
!call ccopy(d*d,a,1,inva,1)	182	182	!call ccopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:)	183	183	inva(:,:)=a(:,:)
	184	184
! LU factorization using LAPACK	185	185	! LU factorization using LAPACK
call cgetrf(d,d,inva,d,ipiv,info)	186	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	187	187	! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then	188	188	if (info /= 0) then
status=0	189	189	status=0
deallocate(ipiv)	190	190	deallocate(ipiv)
return	191	191	return
end if	192	192	end if
! we use the query fonction of xxxtri to obtain the optimal workspace size	193	193	! we use the query fonction of xxxtri to obtain the optimal workspace size
call cgetri(d,inva,d,ipiv,twork,-1,info)	194	194	call cgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1))	195	195	lwork=int(twork(1))
allocate(work(lwork))	196	196	allocate(work(lwork))
! Matrix inversion using LAPACK	197	197	! Matrix inversion using LAPACK
call cgetri(d,inva,d,ipiv,work,lwork,info)	198	198	call cgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0	199	199	! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then	200	200	if (info /= 0) then
status=0	201	201	status=0
end if	202	202	end if
deallocate(work)	203	203	deallocate(work)
deallocate(ipiv)	204	204	deallocate(ipiv)
end subroutine	205	205	end subroutine
	206	206
subroutine fvn_z_matinv(d,a,inva,status)	207	207	subroutine fvn_z_matinv(d,a,inva,status)
!	208	208	!
! Matrix inversion of a double complex matrix using BLAS and LAPACK	209	209	! Matrix inversion of a double complex matrix using BLAS and LAPACK
!	210	210	!
! d (in) : matrix rank	211	211	! d (in) : matrix rank
! a (in) : input matrix	212	212	! a (in) : input matrix
! inva (out) : inversed matrix	213	213	! inva (out) : inversed matrix
! status (ou) : =0 if something failed	214	214	! status (ou) : =0 if something failed
!	215	215	!
implicit none	216	216	implicit none
integer, intent(in) :: d	217	217	integer, intent(in) :: d
double complex, intent(in) :: a(d,d)	218	218	double complex, intent(in) :: a(d,d)
double complex, intent(out) :: inva(d,d)	219	219	double complex, intent(out) :: inva(d,d)
integer, intent(out) :: status	220	220	integer, intent(out) :: status
	221	221
integer, allocatable :: ipiv(:)	222	222	integer, allocatable :: ipiv(:)
double complex, allocatable :: work(:)	223	223	double complex, allocatable :: work(:)
double complex :: twork(1)	224	224	double complex :: twork(1)
integer :: info	225	225	integer :: info
integer :: lwork	226	226	integer :: lwork
	227	227
status=1	228	228	status=1
	229	229
allocate(ipiv(d))	230	230	allocate(ipiv(d))
! copy a into inva using BLAS	231	231	! copy a into inva using BLAS
!call zcopy(d*d,a,1,inva,1)	232	232	!call zcopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:)	233	233	inva(:,:)=a(:,:)
	234	234
! LU factorization using LAPACK	235	235	! LU factorization using LAPACK
call zgetrf(d,d,inva,d,ipiv,info)	236	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	237	237	! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then	238	238	if (info /= 0) then
status=0	239	239	status=0
deallocate(ipiv)	240	240	deallocate(ipiv)
return	241	241	return
end if	242	242	end if
! we use the query fonction of xxxtri to obtain the optimal workspace size	243	243	! we use the query fonction of xxxtri to obtain the optimal workspace size
call zgetri(d,inva,d,ipiv,twork,-1,info)	244	244	call zgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1))	245	245	lwork=int(twork(1))
allocate(work(lwork))	246	246	allocate(work(lwork))
! Matrix inversion using LAPACK	247	247	! Matrix inversion using LAPACK
call zgetri(d,inva,d,ipiv,work,lwork,info)	248	248	call zgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0	249	249	! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then	250	250	if (info /= 0) then
status=0	251	251	status=0
end if	252	252	end if
deallocate(work)	253	253	deallocate(work)
deallocate(ipiv)	254	254	deallocate(ipiv)
end subroutine	255	255	end subroutine
	256	256
	257	257
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	258	258	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	259	259	!
! Determinants	260	260	! Determinants
!	261	261	!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	262	262	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function fvn_s_det(d,a,status)	263	263	function fvn_s_det(d,a,status)
!	264	264	!
! Evaluate the determinant of a square matrix using lapack LU factorization	265	265	! Evaluate the determinant of a square matrix using lapack LU factorization
!	266	266	!
! d (in) : matrix rank	267	267	! d (in) : matrix rank
! a (in) : The Matrix	268	268	! a (in) : The Matrix
! status (out) : =0 if LU factorization failed	269	269	! status (out) : =0 if LU factorization failed
!	270	270	!
implicit none	271	271	implicit none
integer, intent(in) :: d	272	272	integer, intent(in) :: d
real, intent(in) :: a(d,d)	273	273	real, intent(in) :: a(d,d)
integer, intent(out) :: status	274	274	integer, intent(out) :: status
real :: fvn_s_det	275	275	real :: fvn_s_det
	276	276
real, allocatable :: wc_a(:,:)	277	277	real, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:)	278	278	integer, allocatable :: ipiv(:)
integer :: info,i	279	279	integer :: info,i
	280	280
status=1	281	281	status=1
allocate(wc_a(d,d))	282	282	allocate(wc_a(d,d))
allocate(ipiv(d))	283	283	allocate(ipiv(d))
wc_a(:,:)=a(:,:)	284	284	wc_a(:,:)=a(:,:)
call sgetrf(d,d,wc_a,d,ipiv,info)	285	285	call sgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then	286	286	if (info/= 0) then
status=0	287	287	status=0
fvn_s_det=0.e0	288	288	fvn_s_det=0.e0
deallocate(ipiv)	289	289	deallocate(ipiv)
deallocate(wc_a)	290	290	deallocate(wc_a)
return	291	291	return
end if	292	292	end if
fvn_s_det=1.e0	293	293	fvn_s_det=1.e0
do i=1,d	294	294	do i=1,d
if (ipiv(i)==i) then	295	295	if (ipiv(i)==i) then
fvn_s_det=fvn_s_det*wc_a(i,i)	296	296	fvn_s_det=fvn_s_det*wc_a(i,i)
else	297	297	else
fvn_s_det=-fvn_s_det*wc_a(i,i)	298	298	fvn_s_det=-fvn_s_det*wc_a(i,i)
end if	299	299	end if
end do	300	300	end do
deallocate(ipiv)	301	301	deallocate(ipiv)
deallocate(wc_a)	302	302	deallocate(wc_a)
	303	303
end function	304	304	end function
	305	305
function fvn_d_det(d,a,status)	306	306	function fvn_d_det(d,a,status)
!	307	307	!
! Evaluate the determinant of a square matrix using lapack LU factorization	308	308	! Evaluate the determinant of a square matrix using lapack LU factorization
!	309	309	!
! d (in) : matrix rank	310	310	! d (in) : matrix rank
! a (in) : The Matrix	311	311	! a (in) : The Matrix
! status (out) : =0 if LU factorization failed	312	312	! status (out) : =0 if LU factorization failed
!	313	313	!
implicit none	314	314	implicit none
integer, intent(in) :: d	315	315	integer, intent(in) :: d
double precision, intent(in) :: a(d,d)	316	316	double precision, intent(in) :: a(d,d)
integer, intent(out) :: status	317	317	integer, intent(out) :: status
double precision :: fvn_d_det	318	318	double precision :: fvn_d_det
	319	319
double precision, allocatable :: wc_a(:,:)	320	320	double precision, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:)	321	321	integer, allocatable :: ipiv(:)
integer :: info,i	322	322	integer :: info,i
	323	323
status=1	324	324	status=1
allocate(wc_a(d,d))	325	325	allocate(wc_a(d,d))
allocate(ipiv(d))	326	326	allocate(ipiv(d))
wc_a(:,:)=a(:,:)	327	327	wc_a(:,:)=a(:,:)
call dgetrf(d,d,wc_a,d,ipiv,info)	328	328	call dgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then	329	329	if (info/= 0) then
status=0	330	330	status=0
fvn_d_det=0.d0	331	331	fvn_d_det=0.d0
deallocate(ipiv)	332	332	deallocate(ipiv)
deallocate(wc_a)	333	333	deallocate(wc_a)
return	334	334	return
end if	335	335	end if
fvn_d_det=1.d0	336	336	fvn_d_det=1.d0
do i=1,d	337	337	do i=1,d
if (ipiv(i)==i) then	338	338	if (ipiv(i)==i) then
fvn_d_det=fvn_d_det*wc_a(i,i)	339	339	fvn_d_det=fvn_d_det*wc_a(i,i)
else	340	340	else
fvn_d_det=-fvn_d_det*wc_a(i,i)	341	341	fvn_d_det=-fvn_d_det*wc_a(i,i)
end if	342	342	end if
end do	343	343	end do
deallocate(ipiv)	344	344	deallocate(ipiv)
deallocate(wc_a)	345	345	deallocate(wc_a)
	346	346
end function	347	347	end function
	348	348
function fvn_c_det(d,a,status) !	349	349	function fvn_c_det(d,a,status) !
! Evaluate the determinant of a square matrix using lapack LU factorization	350	350	! Evaluate the determinant of a square matrix using lapack LU factorization
!	351	351	!
! d (in) : matrix rank	352	352	! d (in) : matrix rank
! a (in) : The Matrix	353	353	! a (in) : The Matrix
! status (out) : =0 if LU factorization failed	354	354	! status (out) : =0 if LU factorization failed
!	355	355	!
implicit none	356	356	implicit none
integer, intent(in) :: d	357	357	integer, intent(in) :: d
complex, intent(in) :: a(d,d)	358	358	complex, intent(in) :: a(d,d)
integer, intent(out) :: status	359	359	integer, intent(out) :: status
complex :: fvn_c_det	360	360	complex :: fvn_c_det
	361	361
complex, allocatable :: wc_a(:,:)	362	362	complex, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:)	363	363	integer, allocatable :: ipiv(:)
integer :: info,i	364	364	integer :: info,i
	365	365
status=1	366	366	status=1
allocate(wc_a(d,d))	367	367	allocate(wc_a(d,d))
allocate(ipiv(d))	368	368	allocate(ipiv(d))
wc_a(:,:)=a(:,:)	369	369	wc_a(:,:)=a(:,:)
call cgetrf(d,d,wc_a,d,ipiv,info)	370	370	call cgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then	371	371	if (info/= 0) then
status=0	372	372	status=0
fvn_c_det=(0.e0,0.e0)	373	373	fvn_c_det=(0.e0,0.e0)
deallocate(ipiv)	374	374	deallocate(ipiv)
deallocate(wc_a)	375	375	deallocate(wc_a)
return	376	376	return
end if	377	377	end if
fvn_c_det=(1.e0,0.e0)	378	378	fvn_c_det=(1.e0,0.e0)
do i=1,d	379	379	do i=1,d
if (ipiv(i)==i) then	380	380	if (ipiv(i)==i) then
fvn_c_det=fvn_c_det*wc_a(i,i)	381	381	fvn_c_det=fvn_c_det*wc_a(i,i)
else	382	382	else
fvn_c_det=-fvn_c_det*wc_a(i,i)	383	383	fvn_c_det=-fvn_c_det*wc_a(i,i)
end if	384	384	end if
end do	385	385	end do
deallocate(ipiv)	386	386	deallocate(ipiv)
deallocate(wc_a)	387	387	deallocate(wc_a)
	388	388
end function	389	389	end function
	390	390
function fvn_z_det(d,a,status)	391	391	function fvn_z_det(d,a,status)
!	392	392	!
! Evaluate the determinant of a square matrix using lapack LU factorization	393	393	! Evaluate the determinant of a square matrix using lapack LU factorization
!	394	394	!
! d (in) : matrix rank	395	395	! d (in) : matrix rank
! a (in) : The Matrix	396	396	! a (in) : The Matrix
! det (out) : determinant	397	397	! det (out) : determinant
! status (out) : =0 if LU factorization failed	398	398	! status (out) : =0 if LU factorization failed
!	399	399	!
implicit none	400	400	implicit none
integer, intent(in) :: d	401	401	integer, intent(in) :: d
double complex, intent(in) :: a(d,d)	402	402	double complex, intent(in) :: a(d,d)
integer, intent(out) :: status	403	403	integer, intent(out) :: status
double complex :: fvn_z_det	404	404	double complex :: fvn_z_det
	405	405
double complex, allocatable :: wc_a(:,:)	406	406	double complex, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:)	407	407	integer, allocatable :: ipiv(:)
integer :: info,i	408	408	integer :: info,i
	409	409
status=1	410	410	status=1
allocate(wc_a(d,d))	411	411	allocate(wc_a(d,d))
allocate(ipiv(d))	412	412	allocate(ipiv(d))
wc_a(:,:)=a(:,:)	413	413	wc_a(:,:)=a(:,:)
call zgetrf(d,d,wc_a,d,ipiv,info)	414	414	call zgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then	415	415	if (info/= 0) then
status=0	416	416	status=0
fvn_z_det=(0.d0,0.d0)	417	417	fvn_z_det=(0.d0,0.d0)
deallocate(ipiv)	418	418	deallocate(ipiv)
deallocate(wc_a)	419	419	deallocate(wc_a)
return	420	420	return
end if	421	421	end if
fvn_z_det=(1.d0,0.d0)	422	422	fvn_z_det=(1.d0,0.d0)
do i=1,d	423	423	do i=1,d
if (ipiv(i)==i) then	424	424	if (ipiv(i)==i) then
fvn_z_det=fvn_z_det*wc_a(i,i)	425	425	fvn_z_det=fvn_z_det*wc_a(i,i)
else	426	426	else
fvn_z_det=-fvn_z_det*wc_a(i,i)	427	427	fvn_z_det=-fvn_z_det*wc_a(i,i)
end if	428	428	end if
end do	429	429	end do
deallocate(ipiv)	430	430	deallocate(ipiv)
deallocate(wc_a)	431	431	deallocate(wc_a)
	432	432
end function	433	433	end function
	434	434
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	435	435	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	436	436	!
! Condition test	437	437	! Condition test
!	438	438	!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	439	439	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1-norm	440	440	! 1-norm
! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm	441	441	! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm
! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond	442	442	! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond
!	443	443	!
subroutine fvn_s_matcon(d,a,rcond,status)	444	444	subroutine fvn_s_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number)	445	445	! Matrix condition (reciprocal of condition number)
!	446	446	!
! d (in) : matrix rank	447	447	! d (in) : matrix rank
! a (in) : The Matrix	448	448	! a (in) : The Matrix
! rcond (out) : guess what	449	449	! rcond (out) : guess what
! status (out) : =0 if something went wrong	450	450	! status (out) : =0 if something went wrong
!	451	451	!
implicit none	452	452	implicit none
integer, intent(in) :: d	453	453	integer, intent(in) :: d
real, intent(in) :: a(d,d)	454	454	real, intent(in) :: a(d,d)
real, intent(out) :: rcond	455	455	real, intent(out) :: rcond
integer, intent(out) :: status	456	456	integer, intent(out) :: status
	457	457
real, allocatable :: work(:)	458	458	real, allocatable :: work(:)
integer, allocatable :: iwork(:)	459	459	integer, allocatable :: iwork(:)
real :: anorm	460	460	real :: anorm
real, allocatable :: wc_a(:,:) ! working copy of a	461	461	real, allocatable :: wc_a(:,:) ! working copy of a
integer :: info	462	462	integer :: info
integer, allocatable :: ipiv(:)	463	463	integer, allocatable :: ipiv(:)
	464	464
real, external :: slange	465	465	real, external :: slange
	466	466
	467	467
status=1	468	468	status=1
	469	469
anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm	470	470	anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
	471	471
allocate(wc_a(d,d))	472	472	allocate(wc_a(d,d))
!call scopy(d*d,a,1,wc_a,1)	473	473	!call scopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:)	474	474	wc_a(:,:)=a(:,:)
	475	475
allocate(ipiv(d))	476	476	allocate(ipiv(d))
call sgetrf(d,d,wc_a,d,ipiv,info)	477	477	call sgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then	478	478	if (info /= 0) then
status=0	479	479	status=0
deallocate(ipiv)	480	480	deallocate(ipiv)
deallocate(wc_a)	481	481	deallocate(wc_a)
return	482	482	return
end if	483	483	end if
allocate(work(4*d))	484	484	allocate(work(4*d))
allocate(iwork(d))	485	485	allocate(iwork(d))
call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)	486	486	call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
if (info /= 0) then	487	487	if (info /= 0) then
status=0	488	488	status=0
end if	489	489	end if
deallocate(iwork)	490	490	deallocate(iwork)
deallocate(work)	491	491	deallocate(work)
deallocate(ipiv)	492	492	deallocate(ipiv)
deallocate(wc_a)	493	493	deallocate(wc_a)
	494	494
end subroutine	495	495	end subroutine
	496	496
subroutine fvn_d_matcon(d,a,rcond,status)	497	497	subroutine fvn_d_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number)	498	498	! Matrix condition (reciprocal of condition number)
!	499	499	!
! d (in) : matrix rank	500	500	! d (in) : matrix rank
! a (in) : The Matrix	501	501	! a (in) : The Matrix
! rcond (out) : guess what	502	502	! rcond (out) : guess what
! status (out) : =0 if something went wrong	503	503	! status (out) : =0 if something went wrong
!	504	504	!
implicit none	505	505	implicit none
integer, intent(in) :: d	506	506	integer, intent(in) :: d
double precision, intent(in) :: a(d,d)	507	507	double precision, intent(in) :: a(d,d)
double precision, intent(out) :: rcond	508	508	double precision, intent(out) :: rcond
integer, intent(out) :: status	509	509	integer, intent(out) :: status
	510	510
double precision, allocatable :: work(:)	511	511	double precision, allocatable :: work(:)
integer, allocatable :: iwork(:)	512	512	integer, allocatable :: iwork(:)
double precision :: anorm	513	513	double precision :: anorm
double precision, allocatable :: wc_a(:,:) ! working copy of a	514	514	double precision, allocatable :: wc_a(:,:) ! working copy of a
integer :: info	515	515	integer :: info
integer, allocatable :: ipiv(:)	516	516	integer, allocatable :: ipiv(:)
	517	517
double precision, external :: dlange	518	518	double precision, external :: dlange
	519	519
	520	520
status=1	521	521	status=1
	522	522
anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm	523	523	anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
	524	524
allocate(wc_a(d,d))	525	525	allocate(wc_a(d,d))
!call dcopy(d*d,a,1,wc_a,1)	526	526	!call dcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:)	527	527	wc_a(:,:)=a(:,:)
	528	528
allocate(ipiv(d))	529	529	allocate(ipiv(d))
call dgetrf(d,d,wc_a,d,ipiv,info)	530	530	call dgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then	531	531	if (info /= 0) then
status=0	532	532	status=0
deallocate(ipiv)	533	533	deallocate(ipiv)
deallocate(wc_a)	534	534	deallocate(wc_a)
return	535	535	return
end if	536	536	end if
	537	537
allocate(work(4*d))	538	538	allocate(work(4*d))
allocate(iwork(d))	539	539	allocate(iwork(d))
call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)	540	540	call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
if (info /= 0) then	541	541	if (info /= 0) then
status=0	542	542	status=0
end if	543	543	end if
deallocate(iwork)	544	544	deallocate(iwork)
deallocate(work)	545	545	deallocate(work)
deallocate(ipiv)	546	546	deallocate(ipiv)
deallocate(wc_a)	547	547	deallocate(wc_a)
	548	548
end subroutine	549	549	end subroutine
	550	550
subroutine fvn_c_matcon(d,a,rcond,status)	551	551	subroutine fvn_c_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number)	552	552	! Matrix condition (reciprocal of condition number)
!	553	553	!
! d (in) : matrix rank	554	554	! d (in) : matrix rank
! a (in) : The Matrix	555	555	! a (in) : The Matrix
! rcond (out) : guess what	556	556	! rcond (out) : guess what
! status (out) : =0 if something went wrong	557	557	! status (out) : =0 if something went wrong
!	558	558	!
implicit none	559	559	implicit none
integer, intent(in) :: d	560	560	integer, intent(in) :: d
complex, intent(in) :: a(d,d)	561	561	complex, intent(in) :: a(d,d)
real, intent(out) :: rcond	562	562	real, intent(out) :: rcond
integer, intent(out) :: status	563	563	integer, intent(out) :: status
	564	564
real, allocatable :: rwork(:)	565	565	real, allocatable :: rwork(:)
complex, allocatable :: work(:)	566	566	complex, allocatable :: work(:)
integer, allocatable :: iwork(:)	567	567	integer, allocatable :: iwork(:)
real :: anorm	568	568	real :: anorm
complex, allocatable :: wc_a(:,:) ! working copy of a	569	569	complex, allocatable :: wc_a(:,:) ! working copy of a
integer :: info	570	570	integer :: info
integer, allocatable :: ipiv(:)	571	571	integer, allocatable :: ipiv(:)
	572	572
real, external :: clange	573	573	real, external :: clange
	574	574
	575	575
status=1	576	576	status=1
	577	577
anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm	578	578	anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
	579	579
allocate(wc_a(d,d))	580	580	allocate(wc_a(d,d))
!call ccopy(d*d,a,1,wc_a,1)	581	581	!call ccopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:)	582	582	wc_a(:,:)=a(:,:)
	583	583
allocate(ipiv(d))	584	584	allocate(ipiv(d))
call cgetrf(d,d,wc_a,d,ipiv,info)	585	585	call cgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then	586	586	if (info /= 0) then
status=0	587	587	status=0
deallocate(ipiv)	588	588	deallocate(ipiv)
deallocate(wc_a)	589	589	deallocate(wc_a)
return	590	590	return
end if	591	591	end if
allocate(work(2*d))	592	592	allocate(work(2*d))
allocate(rwork(2*d))	593	593	allocate(rwork(2*d))
call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)	594	594	call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
if (info /= 0) then	595	595	if (info /= 0) then
status=0	596	596	status=0
end if	597	597	end if
deallocate(rwork)	598	598	deallocate(rwork)
deallocate(work)	599	599	deallocate(work)
deallocate(ipiv)	600	600	deallocate(ipiv)
deallocate(wc_a)	601	601	deallocate(wc_a)
end subroutine	602	602	end subroutine
	603	603
subroutine fvn_z_matcon(d,a,rcond,status)	604	604	subroutine fvn_z_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number)	605	605	! Matrix condition (reciprocal of condition number)
!	606	606	!
! d (in) : matrix rank	607	607	! d (in) : matrix rank
! a (in) : The Matrix	608	608	! a (in) : The Matrix
! rcond (out) : guess what	609	609	! rcond (out) : guess what
! status (out) : =0 if something went wrong	610	610	! status (out) : =0 if something went wrong
!	611	611	!
implicit none	612	612	implicit none
integer, intent(in) :: d	613	613	integer, intent(in) :: d
double complex, intent(in) :: a(d,d)	614	614	double complex, intent(in) :: a(d,d)
double precision, intent(out) :: rcond	615	615	double precision, intent(out) :: rcond
integer, intent(out) :: status	616	616	integer, intent(out) :: status
	617	617
double complex, allocatable :: work(:)	618	618	double complex, allocatable :: work(:)
double precision, allocatable :: rwork(:)	619	619	double precision, allocatable :: rwork(:)
double precision :: anorm	620	620	double precision :: anorm
double complex, allocatable :: wc_a(:,:) ! working copy of a	621	621	double complex, allocatable :: wc_a(:,:) ! working copy of a
integer :: info	622	622	integer :: info
integer, allocatable :: ipiv(:)	623	623	integer, allocatable :: ipiv(:)
	624	624
double precision, external :: zlange	625	625	double precision, external :: zlange
	626	626
	627	627
status=1	628	628	status=1
	629	629
anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm	630	630	anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
	631	631
allocate(wc_a(d,d))	632	632	allocate(wc_a(d,d))
!call zcopy(d*d,a,1,wc_a,1)	633	633	!call zcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:)	634	634	wc_a(:,:)=a(:,:)
	635	635
allocate(ipiv(d))	636	636	allocate(ipiv(d))
call zgetrf(d,d,wc_a,d,ipiv,info)	637	637	call zgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then	638	638	if (info /= 0) then
status=0	639	639	status=0
deallocate(ipiv)	640	640	deallocate(ipiv)
deallocate(wc_a)	641	641	deallocate(wc_a)
return	642	642	return
end if	643	643	end if
	644	644
allocate(work(2*d))	645	645	allocate(work(2*d))
allocate(rwork(2*d))	646	646	allocate(rwork(2*d))
call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)	647	647	call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
if (info /= 0) then	648	648	if (info /= 0) then
status=0	649	649	status=0
end if	650	650	end if
deallocate(rwork)	651	651	deallocate(rwork)
deallocate(work)	652	652	deallocate(work)
deallocate(ipiv)	653	653	deallocate(ipiv)
deallocate(wc_a)	654	654	deallocate(wc_a)
end subroutine	655	655	end subroutine
	656	656
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	657	657	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	658	658	!
! Valeurs propres/ Vecteurs propre	659	659	! Valeurs propres/ Vecteurs propre
!	660	660	!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	661	661	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	662	662
subroutine fvn_s_matev(d,a,evala,eveca,status)	663	663	subroutine fvn_s_matev(d,a,evala,eveca,status)
!	664	664	!
! integer d (in) : matrice rank	665	665	! integer d (in) : matrice rank
! real a(d,d) (in) : The Matrix	666	666	! real a(d,d) (in) : The Matrix
! complex evala(d) (out) : eigenvalues	667	667	! complex evala(d) (out) : eigenvalues
! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector	668	668	! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong	669	669	! integer (out) : status =0 if something went wrong
!	670	670	!
! interfacing Lapack routine SGEEV	671	671	! interfacing Lapack routine SGEEV
implicit none	672	672	implicit none
integer, intent(in) :: d	673	673	integer, intent(in) :: d
real, intent(in) :: a(d,d)	674	674	real, intent(in) :: a(d,d)
complex, intent(out) :: evala(d)	675	675	complex, intent(out) :: evala(d)
complex, intent(out) :: eveca(d,d)	676	676	complex, intent(out) :: eveca(d,d)
integer, intent(out) :: status	677	677	integer, intent(out) :: status
	678	678
real, allocatable :: wc_a(:,:) ! a working copy of a	679	679	real, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info	680	680	integer :: info
integer :: lwork	681	681	integer :: lwork
real, allocatable :: wr(:),wi(:)	682	682	real, allocatable :: wr(:),wi(:)
real :: vl ! unused but necessary for the call	683	683	real :: vl ! unused but necessary for the call
real, allocatable :: vr(:,:)	684	684	real, allocatable :: vr(:,:)
real, allocatable :: work(:)	685	685	real, allocatable :: work(:)
real :: twork(1)	686	686	real :: twork(1)
integer i	687	687	integer i
integer j	688	688	integer j
	689	689
! making a working copy of a	690	690	! making a working copy of a
allocate(wc_a(d,d))	691	691	allocate(wc_a(d,d))
!call scopy(d*d,a,1,wc_a,1)	692	692	!call scopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:)	693	693	wc_a(:,:)=a(:,:)
	694	694
allocate(wr(d))	695	695	allocate(wr(d))
allocate(wi(d))	696	696	allocate(wi(d))
allocate(vr(d,d))	697	697	allocate(vr(d,d))
! query optimal work size	698	698	! query optimal work size
call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)	699	699	call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
lwork=int(twork(1))	700	700	lwork=int(twork(1))
allocate(work(lwork))	701	701	allocate(work(lwork))
call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)	702	702	call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
	703	703
if (info /= 0) then	704	704	if (info /= 0) then
status=0	705	705	status=0
deallocate(work)	706	706	deallocate(work)
deallocate(vr)	707	707	deallocate(vr)
deallocate(wi)	708	708	deallocate(wi)
deallocate(wr)	709	709	deallocate(wr)
deallocate(wc_a)	710	710	deallocate(wc_a)
return	711	711	return
end if	712	712	end if
	713	713
! now fill in the results	714	714	! now fill in the results
i=1	715	715	i=1
do while(i<=d)	716	716	do while(i<=d)
evala(i)=cmplx(wr(i),wi(i))	717	717	evala(i)=cmplx(wr(i),wi(i))
if (wi(i) == 0.) then ! eigenvalue is real	718	718	if (wi(i) == 0.) then ! eigenvalue is real
eveca(:,i)=cmplx(vr(:,i),0.)	719	719	eveca(:,i)=cmplx(vr(:,i),0.)
else ! eigenvalue is complex	720	720	else ! eigenvalue is complex
evala(i+1)=cmplx(wr(i+1),wi(i+1))	721	721	evala(i+1)=cmplx(wr(i+1),wi(i+1))
eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))	722	722	eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))
eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))	723	723	eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))
i=i+1	724	724	i=i+1
end if	725	725	end if
i=i+1	726	726	i=i+1
enddo	727	727	enddo
deallocate(work)	728	728	deallocate(work)
deallocate(vr)	729	729	deallocate(vr)
deallocate(wi)	730	730	deallocate(wi)
deallocate(wr)	731	731	deallocate(wr)
deallocate(wc_a)	732	732	deallocate(wc_a)
	733	733
end subroutine	734	734	end subroutine
	735	735
subroutine fvn_d_matev(d,a,evala,eveca,status)	736	736	subroutine fvn_d_matev(d,a,evala,eveca,status)
!	737	737	!
! integer d (in) : matrice rank	738	738	! integer d (in) : matrice rank
! double precision a(d,d) (in) : The Matrix	739	739	! double precision a(d,d) (in) : The Matrix
! double complex evala(d) (out) : eigenvalues	740	740	! double complex evala(d) (out) : eigenvalues
! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector	741	741	! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong	742	742	! integer (out) : status =0 if something went wrong
!	743	743	!
! interfacing Lapack routine DGEEV	744	744	! interfacing Lapack routine DGEEV
implicit none	745	745	implicit none
integer, intent(in) :: d	746	746	integer, intent(in) :: d
double precision, intent(in) :: a(d,d)	747	747	double precision, intent(in) :: a(d,d)
double complex, intent(out) :: evala(d)	748	748	double complex, intent(out) :: evala(d)
double complex, intent(out) :: eveca(d,d)	749	749	double complex, intent(out) :: eveca(d,d)
integer, intent(out) :: status	750	750	integer, intent(out) :: status
	751	751
double precision, allocatable :: wc_a(:,:) ! a working copy of a	752	752	double precision, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info	753	753	integer :: info
integer :: lwork	754	754	integer :: lwork
double precision, allocatable :: wr(:),wi(:)	755	755	double precision, allocatable :: wr(:),wi(:)
double precision :: vl ! unused but necessary for the call	756	756	double precision :: vl ! unused but necessary for the call
double precision, allocatable :: vr(:,:)	757	757	double precision, allocatable :: vr(:,:)
double precision, allocatable :: work(:)	758	758	double precision, allocatable :: work(:)
double precision :: twork(1)	759	759	double precision :: twork(1)
integer i	760	760	integer i
integer j	761	761	integer j
	762	762
! making a working copy of a	763	763	! making a working copy of a
allocate(wc_a(d,d))	764	764	allocate(wc_a(d,d))
!call dcopy(d*d,a,1,wc_a,1)	765	765	!call dcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:)	766	766	wc_a(:,:)=a(:,:)
	767	767
allocate(wr(d))	768	768	allocate(wr(d))
allocate(wi(d))	769	769	allocate(wi(d))
allocate(vr(d,d))	770	770	allocate(vr(d,d))
! query optimal work size	771	771	! query optimal work size
call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)	772	772	call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
lwork=int(twork(1))	773	773	lwork=int(twork(1))
allocate(work(lwork))	774	774	allocate(work(lwork))
call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)	775	775	call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
	776	776
if (info /= 0) then	777	777	if (info /= 0) then
status=0	778	778	status=0
deallocate(work)	779	779	deallocate(work)
deallocate(vr)	780	780	deallocate(vr)
deallocate(wi)	781	781	deallocate(wi)
deallocate(wr)	782	782	deallocate(wr)
deallocate(wc_a)	783	783	deallocate(wc_a)
return	784	784	return
end if	785	785	end if
	786	786
! now fill in the results	787	787	! now fill in the results
i=1	788	788	i=1
do while(i<=d)	789	789	do while(i<=d)
evala(i)=dcmplx(wr(i),wi(i))	790	790	evala(i)=dcmplx(wr(i),wi(i))
if (wi(i) == 0.) then ! eigenvalue is real	791	791	if (wi(i) == 0.) then ! eigenvalue is real
eveca(:,i)=dcmplx(vr(:,i),0.)	792	792	eveca(:,i)=dcmplx(vr(:,i),0.)
else ! eigenvalue is complex	793	793	else ! eigenvalue is complex
evala(i+1)=dcmplx(wr(i+1),wi(i+1))	794	794	evala(i+1)=dcmplx(wr(i+1),wi(i+1))
eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))	795	795	eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))
eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))	796	796	eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))
i=i+1	797	797	i=i+1
end if	798	798	end if
i=i+1	799	799	i=i+1
enddo	800	800	enddo
	801	801
deallocate(work)	802	802	deallocate(work)
deallocate(vr)	803	803	deallocate(vr)
deallocate(wi)	804	804	deallocate(wi)
deallocate(wr)	805	805	deallocate(wr)
deallocate(wc_a)	806	806	deallocate(wc_a)
	807	807
end subroutine	808	808	end subroutine
	809	809
subroutine fvn_c_matev(d,a,evala,eveca,status)	810	810	subroutine fvn_c_matev(d,a,evala,eveca,status)
!	811	811	!
! integer d (in) : matrice rank	812	812	! integer d (in) : matrice rank
! complex a(d,d) (in) : The Matrix	813	813	! complex a(d,d) (in) : The Matrix
! complex evala(d) (out) : eigenvalues	814	814	! complex evala(d) (out) : eigenvalues
! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector	815	815	! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong	816	816	! integer (out) : status =0 if something went wrong
!	817	817	!
! interfacing Lapack routine CGEEV	818	818	! interfacing Lapack routine CGEEV
implicit none	819	819	implicit none
integer, intent(in) :: d	820	820	integer, intent(in) :: d
complex, intent(in) :: a(d,d)	821	821	complex, intent(in) :: a(d,d)
complex, intent(out) :: evala(d)	822	822	complex, intent(out) :: evala(d)
complex, intent(out) :: eveca(d,d)	823	823	complex, intent(out) :: eveca(d,d)
integer, intent(out) :: status	824	824	integer, intent(out) :: status
	825	825
complex, allocatable :: wc_a(:,:) ! a working copy of a	826	826	complex, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info	827	827	integer :: info
integer :: lwork	828	828	integer :: lwork
complex, allocatable :: work(:)	829	829	complex, allocatable :: work(:)
complex :: twork(1)	830	830	complex :: twork(1)
real, allocatable :: rwork(:)	831	831	real, allocatable :: rwork(:)
complex :: vl ! unused but necessary for the call	832	832	complex :: vl ! unused but necessary for the call
	833	833
status=1	834	834	status=1
	835	835
! making a working copy of a	836	836	! making a working copy of a
allocate(wc_a(d,d))	837	837	allocate(wc_a(d,d))
!call ccopy(d*d,a,1,wc_a,1)	838	838	!call ccopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:)	839	839	wc_a(:,:)=a(:,:)
	840	840
	841	841
! query optimal work size	842	842	! query optimal work size
call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)	843	843	call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
lwork=int(twork(1))	844	844	lwork=int(twork(1))
allocate(work(lwork))	845	845	allocate(work(lwork))
allocate(rwork(2*d))	846	846	allocate(rwork(2*d))
call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)	847	847	call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
	848	848
if (info /= 0) then	849	849	if (info /= 0) then
status=0	850	850	status=0
end if	851	851	end if
deallocate(rwork)	852	852	deallocate(rwork)
deallocate(work)	853	853	deallocate(work)
deallocate(wc_a)	854	854	deallocate(wc_a)
	855	855
end subroutine	856	856	end subroutine
	857	857
subroutine fvn_z_matev(d,a,evala,eveca,status)	858	858	subroutine fvn_z_matev(d,a,evala,eveca,status)
!	859	859	!
! integer d (in) : matrice rank	860	860	! integer d (in) : matrice rank
! double complex a(d,d) (in) : The Matrix	861	861	! double complex a(d,d) (in) : The Matrix
! double complex evala(d) (out) : eigenvalues	862	862	! double complex evala(d) (out) : eigenvalues
! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector	863	863	! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong	864	864	! integer (out) : status =0 if something went wrong
!	865	865	!
! interfacing Lapack routine ZGEEV	866	866	! interfacing Lapack routine ZGEEV
implicit none	867	867	implicit none
integer, intent(in) :: d	868	868	integer, intent(in) :: d
double complex, intent(in) :: a(d,d)	869	869	double complex, intent(in) :: a(d,d)
double complex, intent(out) :: evala(d)	870	870	double complex, intent(out) :: evala(d)
double complex, intent(out) :: eveca(d,d)	871	871	double complex, intent(out) :: eveca(d,d)
integer, intent(out) :: status	872	872	integer, intent(out) :: status
	873	873
double complex, allocatable :: wc_a(:,:) ! a working copy of a	874	874	double complex, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info	875	875	integer :: info
integer :: lwork	876	876	integer :: lwork
double complex, allocatable :: work(:)	877	877	double complex, allocatable :: work(:)
double complex :: twork(1)	878	878	double complex :: twork(1)
double precision, allocatable :: rwork(:)	879	879	double precision, allocatable :: rwork(:)
double complex :: vl ! unused but necessary for the call	880	880	double complex :: vl ! unused but necessary for the call
	881	881
status=1	882	882	status=1
	883	883
! making a working copy of a	884	884	! making a working copy of a
allocate(wc_a(d,d))	885	885	allocate(wc_a(d,d))
!call zcopy(d*d,a,1,wc_a,1)	886	886	!call zcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:)	887	887	wc_a(:,:)=a(:,:)
	888	888
! query optimal work size	889	889	! query optimal work size
call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)	890	890	call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
lwork=int(twork(1))	891	891	lwork=int(twork(1))
allocate(work(lwork))	892	892	allocate(work(lwork))
allocate(rwork(2*d))	893	893	allocate(rwork(2*d))
call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)	894	894	call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
	895	895
if (info /= 0) then	896	896	if (info /= 0) then
status=0	897	897	status=0
end if	898	898	end if
deallocate(rwork)	899	899	deallocate(rwork)
deallocate(work)	900	900	deallocate(work)
deallocate(wc_a)	901	901	deallocate(wc_a)
	902	902
end subroutine	903	903	end subroutine
	904	904
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	905	905	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	906	906	!
! Akima spline interpolation and spline evaluation	907	907	! Akima spline interpolation and spline evaluation
!	908	908	!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	909	909	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	910	910
! Single precision	911	911	! Single precision
subroutine fvn_s_akima(n,x,y,br,co)	912	912	subroutine fvn_s_akima(n,x,y,br,co)
implicit none	913	913	implicit none
integer, intent(in) :: n	914	914	integer, intent(in) :: n
real, intent(in) :: x(n)	915	915	real, intent(in) :: x(n)
real, intent(in) :: y(n)	916	916	real, intent(in) :: y(n)
real, intent(out) :: br(n)	917	917	real, intent(out) :: br(n)
real, intent(out) :: co(4,n)	918	918	real, intent(out) :: co(4,n)
	919	919
real, allocatable :: var(:),z(:)	920	920	real, allocatable :: var(:),z(:)
real :: wi_1,wi	921	921	real :: wi_1,wi
integer :: i	922	922	integer :: i
real :: dx,a,b	923	923	real :: dx,a,b
	924	924
! br is just a copy of x	925	925	! br is just a copy of x
br(:)=x(:)	926	926	br(:)=x(:)
	927	927
allocate(var(n))	928	928	allocate(var(n))
allocate(z(n))	929	929	allocate(z(n))
! evaluate the variations	930	930	! evaluate the variations
do i=1, n-1	931	931	do i=1, n-1
var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))	932	932	var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
end do	933	933	end do
var(n+2)=2.e0*var(n+1)-var(n)	934	934	var(n+2)=2.e0*var(n+1)-var(n)
var(n+3)=2.e0*var(n+2)-var(n+1)	935	935	var(n+3)=2.e0*var(n+2)-var(n+1)
var(2)=2.e0*var(3)-var(4)	936	936	var(2)=2.e0*var(3)-var(4)
var(1)=2.e0*var(2)-var(3)	937	937	var(1)=2.e0*var(2)-var(3)
	938	938
do i = 1, n	939	939	do i = 1, n
wi_1=abs(var(i+3)-var(i+2))	940	940	wi_1=abs(var(i+3)-var(i+2))
wi=abs(var(i+1)-var(i))	941	941	wi=abs(var(i+1)-var(i))
if ((wi_1+wi).eq.0.e0) then	942	942	if ((wi_1+wi).eq.0.e0) then
z(i)=(var(i+2)+var(i+1))/2.e0	943	943	z(i)=(var(i+2)+var(i+1))/2.e0
else	944	944	else
z(i)=(wi_1var(i+1)+wivar(i+2))/(wi_1+wi)	945	945	z(i)=(wi_1var(i+1)+wivar(i+2))/(wi_1+wi)
end if	946	946	end if
end do	947	947	end do
	948	948
do i=1, n-1	949	949	do i=1, n-1
dx=x(i+1)-x(i)	950	950	dx=x(i+1)-x(i)
a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd	951	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	952	952	b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
co(1,i)=y(i)	953	953	co(1,i)=y(i)
co(2,i)=z(i)	954	954	co(2,i)=z(i)
!co(3,i)=-(a-3.b)/dx*2 ! méthode wd	955	955	!co(3,i)=-(a-3.b)/dx*2 ! méthode wd
!co(4,i)=(a-2.b)/dx*3 ! méthode wd	956	956	!co(4,i)=(a-2.b)/dx*3 ! méthode wd
co(3,i)=(3.e0var(i+2)-2.e0z(i)-z(i+1))/dx ! méthode JP Moreau	957	957	co(3,i)=(3.e0var(i+2)-2.e0z(i)-z(i+1))/dx ! méthode JP Moreau
co(4,i)=(z(i)+z(i+1)-2.e0var(i+2))/dx*2 !	958	958	co(4,i)=(z(i)+z(i+1)-2.e0var(i+2))/dx*2 !
! les coefficients donnés par imsl sont co(3,i)2 et co(4,i)6	959	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 ...	960	960	! etrangement la fonction csval corrige et donne la bonne valeur ...
end do	961	961	end do
co(1,n)=y(n)	962	962	co(1,n)=y(n)
co(2,n)=z(n)	963	963	co(2,n)=z(n)
co(3,n)=0.e0	964	964	co(3,n)=0.e0
co(4,n)=0.e0	965	965	co(4,n)=0.e0
	966	966
deallocate(z)	967	967	deallocate(z)
deallocate(var)	968	968	deallocate(var)
	969	969
end subroutine	970	970	end subroutine
	971	971
! Double precision	972	972	! Double precision
subroutine fvn_d_akima(n,x,y,br,co)	973	973	subroutine fvn_d_akima(n,x,y,br,co)
	974	974
implicit none	975	975	implicit none
integer, intent(in) :: n	976	976	integer, intent(in) :: n
double precision, intent(in) :: x(n)	977	977	double precision, intent(in) :: x(n)
double precision, intent(in) :: y(n)	978	978	double precision, intent(in) :: y(n)
double precision, intent(out) :: br(n)	979	979	double precision, intent(out) :: br(n)
double precision, intent(out) :: co(4,n)	980	980	double precision, intent(out) :: co(4,n)
	981	981
double precision, allocatable :: var(:),z(:)	982	982	double precision, allocatable :: var(:),z(:)
double precision :: wi_1,wi	983	983	double precision :: wi_1,wi
integer :: i	984	984	integer :: i
double precision :: dx,a,b	985	985	double precision :: dx,a,b
	986	986
! br is just a copy of x	987	987	! br is just a copy of x
br(:)=x(:)	988	988	br(:)=x(:)
	989	989
allocate(var(n))	990	990	allocate(var(n))
allocate(z(n))	991	991	allocate(z(n))
! evaluate the variations	992	992	! evaluate the variations
do i=1, n-1	993	993	do i=1, n-1
var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))	994	994	var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
end do	995	995	end do
var(n+2)=2.d0*var(n+1)-var(n)	996	996	var(n+2)=2.d0*var(n+1)-var(n)
var(n+3)=2.d0*var(n+2)-var(n+1)	997	997	var(n+3)=2.d0*var(n+2)-var(n+1)
var(2)=2.d0*var(3)-var(4)	998	998	var(2)=2.d0*var(3)-var(4)
var(1)=2.d0*var(2)-var(3)	999	999	var(1)=2.d0*var(2)-var(3)
	1000	1000
do i = 1, n	1001	1001	do i = 1, n
wi_1=dabs(var(i+3)-var(i+2))	1002	1002	wi_1=dabs(var(i+3)-var(i+2))
wi=dabs(var(i+1)-var(i))	1003	1003	wi=dabs(var(i+1)-var(i))
if ((wi_1+wi).eq.0.d0) then	1004	1004	if ((wi_1+wi).eq.0.d0) then
z(i)=(var(i+2)+var(i+1))/2.d0	1005	1005	z(i)=(var(i+2)+var(i+1))/2.d0
else	1006	1006	else
z(i)=(wi_1var(i+1)+wivar(i+2))/(wi_1+wi)	1007	1007	z(i)=(wi_1var(i+1)+wivar(i+2))/(wi_1+wi)
end if	1008	1008	end if
end do	1009	1009	end do
	1010	1010
do i=1, n-1	1011	1011	do i=1, n-1
dx=x(i+1)-x(i)	1012	1012	dx=x(i+1)-x(i)
a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd	1013	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	1014	1014	b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
co(1,i)=y(i)	1015	1015	co(1,i)=y(i)
co(2,i)=z(i)	1016	1016	co(2,i)=z(i)
!co(3,i)=-(a-3.b)/dx*2 ! méthode wd	1017	1017	!co(3,i)=-(a-3.b)/dx*2 ! méthode wd
!co(4,i)=(a-2.b)/dx*3 ! méthode wd	1018	1018	!co(4,i)=(a-2.b)/dx*3 ! méthode wd
co(3,i)=(3.d0var(i+2)-2.d0z(i)-z(i+1))/dx ! méthode JP Moreau	1019	1019	co(3,i)=(3.d0var(i+2)-2.d0z(i)-z(i+1))/dx ! méthode JP Moreau
co(4,i)=(z(i)+z(i+1)-2.d0var(i+2))/dx*2 !	1020	1020	co(4,i)=(z(i)+z(i+1)-2.d0var(i+2))/dx*2 !
! les coefficients donnés par imsl sont co(3,i)2 et co(4,i)6	1021	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 ...	1022	1022	! etrangement la fonction csval corrige et donne la bonne valeur ...
end do	1023	1023	end do
co(1,n)=y(n)	1024	1024	co(1,n)=y(n)
co(2,n)=z(n)	1025	1025	co(2,n)=z(n)
co(3,n)=0.d0	1026	1026	co(3,n)=0.d0
co(4,n)=0.d0	1027	1027	co(4,n)=0.d0
	1028	1028
deallocate(z)	1029	1029	deallocate(z)
deallocate(var)	1030	1030	deallocate(var)
	1031	1031
end subroutine	1032	1032	end subroutine
	1033	1033
!	1034	1034	!
! Single precision spline evaluation	1035	1035	! Single precision spline evaluation
!	1036	1036	!
function fvn_s_spline_eval(x,n,br,co)	1037	1037	function fvn_s_spline_eval(x,n,br,co)
implicit none	1038	1038	implicit none
real, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated	1039	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	1040	1040	integer, intent(in) :: n ! number of intervals
real, intent(in) :: br(n+1) ! breakpoints	1041	1041	real, intent(in) :: br(n+1) ! breakpoints
real, intent(in) :: co(4,n+1) ! spline coeeficients	1042	1042	real, intent(in) :: co(4,n+1) ! spline coeeficients
real :: fvn_s_spline_eval	1043	1043	real :: fvn_s_spline_eval
	1044	1044
integer :: i	1045	1045	integer :: i
real :: dx	1046	1046	real :: dx
	1047	1047
if (x<=br(1)) then	1048	1048	if (x<=br(1)) then
i=1	1049	1049	i=1
else if (x>=br(n+1)) then	1050	1050	else if (x>=br(n+1)) then
i=n	1051	1051	i=n
else	1052	1052	else
i=1	1053	1053	i=1
do while(x>=br(i))	1054	1054	do while(x>=br(i))
i=i+1	1055	1055	i=i+1
end do	1056	1056	end do
i=i-1	1057	1057	i=i-1
end if	1058	1058	end if
dx=x-br(i)	1059	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	1060	1060	fvn_s_spline_eval=co(1,i)+co(2,i)dx+co(3,i)dx*2+co(4,i)dx**3
	1061	1061
end function	1062	1062	end function
	1063	1063
! Double precision spline evaluation	1064	1064	! Double precision spline evaluation
function fvn_d_spline_eval(x,n,br,co)	1065	1065	function fvn_d_spline_eval(x,n,br,co)
implicit none	1066	1066	implicit none
double precision, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated	1067	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	1068	1068	integer, intent(in) :: n ! number of intervals
double precision, intent(in) :: br(n+1) ! breakpoints	1069	1069	double precision, intent(in) :: br(n+1) ! breakpoints
double precision, intent(in) :: co(4,n+1) ! spline coeeficients	1070	1070	double precision, intent(in) :: co(4,n+1) ! spline coeeficients
double precision :: fvn_d_spline_eval	1071	1071	double precision :: fvn_d_spline_eval
	1072	1072
integer :: i	1073	1073	integer :: i
double precision :: dx	1074	1074	double precision :: dx
	1075	1075
	1076	1076
if (x<=br(1)) then	1077	1077	if (x<=br(1)) then
i=1	1078	1078	i=1
else if (x>=br(n+1)) then	1079	1079	else if (x>=br(n+1)) then
i=n	1080	1080	i=n
else	1081	1081	else
i=1	1082	1082	i=1
do while(x>=br(i))	1083	1083	do while(x>=br(i))
i=i+1	1084	1084	i=i+1
end do	1085	1085	end do
i=i-1	1086	1086	i=i-1
end if	1087	1087	end if
	1088	1088
dx=x-br(i)	1089	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	1090	1090	fvn_d_spline_eval=co(1,i)+co(2,i)dx+co(3,i)dx*2+co(4,i)dx**3
	1091	1091
end function	1092	1092	end function
	1093	1093
	1094	1094
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	1095	1095	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	1096	1096	!
! Least square problem	1097	1097	! Least square problem
!	1098	1098	!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	1099	1099	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	1100	1100	!
!	1101	1101	!
subroutine fvn_d_lspoly(np,x,y,deg,coeff,status)	1102	1102	subroutine fvn_d_lspoly(np,x,y,deg,coeff,status)
!	1103	1103	!
! Least square polynomial fitting	1104	1104	! Least square polynomial fitting
!	1105	1105	!
! Find the coefficients of the least square polynomial of a given degree	1106	1106	! Find the coefficients of the least square polynomial of a given degree
! for a set of coordinates.	1107	1107	! for a set of coordinates.
!	1108	1108	!
! The degree must be lower than the number of points	1109	1109	! The degree must be lower than the number of points
!	1110	1110	!
! np (in) : number of points	1111	1111	! np (in) : number of points
! x(np) (in) : x data	1112	1112	! x(np) (in) : x data
! y(np) (in) : y data	1113	1113	! y(np) (in) : y data
! deg (in) : polynomial's degree	1114	1114	! deg (in) : polynomial's degree
! coeff(deg+1) (out) : polynomial coefficients	1115	1115	! coeff(deg+1) (out) : polynomial coefficients
! status (out) : =0 if a problem occurs	1116	1116	! status (out) : =0 if a problem occurs
!	1117	1117	!
implicit none	1118	1118	implicit none
	1119	1119
integer, intent(in) :: np,deg	1120	1120	integer, intent(in) :: np,deg
real(kind=8), intent(in), dimension(np) :: x,y	1121	1121	real(kind=8), intent(in), dimension(np) :: x,y
real(kind=8), intent(out), dimension(deg+1) :: coeff	1122	1122	real(kind=8), intent(out), dimension(deg+1) :: coeff
integer, intent(out) :: status	1123	1123	integer, intent(out) :: status
	1124	1124
real(kind=8), allocatable, dimension(:,:) :: mat,bmat	1125	1125	real(kind=8), allocatable, dimension(:,:) :: mat,bmat
real(kind=8),dimension(:),allocatable :: work,singval	1126	1126	real(kind=8),dimension(:),allocatable :: work,singval
real(kind=8),dimension(1) :: twork	1127	1127	real(kind=8),dimension(1) :: twork
integer :: lwork,info,rank	1128	1128	integer :: lwork,info,rank
	1129	1129
integer :: i,j	1130	1130	integer :: i,j
	1131	1131
status=1	1132	1132	status=1
allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))	1133	1133	allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))
	1134	1134
! Design matrix valorisation	1135	1135	! Design matrix valorisation
mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )	1136	1136	mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
	1137	1137
! second member valorisation	1138	1138	! second member valorisation
bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))	1139	1139	bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))
	1140	1140
! query workspace size	1141	1141	! query workspace size
call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)	1142	1142	call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)
lwork=twork(1)	1143	1143	lwork=twork(1)
allocate(work(int(lwork)))	1144	1144	allocate(work(int(lwork)))
! real call	1145	1145	! real call
call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)	1146	1146	call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)
	1147	1147
if (info /= 0) then	1148	1148	if (info /= 0) then
status=0	1149	1149	status=0
end if	1150	1150	end if
	1151	1151
coeff = (/ (bmat(i,1),i=1,deg+1) /)	1152	1152	coeff = (/ (bmat(i,1),i=1,deg+1) /)
	1153	1153
deallocate(work)	1154	1154	deallocate(work)
deallocate(mat,bmat,singval)	1155	1155	deallocate(mat,bmat,singval)
end subroutine	1156	1156	end subroutine
	1157	1157
subroutine fvn_s_lspoly(np,x,y,deg,coeff,status)	1158	1158	subroutine fvn_s_lspoly(np,x,y,deg,coeff,status)
!	1159	1159	!
! Least square polynomial fitting	1160	1160	! Least square polynomial fitting
!	1161	1161	!
! Find the coefficients of the least square polynomial of a given degree	1162	1162	! Find the coefficients of the least square polynomial of a given degree
! for a set of coordinates.	1163	1163	! for a set of coordinates.
!	1164	1164	!
! The degree must be lower than the number of points	1165	1165	! The degree must be lower than the number of points
!	1166	1166	!
! np (in) : number of points	1167	1167	! np (in) : number of points
! x(np) (in) : x data	1168	1168	! x(np) (in) : x data
! y(np) (in) : y data	1169	1169	! y(np) (in) : y data
! deg (in) : polynomial's degree	1170	1170	! deg (in) : polynomial's degree
! coeff(deg+1) (out) : polynomial coefficients	1171	1171	! coeff(deg+1) (out) : polynomial coefficients
! status (out) : =0 if a problem occurs	1172	1172	! status (out) : =0 if a problem occurs
!	1173	1173	!
implicit none	1174	1174	implicit none
	1175	1175
integer, intent(in) :: np,deg	1176	1176	integer, intent(in) :: np,deg
real(kind=4), intent(in), dimension(np) :: x,y	1177	1177	real(kind=4), intent(in), dimension(np) :: x,y
real(kind=4), intent(out), dimension(deg+1) :: coeff	1178	1178	real(kind=4), intent(out), dimension(deg+1) :: coeff
integer, intent(out) :: status	1179	1179	integer, intent(out) :: status
	1180	1180
real(kind=4), allocatable, dimension(:,:) :: mat,bmat	1181	1181	real(kind=4), allocatable, dimension(:,:) :: mat,bmat
real(kind=4),dimension(:),allocatable :: work,singval	1182	1182	real(kind=4),dimension(:),allocatable :: work,singval
real(kind=4),dimension(1) :: twork	1183	1183	real(kind=4),dimension(1) :: twork
integer :: lwork,info,rank	1184	1184	integer :: lwork,info,rank
	1185	1185
integer :: i,j	1186	1186	integer :: i,j
	1187	1187
status=1	1188	1188	status=1
allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))	1189	1189	allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))
	1190	1190
! Design matrix valorisation	1191	1191	! Design matrix valorisation
mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )	1192	1192	mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
	1193	1193
! second member valorisation	1194	1194	! second member valorisation
bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))	1195	1195	bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))
	1196	1196
! query workspace size	1197	1197	! query workspace size
call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)	1198	1198	call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)
lwork=twork(1)	1199	1199	lwork=twork(1)
allocate(work(int(lwork)))	1200	1200	allocate(work(int(lwork)))
! real call	1201	1201	! real call
call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)	1202	1202	call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)
	1203	1203
if (info /= 0) then	1204	1204	if (info /= 0) then
status=0	1205	1205	status=0
end if	1206	1206	end if
	1207	1207
coeff = (/ (bmat(i,1),i=1,deg+1) /)	1208	1208	coeff = (/ (bmat(i,1),i=1,deg+1) /)
	1209	1209
deallocate(work)	1210	1210	deallocate(work)
deallocate(mat,bmat,singval)	1211	1211	deallocate(mat,bmat,singval)
end subroutine	1212	1212	end subroutine
	1213	1213
	1214	1214
!	1215	1215	!
! Muller	1216	1216	! Muller
!	1217	1217	!
!	1218	1218	!
!	1219	1219	!
! William Daniau 2007	1220	1220	! William Daniau 2007
!	1221	1221	!
! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f	1222	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	1223	1223	! http://plato.asu.edu/ftp/other_software/muller.f
!	1224	1224	!
! it can be used as a replacement for imsl routine dzanly with minor changes	1225	1225	! it can be used as a replacement for imsl routine dzanly with minor changes
!	1226	1226	!
!-----------------------------------------------------------------------	1227	1227	!-----------------------------------------------------------------------
!	1228	1228	!
! purpose - zeros of an analytic complex function	1229	1229	! purpose - zeros of an analytic complex function
! using the muller method with deflation	1230	1230	! using the muller method with deflation
!	1231	1231	!
! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,	1232	1232	! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,
! infer,ier)	1233	1233	! infer,ier)
!	1234	1234	!
! arguments f - a complex function subprogram, f(z), written	1235	1235	! arguments f - a complex function subprogram, f(z), written
! by the user specifying the equation whose	1236	1236	! by the user specifying the equation whose
! roots are to be found. f must appear in	1237	1237	! roots are to be found. f must appear in
! an external statement in the calling pro-	1238	1238	! an external statement in the calling pro-
! gram.	1239	1239	! gram.
! eps - 1st stopping criterion. let fp(z)=f(z)/p	1240	1240	! eps - 1st stopping criterion. let fp(z)=f(z)/p
! where p = (z-z(1))(z-z(2)),,,*(z-z(k-1))	1241	1241	! where p = (z-z(1))(z-z(2)),,,*(z-z(k-1))
! and z(1),...,z(k-1) are previously found	1242	1242	! and z(1),...,z(k-1) are previously found
! roots. if ((cdabs(f(z)).le.eps) .and.	1243	1243	! roots. if ((cdabs(f(z)).le.eps) .and.
! (cdabs(fp(z)).le.eps)), then z is accepted	1244	1244	! (cdabs(fp(z)).le.eps)), then z is accepted
! as a root. (input)	1245	1245	! as a root. (input)
! eps1 - 2nd stopping criterion. a root is accepted	1246	1246	! eps1 - 2nd stopping criterion. a root is accepted
! if two successive approximations to a given	1247	1247	! if two successive approximations to a given
! root agree within eps1. (input)	1248	1248	! root agree within eps1. (input)
! note. if either or both of the stopping	1249	1249	! note. if either or both of the stopping
! criteria are fulfilled, the root is	1250	1250	! criteria are fulfilled, the root is
! accepted.	1251	1251	! accepted.
! kn - the number of known roots which must be stored	1252	1252	! kn - the number of known roots which must be stored
! in x(1),...,x(kn), prior to entry to muller	1253	1253	! in x(1),...,x(kn), prior to entry to muller
! nguess - the number of initial guesses provided. these	1254	1254	! nguess - the number of initial guesses provided. these
! guesses must be stored in x(kn+1),...,	1255	1255	! guesses must be stored in x(kn+1),...,
! x(kn+nguess). nguess must be set equal	1256	1256	! x(kn+nguess). nguess must be set equal
! to zero if no guesses are provided. (input)	1257	1257	! to zero if no guesses are provided. (input)
! n - the number of new roots to be found by	1258	1258	! n - the number of new roots to be found by
! muller (input)	1259	1259	! muller (input)
! x - a complex vector of length kn+n. x(1),...,	1260	1260	! x - a complex vector of length kn+n. x(1),...,
! x(kn) on input must contain any known	1261	1261	! x(kn) on input must contain any known
! roots. x(kn+1),..., x(kn+n) on input may,	1262	1262	! roots. x(kn+1),..., x(kn+n) on input may,
! on user option, contain initial guesses for	1263	1263	! on user option, contain initial guesses for
! the n new roots which are to be computed.	1264	1264	! the n new roots which are to be computed.
! if the user does not provide an initial	1265	1265	! if the user does not provide an initial
! guess, zero is used.	1266	1266	! guess, zero is used.
! on output, x(kn+1),...,x(kn+n) contain the	1267	1267	! on output, x(kn+1),...,x(kn+n) contain the
! approximate roots found by muller.	1268	1268	! approximate roots found by muller.
! itmax - the maximum allowable number of iterations	1269	1269	! itmax - the maximum allowable number of iterations
! per root (input)	1270	1270	! per root (input)
! infer - an integer vector of length kn+n. on	1271	1271	! infer - an integer vector of length kn+n. on
! output infer(j) contains the number of	1272	1272	! output infer(j) contains the number of
! iterations used in finding the j-th root	1273	1273	! iterations used in finding the j-th root
! when convergence was achieved. if	1274	1274	! when convergence was achieved. if
! convergence was not obtained in itmax	1275	1275	! convergence was not obtained in itmax
! iterations, infer(j) will be greater than	1276	1276	! iterations, infer(j) will be greater than
! itmax (output).	1277	1277	! itmax (output).
! ier - error parameter (output)	1278	1278	! ier - error parameter (output)
! warning error	1279	1279	! warning error
! ier = 33 indicates failure to converge with-	1280	1280	! ier = 33 indicates failure to converge with-
! in itmax iterations for at least one of	1281	1281	! in itmax iterations for at least one of
! the (n) new roots.	1282	1282	! the (n) new roots.
!	1283	1283	!
!	1284	1284	!
! remarks muller always returns the last approximation for root j	1285	1285	! remarks muller always returns the last approximation for root j
! in x(j). if the convergence criterion is satisfied,	1286	1286	! in x(j). if the convergence criterion is satisfied,
! then infer(j) is less than or equal to itmax. if the	1287	1287	! then infer(j) is less than or equal to itmax. if the
! convergence criterion is not satisified, then infer(j)	1288	1288	! convergence criterion is not satisified, then infer(j)
! is set to either itmax+1 or itmax+k, with k greater	1289	1289	! is set to either itmax+1 or itmax+k, with k greater
! than 1. infer(j) = itmax+1 indicates that muller did	1290	1290	! than 1. infer(j) = itmax+1 indicates that muller did
! not obtain convergence in the allowed number of iter-	1291	1291	! not obtain convergence in the allowed number of iter-
! ations. in this case, the user may wish to set itmax	1292	1292	! ations. in this case, the user may wish to set itmax
! to a larger value. infer(j) = itmax+k means that con-	1293	1293	! to a larger value. infer(j) = itmax+k means that con-
! vergence was obtained (on iteration k) for the defla-	1294	1294	! vergence was obtained (on iteration k) for the defla-
! ted function	1295	1295	! ted function
! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))	1296	1296	! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))
!	1297	1297	!
! but failed for f(z). in this case, better initial	1298	1298	! but failed for f(z). in this case, better initial
! guesses might help or, it might be necessary to relax	1299	1299	! guesses might help or, it might be necessary to relax
! the convergence criterion.	1300	1300	! the convergence criterion.
!	1301	1301	!
!-----------------------------------------------------------------------	1302	1302	!-----------------------------------------------------------------------
!	1303	1303	!
subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)	1304	1304	subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
implicit none	1305	1305	implicit none
double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq	1306	1306	double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq
double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &	1307	1307	double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &	1308	1308	tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
zero,p1,one,four,p5	1309	1309	zero,p1,one,four,p5
	1310	1310
double complex, external :: f	1311	1311	double complex, external :: f
integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &	1312	1312	integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
knpng,jk,ick,nn,lm1,errcode	1313	1313	knpng,jk,ick,nn,lm1,errcode
double complex :: x(kn+n)	1314	1314	double complex :: x(kn+n)
integer :: infer(kn+n)	1315	1315	integer :: infer(kn+n)
	1316	1316
	1317	1317
data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &	1318	1318	data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &
one/(1.0,0.0)/,four/(4.0,0.0)/, &	1319	1319	one/(1.0,0.0)/,four/(4.0,0.0)/, &
p5/(0.5,0.0)/, &	1320	1320	p5/(0.5,0.0)/, &
rzero/0.0/,rten/10.0/,rhun/100.0/, &	1321	1321	rzero/0.0/,rten/10.0/,rhun/100.0/, &
ax/0.1/,ickmax/3/,rp01/0.01/	1322	1322	ax/0.1/,ickmax/3/,rp01/0.01/
	1323	1323
ier = 0	1324	1324	ier = 0
if (n .lt. 1) then ! What the hell are doing here then ...	1325	1325	if (n .lt. 1) then ! What the hell are doing here then ...
return	1326	1326	return
end if	1327	1327	end if
!eps1 = rten **(-nsig)	1328	1328	!eps1 = rten **(-nsig)
eps1 = min(eps1,rp01)	1329	1329	eps1 = min(eps1,rp01)
	1330	1330
knp1 = kn+1	1331	1331	knp1 = kn+1
knpn = kn+n	1332	1332	knpn = kn+n
knpng = kn+nguess	1333	1333	knpng = kn+nguess
do i=1,knpn	1334	1334	do i=1,knpn
infer(i) = 0	1335	1335	infer(i) = 0
if (i .gt. knpng) x(i) = zero	1336	1336	if (i .gt. knpng) x(i) = zero
end do	1337	1337	end do
l= knp1	1338	1338	l= knp1
	1339	1339
ic=0	1340	1340	ic=0
rloop: do while (l<=knpn) ! Main loop over new roots	1341	1341	rloop: do while (l<=knpn) ! Main loop over new roots
jk = 0	1342	1342	jk = 0
ick = 0	1343	1343	ick = 0
xl = x(l)	1344	1344	xl = x(l)
icloop: do	1345	1345	icloop: do
ic = 0	1346	1346	ic = 0
h = ax	1347	1347	h = ax
h = p1*h	1348	1348	h = p1*h
if (cdabs(xl) .gt. ax) h = p1*xl	1349	1349	if (cdabs(xl) .gt. ax) h = p1*xl
! first three points are	1350	1350	! first three points are
! xl+h, xl-h, xl	1351	1351	! xl+h, xl-h, xl
rt = xl+h	1352	1352	rt = xl+h
call deflated_work(errcode)	1353	1353	call deflated_work(errcode)
if (errcode == 1) then	1354	1354	if (errcode == 1) then
exit icloop	1355	1355	exit icloop
end if	1356	1356	end if
	1357	1357
z0 = fprt	1358	1358	z0 = fprt
y0 = frt	1359	1359	y0 = frt
x0 = rt	1360	1360	x0 = rt
rt = xl-h	1361	1361	rt = xl-h
call deflated_work(errcode)	1362	1362	call deflated_work(errcode)
if (errcode == 1) then	1363	1363	if (errcode == 1) then
exit icloop	1364	1364	exit icloop
end if	1365	1365	end if
	1366	1366
z1 = fprt	1367	1367	z1 = fprt
y1 = frt	1368	1368	y1 = frt
h = xl-rt	1369	1369	h = xl-rt
d = h/(rt-x0)	1370	1370	d = h/(rt-x0)
rt = xl	1371	1371	rt = xl
	1372	1372
call deflated_work(errcode)	1373	1373	call deflated_work(errcode)
if (errcode == 1) then	1374	1374	if (errcode == 1) then
exit icloop	1375	1375	exit icloop
end if	1376	1376	end if
	1377	1377
	1378	1378
z2 = fprt	1379	1379	z2 = fprt
y2 = frt	1380	1380	y2 = frt
! begin main algorithm	1381	1381	! begin main algorithm
iloop: do	1382	1382	iloop: do
dd = one + d	1383	1383	dd = one + d
t1 = z0dd	1384	1384	t1 = z0dd
t2 = z1dddd	1385	1385	t2 = z1dddd
xx = z2*dd	1386	1386	xx = z2*dd
t3 = z2*d	1387	1387	t3 = z2*d
bi = t1-t2+xx+t3	1388	1388	bi = t1-t2+xx+t3
den = bibi-four(xxt1-t3(t2-xx))	1389	1389	den = bibi-four(xxt1-t3(t2-xx))
! use denominator of maximum amplitude	1390	1390	! use denominator of maximum amplitude
t1 = cdsqrt(den)	1391	1391	t1 = cdsqrt(den)
qz = rhun*max(cdabs(bi),cdabs(t1))	1392	1392	qz = rhun*max(cdabs(bi),cdabs(t1))
t2 = bi + t1	1393	1393	t2 = bi + t1
tpq = cdabs(t2)+qz	1394	1394	tpq = cdabs(t2)+qz
if (tpq .eq. qz) t2 = zero	1395	1395	if (tpq .eq. qz) t2 = zero
t3 = bi - t1	1396	1396	t3 = bi - t1
tpq = cdabs(t3) + qz	1397	1397	tpq = cdabs(t3) + qz
if (tpq .eq. qz) t3 = zero	1398	1398	if (tpq .eq. qz) t3 = zero
den = t2	1399	1399	den = t2
qz = cdabs(t3)-cdabs(t2)	1400	1400	qz = cdabs(t3)-cdabs(t2)
if (qz .gt. rzero) den = t3	1401	1401	if (qz .gt. rzero) den = t3
! test for zero denominator	1402	1402	! test for zero denominator
if (cdabs(den) .eq. rzero) then	1403	1403	if (cdabs(den) .eq. rzero) then
call trans_rt()	1404	1404	call trans_rt()
call deflated_work(errcode)	1405	1405	call deflated_work(errcode)
if (errcode == 1) then	1406	1406	if (errcode == 1) then
exit icloop	1407	1407	exit icloop
end if	1408	1408	end if
z2 = fprt	1409	1409	z2 = fprt
y2 = frt	1410	1410	y2 = frt
cycle iloop	1411	1411	cycle iloop
end if	1412	1412	end if
	1413	1413
	1414	1414
d = -xx/den	1415	1415	d = -xx/den
d = d+d	1416	1416	d = d+d
h = d*h	1417	1417	h = d*h
rt = rt + h	1418	1418	rt = rt + h
! check convergence of the first kind	1419	1419	! check convergence of the first kind
if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then	1420	1420	if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then
if (ic .ne. 0) then	1421	1421	if (ic .ne. 0) then
exit icloop	1422	1422	exit icloop
end if	1423	1423	end if
ic = 1	1424	1424	ic = 1
z0 = y1	1425	1425	z0 = y1
z1 = y2	1426	1426	z1 = y2
z2 = f(rt)	1427	1427	z2 = f(rt)
xl = rt	1428	1428	xl = rt
ick = ick+1	1429	1429	ick = ick+1
if (ick .le. ickmax) then	1430	1430	if (ick .le. ickmax) then
cycle iloop	1431	1431	cycle iloop
end if	1432	1432	end if
! warning error, itmax = maximum	1433	1433	! warning error, itmax = maximum
jk = itmax + jk	1434	1434	jk = itmax + jk
ier = 33	1435	1435	ier = 33
end if	1436	1436	end if
if (ic .ne. 0) then	1437	1437	if (ic .ne. 0) then
cycle icloop	1438	1438	cycle icloop
end if	1439	1439	end if
call deflated_work(errcode)	1440	1440	call deflated_work(errcode)
if (errcode == 1) then	1441	1441	if (errcode == 1) then
exit icloop	1442	1442	exit icloop
end if	1443	1443	end if
	1444	1444
do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)	1445	1445	do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)
! take remedial action to induce	1446	1446	! take remedial action to induce
! convergence	1447	1447	! convergence
d = d*p5	1448	1448	d = d*p5
h = h*p5	1449	1449	h = h*p5
rt = rt-h	1450	1450	rt = rt-h
call deflated_work(errcode)	1451	1451	call deflated_work(errcode)
if (errcode == 1) then	1452	1452	if (errcode == 1) then
exit icloop	1453	1453	exit icloop
end if	1454	1454	end if
end do	1455	1455	end do
z0 = z1	1456	1456	z0 = z1
z1 = z2	1457	1457	z1 = z2
z2 = fprt	1458	1458	z2 = fprt
y0 = y1	1459	1459	y0 = y1
y1 = y2	1460	1460	y1 = y2
y2 = frt	1461	1461	y2 = frt
end do iloop	1462	1462	end do iloop
end do icloop	1463	1463	end do icloop
x(l) = rt	1464	1464	x(l) = rt
infer(l) = jk	1465	1465	infer(l) = jk
l = l+1	1466	1466	l = l+1
end do rloop	1467	1467	end do rloop
	1468	1468
contains	1469	1469	contains
subroutine trans_rt()	1470	1470	subroutine trans_rt()
tem = rten*eps1	1471	1471	tem = rten*eps1
if (cdabs(rt) .gt. ax) tem = tem*rt	1472	1472	if (cdabs(rt) .gt. ax) tem = tem*rt
rt = rt+tem	1473	1473	rt = rt+tem
d = (h+tem)*d/h	1474	1474	d = (h+tem)*d/h
h = h+tem	1475	1475	h = h+tem
end subroutine trans_rt	1476	1476	end subroutine trans_rt
	1477	1477
subroutine deflated_work(errcode)	1478	1478	subroutine deflated_work(errcode)
! errcode=0 => no errors	1479	1479	! errcode=0 => no errors
! errcode=1 => jk>itmax or convergence of second kind achieved	1480	1480	! errcode=1 => jk>itmax or convergence of second kind achieved
integer :: errcode,flag	1481	1481	integer :: errcode,flag
	1482	1482
flag=1	1483	1483	flag=1
loop1: do while(flag==1)	1484	1484	loop1: do while(flag==1)
errcode=0	1485	1485	errcode=0
jk = jk+1	1486	1486	jk = jk+1
if (jk .gt. itmax) then	1487	1487	if (jk .gt. itmax) then
ier=33	1488	1488	ier=33
errcode=1	1489	1489	errcode=1
return	1490	1490	return
end if	1491	1491	end if
frt = f(rt)	1492	1492	frt = f(rt)
fprt = frt	1493	1493	fprt = frt
if (l /= 1) then	1494	1494	if (l /= 1) then
lm1 = l-1	1495	1495	lm1 = l-1
do i=1,lm1	1496	1496	do i=1,lm1
tem = rt - x(i)	1497	1497	tem = rt - x(i)
if (cdabs(tem) .eq. rzero) then	1498	1498	if (cdabs(tem) .eq. rzero) then
!if (ic .ne. 0) go to 15 !! ?? possible?	1499	1499	!if (ic .ne. 0) go to 15 !! ?? possible?
call trans_rt()	1500	1500	call trans_rt()
cycle loop1	1501	1501	cycle loop1
end if	1502	1502	end if
fprt = fprt/tem	1503	1503	fprt = fprt/tem
end do	1504	1504	end do
end if	1505	1505	end if
flag=0	1506	1506	flag=0
end do loop1	1507	1507	end do loop1
	1508	1508
if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then	1509	1509	if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then
errcode=1	1510	1510	errcode=1
return	1511	1511	return
end if	1512	1512	end if
	1513	1513
end subroutine deflated_work	1514	1514	end subroutine deflated_work
	1515	1515
end subroutine	1516	1516	end subroutine
	1517	1517
	1518	1518
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	1519	1519	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	1520	1520	!
! Integration	1521	1521	! Integration
!	1522	1522	!
! Only double precision coded atm	1523	1523	! Only double precision coded atm
!	1524	1524	!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	1525	1525	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1526	1526
	1527	1527
subroutine fvn_d_gauss_legendre(n,qx,qw)	1528	1528	subroutine fvn_d_gauss_legendre(n,qx,qw)
!	1529	1529	!
! This routine compute the n Gauss Legendre abscissas and weights	1530	1530	! This routine compute the n Gauss Legendre abscissas and weights
! Adapted from Numerical Recipes routine gauleg	1531	1531	! Adapted from Numerical Recipes routine gauleg
!	1532	1532	!
! n (in) : number of points	1533	1533	! n (in) : number of points
! qx(out) : abscissas	1534	1534	! qx(out) : abscissas
! qw(out) : weights	1535	1535	! qw(out) : weights
!	1536	1536	!
implicit none	1537	1537	implicit none
double precision,parameter :: pi=3.141592653589793d0	1538	1538	double precision,parameter :: pi=3.141592653589793d0
integer, intent(in) :: n	1539	1539	integer, intent(in) :: n
double precision, intent(out) :: qx(n),qw(n)	1540	1540	double precision, intent(out) :: qx(n),qw(n)
	1541	1541
integer :: m,i,j	1542	1542	integer :: m,i,j
double precision :: z,z1,p1,p2,p3,pp	1543	1543	double precision :: z,z1,p1,p2,p3,pp
m=(n+1)/2	1544	1544	m=(n+1)/2
do i=1,m	1545	1545	do i=1,m
z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))	1546	1546	z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))
iloop: do	1547	1547	iloop: do
p1=1.d0	1548	1548	p1=1.d0
p2=0.d0	1549	1549	p2=0.d0
do j=1,n	1550	1550	do j=1,n
p3=p2	1551	1551	p3=p2
p2=p1	1552	1552	p2=p1
p1=((2.d0dble(j)-1.d0)zp2-(dble(j)-1.d0)p3)/dble(j)	1553	1553	p1=((2.d0dble(j)-1.d0)zp2-(dble(j)-1.d0)p3)/dble(j)
end do	1554	1554	end do
pp=dble(n)(zp1-p2)/(z*z-1.d0)	1555	1555	pp=dble(n)(zp1-p2)/(z*z-1.d0)
z1=z	1556	1556	z1=z
z=z1-p1/pp	1557	1557	z=z1-p1/pp
if (dabs(z-z1)<=epsilon(z)) then	1558	1558	if (dabs(z-z1)<=epsilon(z)) then
exit iloop	1559	1559	exit iloop
end if	1560	1560	end if
end do iloop	1561	1561	end do iloop
qx(i)=-z	1562	1562	qx(i)=-z
qx(n+1-i)=z	1563	1563	qx(n+1-i)=z
qw(i)=2.d0/((1.d0-zz)pp*pp)	1564	1564	qw(i)=2.d0/((1.d0-zz)pp*pp)
qw(n+1-i)=qw(i)	1565	1565	qw(n+1-i)=qw(i)
end do	1566	1566	end do
end subroutine	1567	1567	end subroutine
	1568	1568
	1569	1569
	1570	1570
subroutine fvn_d_gl_integ(f,a,b,n,res)	1571	1571	subroutine fvn_d_gl_integ(f,a,b,n,res)
!	1572	1572	!
! This is a simple non adaptative integration routine	1573	1573	! This is a simple non adaptative integration routine
! using n gauss legendre abscissas and weights	1574	1574	! using n gauss legendre abscissas and weights
!	1575	1575	!
! f(in) : the function to integrate	1576	1576	! f(in) : the function to integrate
! a(in) : lower bound	1577	1577	! a(in) : lower bound
! b(in) : higher bound	1578	1578	! b(in) : higher bound
! n(in) : number of gauss legendre pairs	1579	1579	! n(in) : number of gauss legendre pairs
! res(out): the evaluation of the integral	1580	1580	! res(out): the evaluation of the integral
!	1581	1581	!
double precision,external :: f	1582	1582	double precision,external :: f
double precision, intent(in) :: a,b	1583	1583	double precision, intent(in) :: a,b
integer, intent(in):: n	1584	1584	integer, intent(in):: n
double precision, intent(out) :: res	1585	1585	double precision, intent(out) :: res
	1586	1586
double precision, allocatable :: qx(:),qw(:)	1587	1587	double precision, allocatable :: qx(:),qw(:)
double precision :: xm,xr	1588	1588	double precision :: xm,xr
integer :: i	1589	1589	integer :: i
	1590	1590
! First compute n gauss legendre abs and weight	1591	1591	! First compute n gauss legendre abs and weight
allocate(qx(n))	1592	1592	allocate(qx(n))
allocate(qw(n))	1593	1593	allocate(qw(n))
call fvn_d_gauss_legendre(n,qx,qw)	1594	1594	call fvn_d_gauss_legendre(n,qx,qw)
	1595	1595
xm=0.5d0*(b+a)	1596	1596	xm=0.5d0*(b+a)
xr=0.5d0*(b-a)	1597	1597	xr=0.5d0*(b-a)
	1598	1598
res=0.d0	1599	1599	res=0.d0
	1600	1600
do i=1,n	1601	1601	do i=1,n
res=res+qw(i)f(xm+xrqx(i))	1602	1602	res=res+qw(i)f(xm+xrqx(i))
end do	1603	1603	end do
	1604	1604
res=xr*res	1605	1605	res=xr*res
	1606	1606
deallocate(qw)	1607	1607	deallocate(qw)
deallocate(qx)	1608	1608	deallocate(qx)
	1609	1609
end subroutine	1610	1610	end subroutine
	1611	1611
!!!!!!!!!!!!!!!!!!!!!!!!	1612	1612	!!!!!!!!!!!!!!!!!!!!!!!!
!	1613	1613	!
! Simple and double adaptative Gauss Kronrod integration based on	1614	1614	! Simple and double adaptative Gauss Kronrod integration based on
! a modified version of quadpack ( http://www.netlib.org/quadpack	1615	1615	! a modified version of quadpack ( http://www.netlib.org/quadpack
!	1616	1616	!
! Common parameters :	1617	1617	! Common parameters :
!	1618	1618	!
! key (in)	1619	1619	! key (in)
! epsabs	1620	1620	! epsabs
! epsrel	1621	1621	! epsrel
!	1622	1622	!
!	1623	1623	!
!!!!!!!!!!!!!!!!!!!!!!!!	1624	1624	!!!!!!!!!!!!!!!!!!!!!!!!
	1625	1625
subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)	1626	1626	subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
!	1627	1627	!
! Evaluate the integral of function f(x) between a and b	1628	1628	! Evaluate the integral of function f(x) between a and b
!	1629	1629	!
! f(in) : the function	1630	1630	! f(in) : the function
! a(in) : lower bound	1631	1631	! a(in) : lower bound
! b(in) : higher bound	1632	1632	! b(in) : higher bound
! epsabs(in) : desired absolute error	1633	1633	! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error	1634	1634	! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule	1635	1635	! key(in) : gauss kronrod rule
! 1: 7 - 15 points	1636	1636	! 1: 7 - 15 points
! 2: 10 - 21 points	1637	1637	! 2: 10 - 21 points
! 3: 15 - 31 points	1638	1638	! 3: 15 - 31 points
! 4: 20 - 41 points	1639	1639	! 4: 20 - 41 points
! 5: 25 - 51 points	1640	1640	! 5: 25 - 51 points
! 6: 30 - 61 points	1641	1641	! 6: 30 - 61 points
!	1642	1642	!
! limit(in) : maximum number of subintervals in the partition of the	1643	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	1644	1644	! given integration interval (a,b). A value of 500 will give the same
! behaviour as the imsl routine dqdag	1645	1645	! behaviour as the imsl routine dqdag
!	1646	1646	!
! res(out) : estimated integral value	1647	1647	! res(out) : estimated integral value
! abserr(out) : estimated absolute error	1648	1648	! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines	1649	1649	! ier(out) : error flag from quadpack routines
! 0 : no error	1650	1650	! 0 : no error
! 1 : maximum number of subdivisions allowed	1651	1651	! 1 : maximum number of subdivisions allowed
! has been achieved. one can allow more	1652	1652	! has been achieved. one can allow more
! subdivisions by increasing the value of	1653	1653	! subdivisions by increasing the value of
! limit (and taking the according dimension	1654	1654	! limit (and taking the according dimension
! adjustments into account). however, if	1655	1655	! adjustments into account). however, if
! this yield no improvement it is advised	1656	1656	! this yield no improvement it is advised
! to analyze the integrand in order to	1657	1657	! to analyze the integrand in order to
! determine the integration difficulaties.	1658	1658	! determine the integration difficulaties.
! if the position of a local difficulty can	1659	1659	! if the position of a local difficulty can
! be determined (i.e.singularity,	1660	1660	! be determined (i.e.singularity,
! discontinuity within the interval) one	1661	1661	! discontinuity within the interval) one
! will probably gain from splitting up the	1662	1662	! will probably gain from splitting up the
! interval at this point and calling the	1663	1663	! interval at this point and calling the
! integrator on the subranges. if possible,	1664	1664	! integrator on the subranges. if possible,
! an appropriate special-purpose integrator	1665	1665	! an appropriate special-purpose integrator
! should be used which is designed for	1666	1666	! should be used which is designed for
! handling the type of difficulty involved.	1667	1667	! handling the type of difficulty involved.
! 2 : the occurrence of roundoff error is	1668	1668	! 2 : the occurrence of roundoff error is
! detected, which prevents the requested	1669	1669	! detected, which prevents the requested
! tolerance from being achieved.	1670	1670	! tolerance from being achieved.
! 3 : extremely bad integrand behaviour occurs	1671	1671	! 3 : extremely bad integrand behaviour occurs
! at some points of the integration	1672	1672	! at some points of the integration
! interval.	1673	1673	! interval.
! 6 : the input is invalid, because	1674	1674	! 6 : the input is invalid, because
! (epsabs.le.0 and	1675	1675	! (epsabs.le.0 and
! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))	1676	1676	! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
! or limit.lt.1 or lenw.lt.limit*4.	1677	1677	! or limit.lt.1 or lenw.lt.limit*4.
! result, abserr, neval, last are set	1678	1678	! result, abserr, neval, last are set
! to zero.	1679	1679	! to zero.
! except when lenw is invalid, iwork(1),	1680	1680	! except when lenw is invalid, iwork(1),
! work(limit2+1) and work(limit3+1) are	1681	1681	! work(limit2+1) and work(limit3+1) are
! set to zero, work(1) is set to a and	1682	1682	! set to zero, work(1) is set to a and
! work(limit+1) to b.	1683	1683	! work(limit+1) to b.
	1684	1684
implicit none	1685	1685	implicit none
double precision, external :: f	1686	1686	double precision, external :: f
double precision, intent(in) :: a,b,epsabs,epsrel	1687	1687	double precision, intent(in) :: a,b,epsabs,epsrel
integer, intent(in) :: key	1688	1688	integer, intent(in) :: key
integer, intent(in) :: limit	1689	1689	integer, intent(in) :: limit
double precision, intent(out) :: res,abserr	1690	1690	double precision, intent(out) :: res,abserr
integer, intent(out) :: ier	1691	1691	integer, intent(out) :: ier
	1692	1692
double precision, allocatable :: work(:)	1693	1693	double precision, allocatable :: work(:)
integer, allocatable :: iwork(:)	1694	1694	integer, allocatable :: iwork(:)
integer :: lenw,neval,last	1695	1695	integer :: lenw,neval,last
	1696	1696
! imsl value for limit is 500	1697	1697	! imsl value for limit is 500
lenw=limit*4	1698	1698	lenw=limit*4
	1699	1699
allocate(iwork(limit))	1700	1700	allocate(iwork(limit))
allocate(work(lenw))	1701	1701	allocate(work(lenw))
	1702	1702
call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)	1703	1703	call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
	1704	1704
deallocate(work)	1705	1705	deallocate(work)
deallocate(iwork)	1706	1706	deallocate(iwork)
	1707	1707
end subroutine	1708	1708	end subroutine
	1709	1709
	1710	1710
	1711	1711
subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)	1712	1712	subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)
!	1713	1713	!
! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)	1714	1714	! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)
!	1715	1715	!
! f(in) : the function	1716	1716	! f(in) : the function
! a(in) : lower bound	1717	1717	! a(in) : lower bound
! b(in) : higher bound	1718	1718	! b(in) : higher bound
! g(in) : external function describing lower bound for y	1719	1719	! g(in) : external function describing lower bound for y
! h(in) : external function describing higher bound for y	1720	1720	! h(in) : external function describing higher bound for y
! epsabs(in) : desired absolute error	1721	1721	! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error	1722	1722	! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule	1723	1723	! key(in) : gauss kronrod rule
! 1: 7 - 15 points	1724	1724	! 1: 7 - 15 points
! 2: 10 - 21 points	1725	1725	! 2: 10 - 21 points
! 3: 15 - 31 points	1726	1726	! 3: 15 - 31 points
! 4: 20 - 41 points	1727	1727	! 4: 20 - 41 points
! 5: 25 - 51 points	1728	1728	! 5: 25 - 51 points
! 6: 30 - 61 points	1729	1729	! 6: 30 - 61 points
!	1730	1730	!
! limit(in) : maximum number of subintervals in the partition of the	1731	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	1732	1732	! given integration interval (a,b). A value of 500 will give the same
! behaviour as the imsl routine dqdag	1733	1733	! behaviour as the imsl routine dqdag
!	1734	1734	!
! res(out) : estimated integral value	1735	1735	! res(out) : estimated integral value
! abserr(out) : estimated absolute error	1736	1736	! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines	1737	1737	! ier(out) : error flag from quadpack routines
! 0 : no error	1738	1738	! 0 : no error
! 1 : maximum number of subdivisions allowed	1739	1739	! 1 : maximum number of subdivisions allowed
! has been achieved. one can allow more	1740	1740	! has been achieved. one can allow more
! subdivisions by increasing the value of	1741	1741	! subdivisions by increasing the value of
! limit (and taking the according dimension	1742	1742	! limit (and taking the according dimension
! adjustments into account). however, if	1743	1743	! adjustments into account). however, if
! this yield no improvement it is advised	1744	1744	! this yield no improvement it is advised
! to analyze the integrand in order to	1745	1745	! to analyze the integrand in order to
! determine the integration difficulaties.	1746	1746	! determine the integration difficulaties.
! if the position of a local difficulty can	1747	1747	! if the position of a local difficulty can
! be determined (i.e.singularity,	1748	1748	! be determined (i.e.singularity,
! discontinuity within the interval) one	1749	1749	! discontinuity within the interval) one
! will probably gain from splitting up the	1750	1750	! will probably gain from splitting up the
! interval at this point and calling the	1751	1751	! interval at this point and calling the
! integrator on the subranges. if possible,	1752	1752	! integrator on the subranges. if possible,
! an appropriate special-purpose integrator	1753	1753	! an appropriate special-purpose integrator
! should be used which is designed for	1754	1754	! should be used which is designed for
! handling the type of difficulty involved.	1755	1755	! handling the type of difficulty involved.
! 2 : the occurrence of roundoff error is	1756	1756	! 2 : the occurrence of roundoff error is
! detected, which prevents the requested	1757	1757	! detected, which prevents the requested
! tolerance from being achieved.	1758	1758	! tolerance from being achieved.
! 3 : extremely bad integrand behaviour occurs	1759	1759	! 3 : extremely bad integrand behaviour occurs
! at some points of the integration	1760	1760	! at some points of the integration
! interval.	1761	1761	! interval.
! 6 : the input is invalid, because	1762	1762	! 6 : the input is invalid, because
! (epsabs.le.0 and	1763	1763	! (epsabs.le.0 and
! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))	1764	1764	! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
! or limit.lt.1 or lenw.lt.limit*4.	1765	1765	! or limit.lt.1 or lenw.lt.limit*4.
! result, abserr, neval, last are set	1766	1766	! result, abserr, neval, last are set
! to zero.	1767	1767	! to zero.
! except when lenw is invalid, iwork(1),	1768	1768	! except when lenw is invalid, iwork(1),
! work(limit2+1) and work(limit3+1) are	1769	1769	! work(limit2+1) and work(limit3+1) are
! set to zero, work(1) is set to a and	1770	1770	! set to zero, work(1) is set to a and
! work(limit+1) to b.	1771	1771	! work(limit+1) to b.
	1772	1772
implicit none	1773	1773	implicit none
double precision, external:: f,g,h	1774	1774	double precision, external:: f,g,h
double precision, intent(in) :: a,b,epsabs,epsrel	1775	1775	double precision, intent(in) :: a,b,epsabs,epsrel
integer, intent(in) :: key,limit	1776	1776	integer, intent(in) :: key,limit
integer, intent(out) :: ier	1777	1777	integer, intent(out) :: ier
double precision, intent(out) :: res,abserr	1778	1778	double precision, intent(out) :: res,abserr
	1779	1779
	1780	1780
double precision, allocatable :: work(:)	1781	1781	double precision, allocatable :: work(:)
integer, allocatable :: iwork(:)	1782	1782	integer, allocatable :: iwork(:)
integer :: lenw,neval,last	1783	1783	integer :: lenw,neval,last
	1784	1784
! imsl value for limit is 500	1785	1785	! imsl value for limit is 500
lenw=limit*4	1786	1786	lenw=limit*4
allocate(work(lenw))	1787	1787	allocate(work(lenw))
allocate(iwork(limit))	1788	1788	allocate(iwork(limit))
	1789	1789
call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)	1790	1790	call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
	1791	1791
deallocate(iwork)	1792	1792	deallocate(iwork)
deallocate(work)	1793	1793	deallocate(work)
end subroutine	1794	1794	end subroutine
	1795	1795
	1796	1796
	1797	1797
subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)	1798	1798	subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
!	1799	1799	!
! Evaluate the single integral of function f(x,y) for y between a and b with a	1800	1800	! Evaluate the single integral of function f(x,y) for y between a and b with a
! given x value	1801	1801	! given x value
!	1802	1802	!
! This function is used for the evaluation of the double integral fvn_d_integ_2_gk	1803	1803	! This function is used for the evaluation of the double integral fvn_d_integ_2_gk
!	1804	1804	!
! f(in) : the function	1805	1805	! f(in) : the function
! x(in) : x	1806	1806	! x(in) : x
! a(in) : lower bound	1807	1807	! a(in) : lower bound
! b(in) : higher bound	1808	1808	! b(in) : higher bound
! epsabs(in) : desired absolute error	1809	1809	! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error	1810	1810	! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule	1811	1811	! key(in) : gauss kronrod rule
! 1: 7 - 15 points	1812	1812	! 1: 7 - 15 points
! 2: 10 - 21 points	1813	1813	! 2: 10 - 21 points
! 3: 15 - 31 points	1814	1814	! 3: 15 - 31 points
! 4: 20 - 41 points	1815	1815	! 4: 20 - 41 points
! 5: 25 - 51 points	1816	1816	! 5: 25 - 51 points
! 6: 30 - 61 points	1817	1817	! 6: 30 - 61 points
!	1818	1818	!
! limit(in) : maximum number of subintervals in the partition of the	1819	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	1820	1820	! given integration interval (a,b). A value of 500 will give the same
! behaviour as the imsl routine dqdag	1821	1821	! behaviour as the imsl routine dqdag
!	1822	1822	!
! res(out) : estimated integral value	1823	1823	! res(out) : estimated integral value
! abserr(out) : estimated absolute error	1824	1824	! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines	1825	1825	! ier(out) : error flag from quadpack routines
! 0 : no error	1826	1826	! 0 : no error
! 1 : maximum number of subdivisions allowed	1827	1827	! 1 : maximum number of subdivisions allowed
! has been achieved. one can allow more	1828	1828	! has been achieved. one can allow more
! subdivisions by increasing the value of	1829	1829	! subdivisions by increasing the value of
! limit (and taking the according dimension	1830	1830	! limit (and taking the according dimension
! adjustments into account). however, if	1831	1831	! adjustments into account). however, if
! this yield no improvement it is advised	1832	1832	! this yield no improvement it is advised
! to analyze the integrand in order to	1833	1833	! to analyze the integrand in order to
! determine the integration difficulaties.	1834	1834	! determine the integration difficulaties.
! if the position of a local difficulty can	1835	1835	! if the position of a local difficulty can
! be determined (i.e.singularity,	1836	1836	! be determined (i.e.singularity,
! discontinuity within the interval) one	1837	1837	! discontinuity within the interval) one
! will probably gain from splitting up the	1838	1838	! will probably gain from splitting up the
! interval at this point and calling the	1839	1839	! interval at this point and calling the
! integrator on the subranges. if possible,	1840	1840	! integrator on the subranges. if possible,
! an appropriate special-purpose integrator	1841	1841	! an appropriate special-purpose integrator
! should be used which is designed for	1842	1842	! should be used which is designed for
! handling the type of difficulty involved.	1843	1843	! handling the type of difficulty involved.
! 2 : the occurrence of roundoff error is	1844	1844	! 2 : the occurrence of roundoff error is
! detected, which prevents the requested	1845	1845	! detected, which prevents the requested
! tolerance from being achieved.	1846	1846	! tolerance from being achieved.
! 3 : extremely bad integrand behaviour occurs	1847	1847	! 3 : extremely bad integrand behaviour occurs
! at some points of the integration	1848	1848	! at some points of the integration
! interval.	1849	1849	! interval.
! 6 : the input is invalid, because	1850	1850	! 6 : the input is invalid, because
! (epsabs.le.0 and	1851	1851	! (epsabs.le.0 and
! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))	1852	1852	! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
! or limit.lt.1 or lenw.lt.limit*4.	1853	1853	! or limit.lt.1 or lenw.lt.limit*4.
! result, abserr, neval, last are set	1854	1854	! result, abserr, neval, last are set
! to zero.	1855	1855	! to zero.
! except when lenw is invalid, iwork(1),	1856	1856	! except when lenw is invalid, iwork(1),
! work(limit2+1) and work(limit3+1) are	1857	1857	! work(limit2+1) and work(limit3+1) are
! set to zero, work(1) is set to a and	1858	1858	! set to zero, work(1) is set to a and
! work(limit+1) to b.	1859	1859	! work(limit+1) to b.
	1860	1860
implicit none	1861	1861	implicit none
double precision, external:: f	1862	1862	double precision, external:: f
double precision, intent(in) :: x,a,b,epsabs,epsrel	1863	1863	double precision, intent(in) :: x,a,b,epsabs,epsrel
integer, intent(in) :: key,limit	1864	1864	integer, intent(in) :: key,limit
integer, intent(out) :: ier	1865	1865	integer, intent(out) :: ier
double precision, intent(out) :: res,abserr	1866	1866	double precision, intent(out) :: res,abserr
	1867	1867
	1868	1868
double precision, allocatable :: work(:)	1869	1869	double precision, allocatable :: work(:)
integer, allocatable :: iwork(:)	1870	1870	integer, allocatable :: iwork(:)
integer :: lenw,neval,last	1871	1871	integer :: lenw,neval,last
	1872	1872
! imsl value for limit is 500	1873	1873	! imsl value for limit is 500
lenw=limit*4	1874	1874	lenw=limit*4
allocate(work(lenw))	1875	1875	allocate(work(lenw))
allocate(iwork(limit))	1876	1876	allocate(iwork(limit))
	1877	1877
call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)	1878	1878	call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
	1879	1879
deallocate(iwork)	1880	1880	deallocate(iwork)
deallocate(work)	1881	1881	deallocate(work)
end subroutine	1882	1882	end subroutine
	1883	1883
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	1884	1884	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Include the modified quadpack files	1885	1885	! Include the modified quadpack files
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	1886	1886	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
include "fvn_quadpack/dqag_2d_inner.f"	1887	1887	include "fvn_quadpack/dqag_2d_inner.f"
include "fvn_quadpack/dqk15_2d_inner.f"	1888	1888	include "fvn_quadpack/dqk15_2d_inner.f"
include "fvn_quadpack/dqk31_2d_outer.f"	1889	1889	include "fvn_quadpack/dqk31_2d_outer.f"
include "fvn_quadpack/d1mach.f"	1890	1890	include "fvn_quadpack/d1mach.f"
include "fvn_quadpack/dqk31_2d_inner.f"	1891	1891	include "fvn_quadpack/dqk31_2d_inner.f"
include "fvn_quadpack/dqage.f"	1892	1892	include "fvn_quadpack/dqage.f"
include "fvn_quadpack/dqk15.f"	1893	1893	include "fvn_quadpack/dqk15.f"
include "fvn_quadpack/dqk21.f"	1894	1894	include "fvn_quadpack/dqk21.f"
include "fvn_quadpack/dqk31.f"	1895	1895	include "fvn_quadpack/dqk31.f"
include "fvn_quadpack/dqk41.f"	1896	1896	include "fvn_quadpack/dqk41.f"
include "fvn_quadpack/dqk51.f"	1897	1897	include "fvn_quadpack/dqk51.f"
include "fvn_quadpack/dqk61.f"	1898	1898	include "fvn_quadpack/dqk61.f"
include "fvn_quadpack/dqk41_2d_outer.f"	1899	1899	include "fvn_quadpack/dqk41_2d_outer.f"
include "fvn_quadpack/dqk41_2d_inner.f"	1900	1900	include "fvn_quadpack/dqk41_2d_inner.f"
include "fvn_quadpack/dqag_2d_outer.f"	1901	1901	include "fvn_quadpack/dqag_2d_outer.f"
include "fvn_quadpack/dqpsrt.f"	1902	1902	include "fvn_quadpack/dqpsrt.f"
include "fvn_quadpack/dqag.f"	1903	1903	include "fvn_quadpack/dqag.f"
include "fvn_quadpack/dqage_2d_outer.f"	1904	1904	include "fvn_quadpack/dqage_2d_outer.f"
include "fvn_quadpack/dqage_2d_inner.f"	1905	1905	include "fvn_quadpack/dqage_2d_inner.f"
include "fvn_quadpack/dqk51_2d_outer.f"	1906	1906	include "fvn_quadpack/dqk51_2d_outer.f"
include "fvn_quadpack/dqk51_2d_inner.f"	1907	1907	include "fvn_quadpack/dqk51_2d_inner.f"
include "fvn_quadpack/dqk61_2d_outer.f"	1908	1908	include "fvn_quadpack/dqk61_2d_outer.f"
include "fvn_quadpack/dqk21_2d_outer.f"	1909	1909	include "fvn_quadpack/dqk21_2d_outer.f"
include "fvn_quadpack/dqk61_2d_inner.f"	1910	1910	include "fvn_quadpack/dqk61_2d_inner.f"
include "fvn_quadpack/dqk21_2d_inner.f"	1911	1911	include "fvn_quadpack/dqk21_2d_inner.f"
include "fvn_quadpack/dqk15_2d_outer.f"	1912	1912	include "fvn_quadpack/dqk15_2d_outer.f"
	1913	1913
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	1914	1914	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	1915	1915	!
! Trigonometric functions	1916	1916	! Trigonometric functions
!	1917	1917	!
! fvn_z_acos, fvn_z_asin : complex arc cosine and sine	1918	1918	! fvn_z_acos, fvn_z_asin : complex arc cosine and sine
! fvn_d_acosh : arc cosinus hyperbolic	1919	1919	! fvn_d_acosh : arc cosinus hyperbolic
! fvn_d_asinh : arc sinus hyperbolic	1920	1920	! fvn_d_asinh : arc sinus hyperbolic
!	1921	1921	!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	1922	1922	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function fvn_z_acos(z)	1923	1923	function fvn_z_acos(z)
! double complex arccos function adapted from	1924	1924	! double complex arccos function adapted from
! the c gsl library	1925	1925	! the c gsl library
! http://www.gnu.org/software/gsl/	1926	1926	! http://www.gnu.org/software/gsl/
implicit none	1927	1927	implicit none
complex(kind=8) :: fvn_z_acos	1928	1928	complex(kind=8) :: fvn_z_acos
complex(kind=8) :: z	1929	1929	complex(kind=8) :: z
real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1	1930	1930	real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1
real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8	1931	1931	real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8
complex(kind=8),parameter :: i=(0._8,1._8)	1932	1932	complex(kind=8),parameter :: i=(0._8,1._8)
real(kind=8) :: r_res,i_res	1933	1933	real(kind=8) :: r_res,i_res
	1934	1934
rz=dreal(z)	1935	1935	rz=dreal(z)
iz=dimag(z)	1936	1936	iz=dimag(z)
if ( iz == 0._8 ) then	1937	1937	if ( iz == 0._8 ) then
fvn_z_acos=fvn_z_acos_real(rz)	1938	1938	fvn_z_acos=fvn_z_acos_real(rz)
return	1939	1939	return
end if	1940	1940	end if
	1941	1941
x=dabs(rz)	1942	1942	x=dabs(rz)
y=dabs(iz)	1943	1943	y=dabs(iz)
r=fvn_d_hypot(x+1.,y)	1944	1944	r=fvn_d_hypot(x+1.,y)
s=fvn_d_hypot(x-1.,y)	1945	1945	s=fvn_d_hypot(x-1.,y)
a=0.5*(r + s)	1946	1946	a=0.5*(r + s)
b=x/a	1947	1947	b=x/a
y2=y*y	1948	1948	y2=y*y
	1949	1949
if (b <= b_crossover) then	1950	1950	if (b <= b_crossover) then
r_res=dacos(b)	1951	1951	r_res=dacos(b)
else	1952	1952	else
if (x <= 1.) then	1953	1953	if (x <= 1.) then
d=0.5(a+x)(y2/(r+x+1)+(s + (1 - x)))	1954	1954	d=0.5(a+x)(y2/(r+x+1)+(s + (1 - x)))
r_res=datan(dsqrt(d)/x)	1955	1955	r_res=datan(dsqrt(d)/x)
else	1956	1956	else
apx=a+x	1957	1957	apx=a+x
d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))	1958	1958	d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))
r_res=datan((y*dsqrt(d))/x);	1959	1959	r_res=datan((y*dsqrt(d))/x);
end if	1960	1960	end if
end if	1961	1961	end if
	1962	1962
if (a <= a_crossover) then	1963	1963	if (a <= a_crossover) then
if (x < 1.) then	1964	1964	if (x < 1.) then
am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))	1965	1965	am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))
else	1966	1966	else
am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))	1967	1967	am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))
end if	1968	1968	end if
i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));	1969	1969	i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));
else	1970	1970	else
i_res = dlog (a + dsqrt (a*a - 1.));	1971	1971	i_res = dlog (a + dsqrt (a*a - 1.));
end if	1972	1972	end if
if (rz <0.) then	1973	1973	if (rz <0.) then
r_res=fvn_pi-r_res	1974	1974	r_res=fvn_pi-r_res
end if	1975	1975	end if
i_res=-sign(1._8,iz)*i_res	1976	1976	i_res=-sign(1._8,iz)*i_res
fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res)	1977	1977	fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res)
	1978	1978
end function fvn_z_acos	1979	1979	end function fvn_z_acos
	1980	1980
function fvn_z_acos_real(r)	1981	1981	function fvn_z_acos_real(r)
! return the double complex arc cosinus for a	1982	1982	! return the double complex arc cosinus for a
! double precision argument	1983	1983	! double precision argument
implicit none	1984	1984	implicit none
real(kind=8) :: r	1985	1985	real(kind=8) :: r
complex(kind=8) :: fvn_z_acos_real	1986	1986	complex(kind=8) :: fvn_z_acos_real
	1987	1987
if (dabs(r)<=1._8) then	1988	1988	if (dabs(r)<=1._8) then
fvn_z_acos_real=dcmplx(dacos(r))	1989	1989	fvn_z_acos_real=dcmplx(dacos(r))
return	1990	1990	return
end if	1991	1991	end if
if (r < 0._8) then	1992	1992	if (r < 0._8) then
fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r))	1993	1993	fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r))
else	1994	1994	else
fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r))	1995	1995	fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r))
end if	1996	1996	end if
end function	1997	1997	end function
	1998	1998
	1999	1999
function fvn_z_asin(z)	2000	2000	function fvn_z_asin(z)
! double complex arcsin function derived from	2001	2001	! double complex arcsin function derived from
! the c gsl library	2002	2002	! the c gsl library
! http://www.gnu.org/software/gsl/	2003	2003	! http://www.gnu.org/software/gsl/
implicit none	2004	2004	implicit none
complex(kind=8) :: fvn_z_asin	2005	2005	complex(kind=8) :: fvn_z_asin
complex(kind=8) :: z	2006	2006	complex(kind=8) :: z
real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1	2007	2007	real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1
real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8	2008	2008	real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8
real(kind=8) :: r_res,i_res	2009	2009	real(kind=8) :: r_res,i_res
	2010	2010
rz=dreal(z)	2011	2011	rz=dreal(z)
iz=dimag(z)	2012	2012	iz=dimag(z)
if ( iz == 0._8 ) then	2013	2013	if ( iz == 0._8 ) then
! z is real	2014	2014	! z is real
fvn_z_asin=fvn_z_asin_real(rz)	2015	2015	fvn_z_asin=fvn_z_asin_real(rz)
return	2016	2016	return
end if	2017	2017	end if
	2018	2018
x=dabs(rz)	2019	2019	x=dabs(rz)
y=dabs(iz)	2020	2020	y=dabs(iz)
r=fvn_d_hypot(x+1.,y)	2021	2021	r=fvn_d_hypot(x+1.,y)
s=fvn_d_hypot(x-1.,y)	2022	2022	s=fvn_d_hypot(x-1.,y)
a=0.5*(r + s)	2023	2023	a=0.5*(r + s)
b=x/a	2024	2024	b=x/a
y2=y*y	2025	2025	y2=y*y
	2026	2026
if (b <= b_crossover) then	2027	2027	if (b <= b_crossover) then
r_res=dasin(b)	2028	2028	r_res=dasin(b)
else	2029	2029	else
if (x <= 1.) then	2030	2030	if (x <= 1.) then
d=0.5(a+x)(y2/(r+x+1)+(s + (1 - x)))	2031	2031	d=0.5(a+x)(y2/(r+x+1)+(s + (1 - x)))
r_res=datan(x/dsqrt(d))	2032	2032	r_res=datan(x/dsqrt(d))
else	2033	2033	else
apx=a+x	2034	2034	apx=a+x
d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))	2035	2035	d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))
r_res=datan(x/(y*dsqrt(d)));	2036	2036	r_res=datan(x/(y*dsqrt(d)));
end if	2037	2037	end if
end if	2038	2038	end if
	2039	2039
if (a <= a_crossover) then	2040	2040	if (a <= a_crossover) then
if (x < 1.) then	2041	2041	if (x < 1.) then
am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))	2042	2042	am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))
else	2043	2043	else
am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))	2044	2044	am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))
end if	2045	2045	end if
i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));	2046	2046	i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));
else	2047	2047	else
i_res = dlog (a + dsqrt (a*a - 1.));	2048	2048	i_res = dlog (a + dsqrt (a*a - 1.));
end if	2049	2049	end if
r_res=sign(1._8,rz)*r_res	2050	2050	r_res=sign(1._8,rz)*r_res
i_res=sign(1._8,iz)*i_res	2051	2051	i_res=sign(1._8,iz)*i_res
fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res)	2052	2052	fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res)
	2053	2053
end function fvn_z_asin	2054	2054	end function fvn_z_asin
	2055	2055
function fvn_z_asin_real(r)	2056	2056	function fvn_z_asin_real(r)
! return the double complex arc sinus for a	2057	2057	! return the double complex arc sinus for a
! double precision argument	2058	2058	! double precision argument
implicit none	2059	2059	implicit none
real(kind=8) :: r	2060	2060	real(kind=8) :: r
complex(kind=8) :: fvn_z_asin_real	2061	2061	complex(kind=8) :: fvn_z_asin_real
	2062	2062
if (dabs(r)<=1._8) then	2063	2063	if (dabs(r)<=1._8) then
fvn_z_asin_real=dcmplx(dasin(r))	2064	2064	fvn_z_asin_real=dcmplx(dasin(r))
return	2065	2065	return
end if	2066	2066	end if
if (r < 0._8) then	2067	2067	if (r < 0._8) then
fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r))	2068	2068	fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r))
else	2069	2069	else
fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r))	2070	2070	fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r))
end if	2071	2071	end if
end function fvn_z_asin_real	2072	2072	end function fvn_z_asin_real
	2073	2073
function fvn_d_acosh(r)	2074	2074	function fvn_d_acosh(r)
! return the arc hyperbolic cosine	2075	2075	! return the arc hyperbolic cosine
implicit none	2076	2076	implicit none
real(kind=8) :: r	2077	2077	real(kind=8) :: r
real(kind=8) :: fvn_d_acosh	2078	2078	real(kind=8) :: fvn_d_acosh
if (r >=1) then	2079	2079	if (r >=1) then
fvn_d_acosh=dlog(r+dsqrt(r*r-1))	2080	2080	fvn_d_acosh=dlog(r+dsqrt(r*r-1))
else	2081	2081	else
!! TODO : Better error handling!!!!!!	2082	2082	!! TODO : Better error handling!!!!!!
stop "Argument to fvn_d_acosh lesser than 1"	2083	2083	stop "Argument to fvn_d_acosh lesser than 1"
end if	2084	2084	end if
end function fvn_d_acosh	2085	2085	end function fvn_d_acosh
	2086	2086
function fvn_d_asinh(r)	2087	2087	function fvn_d_asinh(r)
! return the arc hyperbolic sine	2088	2088	! return the arc hyperbolic sine
implicit none	2089	2089	implicit none
real(kind=8) :: r	2090	2090	real(kind=8) :: r
real(kind=8) :: fvn_d_asinh	2091	2091	real(kind=8) :: fvn_d_asinh
fvn_d_asinh=dlog(r+dsqrt(r*r+1))	2092	2092	fvn_d_asinh=dlog(r+dsqrt(r*r+1))
end function fvn_d_asinh	2093	2093	end function fvn_d_asinh
	2094	2094
function fvn_d_hypot(a,b)	2095	2095	function fvn_d_hypot(a,b)
implicit none	2096	2096	implicit none
! return the euclidian norm of vector(a,b)	2097	2097	! return the euclidian norm of vector(a,b)
real(kind=8) :: a,b	2098	2098	real(kind=8) :: a,b
real(kind=8) :: fvn_d_hypot	2099	2099	real(kind=8) :: fvn_d_hypot
fvn_d_hypot=dsqrt(aa+bb)	2100	2100	fvn_d_hypot=dsqrt(aa+bb)
end function	2101	2101	end function
	2102	2102
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	2103	2103	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	2104	2104	!
! SPARSE RESOLUTION	2105	2105	! SPARSE RESOLUTION
!	2106	2106	!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!	2107	2107	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	2108	2108	!
! Sparse resolution is done by interfaçing Tim Davi's UMFPACK	2109	2109	! Sparse resolution is done by interfaçing Tim Davi's UMFPACK
! http://www.cise.ufl.edu/research/sparse/SuiteSparse/	2110	2110	! http://www.cise.ufl.edu/research/sparse/SuiteSparse/
! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig	2111	2111	! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig
!	2112	2112	!
! Solve Ax=B using UMFPACK	2113	2113	! Solve Ax=B using UMFPACK
!	2114	2114	!
! Where A is a sparse matrix given in its triplet form	2115	2115	! Where A is a sparse matrix given in its triplet form
! T -> non zero elements	2116	2116	! T -> non zero elements
! Ti,Tj -> row and column index (1-based) of the given elt	2117	2117	! Ti,Tj -> row and column index (1-based) of the given elt
! n : rank of matrix A	2118	2118	! n : rank of matrix A
! nz : number of non zero elts	2119	2119	! nz : number of non zero elts
!	2120	2120	!
! fvn_*_sparse_solve	2121	2121	! fvn_*_sparse_solve
! * = zl : double complex + integer(8)	2122	2122	! * = zl : double complex + integer(8)
! * = zi : double complex + integer(4)	2123	2123	! * = zi : double complex + integer(4)
!	2124	2124	!
subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status)	2125	2125	subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
implicit none	2126	2126	implicit none
integer(8), intent(in) :: n,nz	2127	2127	integer(8), intent(in) :: n,nz
complex(8),dimension(nz),intent(in) :: T	2128	2128	complex(8),dimension(nz),intent(in) :: T
integer(8),dimension(nz),intent(in) :: Ti,Tj	2129	2129	integer(8),dimension(nz),intent(in) :: Ti,Tj
complex(8),dimension(n),intent(in) :: B	2130	2130	complex(8),dimension(n),intent(in) :: B
complex(8),dimension(n),intent(out) :: x	2131	2131	complex(8),dimension(n),intent(out) :: x
integer(8), intent(out) :: status	2132	2132	integer(8), intent(out) :: status
	2133	2133
integer(8),dimension(:),allocatable :: wTi,wTj	2134	2134	integer(8),dimension(:),allocatable :: wTi,wTj
real(8),dimension(:),allocatable :: Tx,Tz	2135	2135	real(8),dimension(:),allocatable :: Tx,Tz
real(8),dimension(:),allocatable :: Ax,Az	2136	2136	real(8),dimension(:),allocatable :: Ax,Az
integer(8),dimension(:),allocatable :: Ap,Ai	2137	2137	integer(8),dimension(:),allocatable :: Ap,Ai
integer(8) :: symbolic,numeric	2138	2138	integer(8) :: symbolic,numeric
real(8),dimension(:),allocatable :: xx,xz,bx,bz	2139	2139	real(8),dimension(:),allocatable :: xx,xz,bx,bz
real(8),dimension(90) :: info	2140	2140	real(8),dimension(90) :: info
real(8),dimension(20) :: control	2141	2141	real(8),dimension(20) :: control
integer(8) :: sys	2142	2142	integer(8) :: sys
	2143	2143
	2144	2144
status=0	2145	2145	status=0
	2146	2146
! we use a working copy of Ti and Tj to perform 1-based to 0-based translation	2147	2147	! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
! Tx and Tz are the real and imaginary parts of T	2148	2148	! Tx and Tz are the real and imaginary parts of T
allocate(wTi(nz),wTj(nz))	2149	2149	allocate(wTi(nz),wTj(nz))
allocate(Tx(nz),Tz(nz))	2150	2150	allocate(Tx(nz),Tz(nz))
Tx=dble(T)	2151	2151	Tx=dble(T)
Tz=aimag(T)	2152	2152	Tz=aimag(T)
wTi=Ti-1	2153	2153	wTi=Ti-1
wTj=Tj-1	2154	2154	wTj=Tj-1
allocate(Ax(nz),Az(nz))	2155	2155	allocate(Ax(nz),Az(nz))
allocate(Ap(n+1),Ai(nz))	2156	2156	allocate(Ap(n+1),Ai(nz))
	2157	2157
! perform the triplet to compressed column form -> Ap,Ai,Ax,Az	2158	2158	! perform the triplet to compressed column form -> Ap,Ai,Ax,Az
call umfpack_zl_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status)	2159	2159	call umfpack_zl_triplet_to_col(n,n,nz,wTi,wTj,Tx,Tz,Ap,Ai,Ax,Az,status)
! if status is not zero a problem has occured	2160	2160	! if status is not zero a problem has occured
if (status /= 0) then	2161	2161	if (status /= 0) then
write(,) "Problem during umfpack_zl_triplet_to_col"	2162	2162	write(,) "Problem during umfpack_zl_triplet_to_col"
endif	2163	2163	endif
	2164	2164
! Define defaults control values	2165	2165	! Define defaults control values
call umfpack_zl_defaults(control)	2166	2166	call umfpack_zl_defaults(control)
	2167	2167
! Symbolic analysis	2168	2168	! Symbolic analysis
call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info)	2169	2169	call umfpack_zl_symbolic(n,n,Ap,Ai,Ax,Az,symbolic, control, info)
! info(1) should be zero	2170	2170	! info(1) should be zero
if (info(1) /= 0) then	2171	2171	if (info(1) /= 0) then
write(,) "Problem during symbolic analysis"	2172	2172	write(,) "Problem during symbolic analysis"
status=info(1)	2173	2173	status=info(1)
endif	2174	2174	endif
	2175	2175
! Numerical factorization	2176	2176	! Numerical factorization
call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info)	2177	2177	call umfpack_zl_numeric (Ap, Ai, Ax, Az, symbolic, numeric, control, info)
! info(1) should be zero	2178	2178	! info(1) should be zero
if (info(1) /= 0) then	2179	2179	if (info(1) /= 0) then
write(,) "Problem during numerical factorization"	2180	2180	write(,) "Problem during numerical factorization"
status=info(1)	2181	2181	status=info(1)
endif	2182	2182	endif
	2183	2183
! free the C symbolic pointer	2184	2184	! free the C symbolic pointer
call umfpack_zl_free_symbolic (symbolic)	2185	2185	call umfpack_zl_free_symbolic (symbolic)
	2186	2186
allocate(bx(n),bz(n),xx(n),xz(n))	2187	2187	allocate(bx(n),bz(n),xx(n),xz(n))
bx=dble(B)	2188	2188	bx=dble(B)
bz=aimag(B)	2189	2189	bz=aimag(B)
sys=0	2190	2190	sys=0