🚀 go-pugleaf

/set/lang/valid7/gcc-intel-S2/g271-0517/bin/gcc -v 
Reading specs from /set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/specs
gcc version 2.7.1

The following program is OK compiled 

/set/lang/valid7/gcc-intel-S2/g271-0517/bin/gcc floquet.c -lm -O0

or with -O2 (O0 and O2 are slightly different), 
but with -O1 the output is significantly different:

diff O1.output O2.output 

< 3.000000 -3580.563852
< 3.100000 -3308.595397
< 3.200000 542.018952
< 3.300000 342.018952
< Numerical Recipes run-time error...
< Singular matrix in routine ludcmp
< ...now exiting to system...

> 3.000000 -2442.598801
> 3.100000 -2454.392663
> 3.200000 -2462.952770
> 3.300000 -2468.893193
> 3.400000 -2472.921185
> 3.500000 -2475.868972
> 3.600000 -2478.731600
> 3.700000 -2482.706057
> 3.800000 -2489.218957
> 3.900000 -2499.926271
> 4.000000 -2516.661841


The program floquet.c, plus the "data" file that needs to reside in the
same directory, follow.

There seems to be no problem at -O1 or -O2 with SPARC.

===== data (16 lines) =====
18
0.75
1.5
4000
4000
-2.0
1.2
12.5663706143592

0
0.05
0.0001
1e-30
1e-4


===== floquet.c =====

# 1 "/set/lang/valid0/src/floquet/floquet.c"
# 1 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/stdio.h" 1 3
 
 

 
 
 

 






#pragma ident	"@(#)stdio.h	1.37	95/08/29 SMI"	

# 1 "/usr/include/sys/feature_tests.h" 1 3
 
 

 
 
 




#pragma ident	"@(#)feature_tests.h	1.7	94/12/06 SMI"





 
















# 17 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/stdio.h" 2 3








typedef unsigned int	size_t;


typedef long	fpos_t;













 


















































































typedef struct	 
{




	int		_cnt;	 
	unsigned char	*_ptr;	 

	unsigned char	*_base;	 
	unsigned char	_flag;	 
	unsigned char	_file;	 
} FILE;


extern FILE		__iob[60 ];



extern FILE		*_lastbuf;
extern unsigned char	*_bufendtab[];

extern unsigned char	 _sibuf[], _sobuf[];




extern int	remove(const char *);
extern int	rename(const char *, const char *);
extern FILE	*tmpfile(void);
extern char	*tmpnam(char *);



extern int	fclose(FILE *);
extern int	fflush(FILE *);
extern FILE	*fopen(const char *, const char *);
extern FILE	*freopen(const char *, const char *, FILE *);
extern void	setbuf(FILE *, char *);
extern int	setvbuf(FILE *, char *, int, size_t);
 
extern int	fprintf(FILE *, const char *, ...);
 
extern int	fscanf(FILE *, const char *, ...);
 
extern int	printf(const char *, ...);
 
extern int	scanf(const char *, ...);
 
extern int	sprintf(char *, const char *, ...);
 
extern int	sscanf(const char *, const char *, ...);





extern int	vfprintf(FILE *, const char *, void *);
extern int	vprintf(const char *, void *);
extern int	vsprintf(char *, const char *, void *);

extern int	fgetc(FILE *);
extern char	*fgets(char *, int, FILE *);
extern int	fputc(int, FILE *);
extern int	fputs(const char *, FILE *);
extern int	getc(FILE *);
extern int	getchar(void);
extern char	*gets(char *);
extern int	putc(int, FILE *);
extern int	putchar(int);
extern int	puts(const char *);
extern int	ungetc(int, FILE *);
extern size_t	fread(void *, size_t, size_t, FILE *);
extern size_t	fwrite(const void *, size_t, size_t, FILE *);
extern int	fgetpos(FILE *, fpos_t *);
extern int	fseek(FILE *, long, int);
extern int	fsetpos(FILE *, const fpos_t *);
extern long	ftell(FILE *);
extern void	rewind(FILE *);
extern void	clearerr(FILE *);
extern int	feof(FILE *);
extern int	ferror(FILE *);
extern void	perror(const char *);

extern int	__filbuf(FILE *);
extern int	__flsbuf(int, FILE *);

 





extern FILE	*fdopen(int, const char *);
extern char	*ctermid(char *);
extern int	fileno(FILE *);



 


# 237 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/stdio.h" 3


 




extern FILE	*popen(const char *, const char *);
extern char	*cuserid(char *);
extern char	*tempnam(const char *, const char *);
extern int	getopt(int, char *const *, const char *);

extern int	getsubopt(char **, char *const *, char **);

extern char	*optarg;
extern int	optind, opterr, optopt;
extern int	getw(FILE *);
extern int	putw(int, FILE *);
extern int	pclose(FILE *);



# 341 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/stdio.h" 3


 
































# 395 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/stdio.h" 3









# 1 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/math.h" 1 3
 










#pragma ident	"@(#)math.h	2.5	95/02/07"















 


typedef union _h_val {
  	unsigned long _i[2];
	double _d;
} _h_val;


extern const _h_val __huge_val;










 
















extern int signgam;




 


enum version {libm_ieee = -1, c_issue_4, ansi_1, strict_ansi};


extern const enum version _lib_version;




struct exception {
	int type;
	char *name;
	double arg1;
	double arg2;
	double retval;
};




























 


extern double acos  (double)  ;
extern double asin  (double)  ;
extern double atan  (double)  ;
extern double atan2  (double, double)  ;
extern double cos  (double)  ;
extern double sin  (double)  ;
extern double tan  (double)  ;

extern double cosh  (double)  ;
extern double sinh  (double)  ;
extern double tanh  (double)  ;

extern double exp  (double)  ;
extern double frexp  (double, int *)  ;
extern double ldexp  (double, int)  ;
extern double log  (double)  ;
extern double log10  (double)  ;
extern double modf  (double, double *)  ;

extern double pow  (double, double)  ;
extern double sqrt  (double)  ;

extern double ceil  (double)  ;
extern double fabs  (double)  ;
extern double floor  (double)  ;
extern double fmod  (double, double)  ;



 


extern double erf  (double)  ;
extern double erfc  (double)  ;
extern double gamma  (double)  ;
extern double hypot  (double, double)  ;
extern int isnan  (double)  ;
extern double j0  (double)  ;
extern double j1  (double)  ;
extern double jn  (int, double)  ;
extern double lgamma  (double)  ;
extern double y0  (double)  ;
extern double y1  (double)  ;
extern double yn  (int, double)  ;



 


extern double acosh  (double)  ;
extern double asinh  (double)  ;
extern double atanh  (double)  ;
extern double cbrt  (double)  ;
extern double logb  (double)  ;
extern double nextafter  (double, double)  ;
extern double remainder  (double, double)  ;
extern double scalb  (double, double)  ;

 


extern double expm1  (double)  ;
extern int ilogb  (double)  ;
extern double log1p  (double)  ;
extern double rint  (double)  ;



 


extern int matherr  (struct exception *)  ;

 


extern double significand  (double)  ;

 


extern double copysign  (double, double)  ;
extern double scalbn  (double, int)  ;

 








 


extern float modff  (float, float *)  ;

# 1 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/floatingpoint.h" 1 3
 
 

 
 
 








#pragma ident	"@(#)floatingpoint.h	2.4 94/06/09"

 




 












# 1 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/sys/ieeefp.h" 1 3
 






#pragma ident	"@(#)ieeefp.h	2.7 94/11/09"





 


# 47 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/sys/ieeefp.h" 3



enum fp_direction_type {	 
	fp_nearest	= 0,
	fp_negative	= 1,
	fp_positive	= 2,
	fp_tozero	= 3
};

enum fp_precision_type {	 
	fp_single	= 0,
	fp_precision_3	= 1,
	fp_double	= 2,
	fp_extended	= 3
};

enum fp_exception_type {	 
	fp_invalid	= 0,
	fp_denormalized	= 1,
	fp_division	= 2,
	fp_overflow	= 3,
	fp_underflow	= 4,
	fp_inexact	= 5
};

enum fp_trap_enable_type {	 
	fp_trap_invalid	= 0,
	fp_trap_denormalized	= 1,
	fp_trap_division	= 2,
	fp_trap_overflow	= 3,
	fp_trap_underflow	= 4,
	fp_trap_inexact	= 5
};


# 122 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/sys/ieeefp.h" 3


enum fp_class_type {		 
	fp_zero		= 0,
	fp_subnormal	= 1,
	fp_normal	= 2,
	fp_infinity   	= 3,
	fp_quiet	= 4,
	fp_signaling	= 5
};






# 35 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/floatingpoint.h" 2 3














typedef int sigfpe_code_type;	 

typedef void (*sigfpe_handler_type)();	 





extern sigfpe_handler_type sigfpe  (sigfpe_code_type, sigfpe_handler_type)  ;

 


typedef float single;			



typedef unsigned extended[3];


typedef long double quadruple;	 

typedef unsigned fp_exception_field_type;
				 



 




typedef char decimal_string[512 ];	
				 

typedef struct {
	enum fp_class_type fpclass;
	int	sign;
	int	exponent;
	decimal_string ds;	 


	int	more;		 


	int 	ndigits;	 


} decimal_record;

enum decimal_form {
	fixed_form,		 


	floating_form		 

};

typedef struct {
	enum fp_direction_type rd;	
				 
	enum decimal_form df;	 

	int ndigits;		 
} decimal_mode;

enum decimal_string_form {	 
	invalid_form,		 
	whitespace_form,	 
	fixed_int_form,		 
	fixed_intdot_form,	 
	fixed_dotfrac_form,	 
	fixed_intdotfrac_form,	 
	floating_int_form,	 	
	floating_intdot_form,	 
	floating_dotfrac_form,	 
	floating_intdotfrac_form, 
	inf_form,		 
	infinity_form,		 
	nan_form,		 
	nanstring_form		 
};

extern void single_to_decimal  (single *, decimal_mode *, decimal_record *,
				fp_exception_field_type *)  ;
extern void double_to_decimal  (double *, decimal_mode *, decimal_record *,
				fp_exception_field_type *)  ;
extern void extended_to_decimal  (extended *, decimal_mode *,
				decimal_record *, fp_exception_field_type *)  ;
extern void quadruple_to_decimal  (quadruple *, decimal_mode *,
				decimal_record *, fp_exception_field_type *)  ;

extern void decimal_to_single  (single *, decimal_mode *, decimal_record *,
				fp_exception_field_type *)  ;
extern void decimal_to_double  (double *, decimal_mode *, decimal_record *,
				fp_exception_field_type *)  ;
extern void decimal_to_extended  (extended *, decimal_mode *,
				decimal_record *, fp_exception_field_type *)  ;
extern void decimal_to_quadruple  (quadruple *, decimal_mode *,
				decimal_record *, fp_exception_field_type *)  ;

extern void string_to_decimal  (char **, int, int, decimal_record *,
				enum decimal_string_form *, char **)  ;
extern void func_to_decimal  (char **, int, int, decimal_record *,
				enum decimal_string_form *, char **,
				int (*)(void), int *, int (*)(int))  ;
extern void file_to_decimal  (char **, int, int, decimal_record *,
				enum decimal_string_form *, char **,
				FILE *, int *)  ;

extern char *seconvert  (single *, int, int *, int *, char *)  ;
extern char *sfconvert  (single *, int, int *, int *, char *)  ;
extern char *sgconvert  (single *, int, int, char *)  ;
extern char *econvert  (double, int, int *, int *, char *)  ;
extern char *fconvert  (double, int, int *, int *, char *)  ;
extern char *gconvert  (double, int, int, char *)  ;
extern char *qeconvert  (quadruple *, int, int *, int *, char *)  ;
extern char *qfconvert  (quadruple *, int, int *, int *, char *)  ;
extern char *qgconvert  (quadruple *, int, int, char *)  ;

extern char *ecvt  (double, int, int *, int *)  ;
extern char *fcvt  (double, int, int *, int *)  ;
extern char *gcvt  (double, int, char *)  ;

 



extern double atof  (const char *)  ;
extern double strtod  (const char *, char **)  ;






# 218 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/math.h" 2 3









# 2 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/stdlib.h" 1 3
 
 

 
 
 




#pragma ident	"@(#)stdlib.h	1.27	95/08/28 SMI"	












typedef	struct {
	int	quot;
	int	rem;
} div_t;

typedef struct {
	long	quot;
	long	rem;
} ldiv_t;


typedef struct {
	long long	quot;
	long long	rem;
} lldiv_t;












typedef long	uid_t;















typedef long int wchar_t;




extern unsigned char	__ctype[];



extern double atof(const char *);
extern int atoi(const char *);
extern long int atol(const char *);
extern double strtod(const char *, char **);
extern long int strtol(const char *, char **, int);
extern unsigned long int strtoul(const char *, char **, int);

extern int rand(void);
extern void srand(unsigned int);





extern void *calloc(size_t, size_t);
extern void free(void *);
extern void *malloc(size_t);
extern void *realloc(void *, size_t);

extern void abort(void);
extern int atexit(void (*)(void));
extern void exit(int);
extern char *getenv(const char *);
extern int system(const char *);

extern void *bsearch(const void *, const void *, size_t, size_t,
	int (*)(const void *, const void *));
extern void qsort(void *, size_t, size_t,
	int (*)(const void *, const void *));

extern int abs(int);
extern div_t div(int, int);
extern long int labs(long);
extern ldiv_t ldiv(long, long);

extern int mbtowc(wchar_t *, const char *, size_t);
extern int mblen(const char *, size_t);
extern int wctomb(char *, wchar_t);

extern size_t mbstowcs(wchar_t *, const char *, size_t);
extern size_t wcstombs(char *, const wchar_t *, size_t);





extern double drand48(void);
extern double erand48(unsigned short *);
extern long jrand48(unsigned short *);
extern void lcong48(unsigned short *);
extern long lrand48(void);
extern long mrand48(void);
extern long nrand48(unsigned short *);
extern unsigned short *seed48(unsigned short *);
extern void srand48(long);
extern int putenv(const char *);
extern void setkey(const char *);






extern void swab(const char *, char *, int);





extern long a64l(const char *);
extern int dup2(int, int);
extern char *ecvt(double, int, int *, int *);
extern char *fcvt(double, int, int *, int *);
extern char *qecvt(long double, int, int *, int *);
extern char *qfcvt(long double, int, int *, int *);
extern char *qgcvt(long double, int, char *);
extern char *getcwd(char *, size_t);
extern char *getlogin(void);
extern int getopt(int, char *const *, const char *);
extern int getsubopt(char **, char *const *, char **);
extern char *optarg;
extern int optind, opterr, optopt;
extern char *getpass(const char *);
extern int getpw(uid_t, char *);
extern char *gcvt(double, int, char *);
extern int isatty(int);
extern char *l64a(long);
extern void *memalign(size_t, size_t);
extern char *mktemp(char *);
extern char *realpath(char *, char *);
extern char *ttyname(int);
extern int ttyslot(void);
extern void *valloc(size_t);
extern char *ptsname(int);
extern int  grantpt(int);
extern int  unlockpt(int);


extern long long atoll(const char *);
extern long long llabs(long long);
extern lldiv_t lldiv(long long, long long);
extern char *lltostr(long long, char *);
extern long long strtoll(const char *, char **, int);
extern unsigned long long strtoull(const char *, char **, int);
extern char *ulltostr(unsigned long long, char *);




# 282 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/stdlib.h" 3







# 3 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid0/src/floquet/nrutil.c" 1



# 1 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/stddef.h" 1 3





 






 

 




 


 





 


# 59 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/stddef.h" 3


 





 


















 





 

 





















typedef int ptrdiff_t;









 




 

 


# 180 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/stddef.h" 3





 




 


# 254 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/stddef.h" 3




 

# 288 "/set/lang/valid7/gcc-intel-S2/g271-0517/lib/gcc-lib/i386-unknown-solaris2.5/2.7.1/include/stddef.h" 3




 









 









# 4 "/set/lang/valid0/src/floquet/nrutil.c" 2





void nrerror(char error_text[])
 
{
	fprintf((&__iob[2]) ,"Numerical Recipes run-time error...\n");
	fprintf((&__iob[2]) ,"%s\n",error_text);
	fprintf((&__iob[2]) ,"...now exiting to system...\n");
	exit(1);
}

float *vector(long nl, long nh)
 
{
	float *v;

	v=(float *)malloc((size_t) ((nh-nl+1+ 1 )*sizeof(float)));
	if (!v) nrerror("allocation failure in vector()");
	return v-nl+ 1 ;
}

int *ivector(long nl, long nh)
 
{
	int *v;

	v=(int *)malloc((size_t) ((nh-nl+1+ 1 )*sizeof(int)));
	if (!v) nrerror("allocation failure in ivector()");
	return v-nl+ 1 ;
}

unsigned char *cvector(long nl, long nh)
 
{
	unsigned char *v;

	v=(unsigned char *)malloc((size_t) ((nh-nl+1+ 1 )*sizeof(unsigned char)));
	if (!v) nrerror("allocation failure in cvector()");
	return v-nl+ 1 ;
}

unsigned long *lvector(long nl, long nh)
 
{
	unsigned long *v;

	v=(unsigned long *)malloc((size_t) ((nh-nl+1+ 1 )*sizeof(long)));
	if (!v) nrerror("allocation failure in lvector()");
	return v-nl+ 1 ;
}

double *dvector(long nl, long nh)
 
{
	double *v;

	v=(double *)malloc((size_t) ((nh-nl+1+ 1 )*sizeof(double)));
	if (!v) nrerror("allocation failure in dvector()");
	return v-nl+ 1 ;
}

float **matrix(long nrl, long nrh, long ncl, long nch)
 
{
	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
	float **m;

	 
	m=(float **) malloc((size_t)((nrow+ 1 )*sizeof(float*)));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m += 1 ;
	m -= nrl;

	 
	m[nrl]=(float *) malloc((size_t)((nrow*ncol+ 1 )*sizeof(float)));
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
	m[nrl] += 1 ;
	m[nrl] -= ncl;

	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

	 
	return m;
}

double **dmatrix(long nrl, long nrh, long ncl, long nch)
 
{
	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
	double **m;

	 
	m=(double **) malloc((size_t)((nrow+ 1 )*sizeof(double*)));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m += 1 ;
	m -= nrl;

	 
	m[nrl]=(double *) malloc((size_t)((nrow*ncol+ 1 )*sizeof(double)));
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
	m[nrl] += 1 ;
	m[nrl] -= ncl;

	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

	 
	return m;
}

int **imatrix(long nrl, long nrh, long ncl, long nch)
 
{
	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
	int **m;

	 
	m=(int **) malloc((size_t)((nrow+ 1 )*sizeof(int*)));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m += 1 ;
	m -= nrl;


	 
	m[nrl]=(int *) malloc((size_t)((nrow*ncol+ 1 )*sizeof(int)));
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
	m[nrl] += 1 ;
	m[nrl] -= ncl;

	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

	 
	return m;
}

float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
	long newrl, long newcl)
 
{
	long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
	float **m;

	 
	m=(float **) malloc((size_t) ((nrow+ 1 )*sizeof(float*)));
	if (!m) nrerror("allocation failure in submatrix()");
	m += 1 ;
	m -= newrl;

	 
	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;

	 
	return m;
}

float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
 



{
	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
	float **m;

	 
	m=(float **) malloc((size_t) ((nrow+ 1 )*sizeof(float*)));
	if (!m) nrerror("allocation failure in convert_matrix()");
	m += 1 ;
	m -= nrl;

	 
	m[nrl]=a-ncl;
	for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
	 
	return m;
}

float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
 
{
	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
	float ***t;

	 
	t=(float ***) malloc((size_t)((nrow+ 1 )*sizeof(float**)));
	if (!t) nrerror("allocation failure 1 in f3tensor()");
	t += 1 ;
	t -= nrl;

	 
	t[nrl]=(float **) malloc((size_t)((nrow*ncol+ 1 )*sizeof(float*)));
	if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
	t[nrl] += 1 ;
	t[nrl] -= ncl;

	 
	t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+ 1 )*sizeof(float)));
	if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
	t[nrl][ncl] += 1 ;
	t[nrl][ncl] -= ndl;

	for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
	for(i=nrl+1;i<=nrh;i++) {
		t[i]=t[i-1]+ncol;
		t[i][ncl]=t[i-1][ncl]+ncol*ndep;
		for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
	}

	 
	return t;
}

void free_vector(float *v, long nl, long nh)
 
{
	free((char* ) (v+nl- 1 ));
}

void free_ivector(int *v, long nl, long nh)
 
{
	free((char* ) (v+nl- 1 ));
}

void free_cvector(unsigned char *v, long nl, long nh)
 
{
	free((char* ) (v+nl- 1 ));
}

void free_lvector(unsigned long *v, long nl, long nh)
 
{
	free((char* ) (v+nl- 1 ));
}

void free_dvector(double *v, long nl, long nh)
 
{
	free((char* ) (v+nl- 1 ));
}

void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
 
{
	free((char* ) (m[nrl]+ncl- 1 ));
	free((char* ) (m+nrl- 1 ));
}

void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
 
{
	free((char* ) (m[nrl]+ncl- 1 ));
	free((char* ) (m+nrl- 1 ));
}

void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
 
{
	free((char* ) (m[nrl]+ncl- 1 ));
	free((char* ) (m+nrl- 1 ));
}

void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
 
{
	free((char* ) (b+nrl- 1 ));
}

void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
 
{
	free((char* ) (b+nrl- 1 ));
}

void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
	long ndl, long ndh)
 
{
	free((char* ) (t[nrl][ncl]+ndl- 1 ));
	free((char* ) (t[nrl]+ncl- 1 ));
	free((char* ) (t+nrl- 1 ));
}

# 614 "/set/lang/valid0/src/floquet/nrutil.c"

# 4 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid0/src/floquet/gaussj.c" 1


# 1 "/set/lang/valid0/src/floquet/nrutil.h" 1



static float sqrarg;


static double dsqrarg;


static double dmaxarg1,dmaxarg2;



static double dminarg1,dminarg2;



static float maxarg1,maxarg2;



static float minarg1,minarg2;



static long lmaxarg1,lmaxarg2;



static long lminarg1,lminarg2;



static int imaxarg1,imaxarg2;



static int iminarg1,iminarg2;







void nrerror(char error_text[]);
float *vector(long nl, long nh);
int *ivector(long nl, long nh);
unsigned char *cvector(long nl, long nh);
unsigned long *lvector(long nl, long nh);
double *dvector(long nl, long nh);
float **matrix(long nrl, long nrh, long ncl, long nch);
double **dmatrix(long nrl, long nrh, long ncl, long nch);
int **imatrix(long nrl, long nrh, long ncl, long nch);
float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
	long newrl, long newcl);
float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);
float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
void free_vector(float *v, long nl, long nh);
void free_ivector(int *v, long nl, long nh);
void free_cvector(unsigned char *v, long nl, long nh);
void free_lvector(unsigned long *v, long nl, long nh);
void free_dvector(double *v, long nl, long nh);
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);
void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);
void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
	long ndl, long ndh);

# 99 "/set/lang/valid0/src/floquet/nrutil.h"



# 3 "/set/lang/valid0/src/floquet/gaussj.c" 2



void gaussj(double **a, int n, double **b, int m)
{
	int *indxc,*indxr,*ipiv;
	int i,icol,irow,j,k,l,ll;
	double big,dum,pivinv,temp;

	indxc=ivector(1,n);
	indxr=ivector(1,n);
	ipiv=ivector(1,n);
	for (j=1;j<=n;j++) ipiv[j]=0;
	for (i=1;i<=n;i++) {
		big=0.0;
		for (j=1;j<=n;j++)
			if (ipiv[j] != 1)
				for (k=1;k<=n;k++) {
					if (ipiv[k] == 0) {
						if (fabs(a[j][k]) >= big) {
							big=fabs(a[j][k]);
							irow=j;
							icol=k;
						}
					} else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1");
				}
		++(ipiv[icol]);
		if (irow != icol) {
			for (l=1;l<=n;l++) {temp=( a[irow][l] );( a[irow][l] )=( a[icol][l] );( a[icol][l] )=temp;} 
			for (l=1;l<=m;l++) {temp=( b[irow][l] );( b[irow][l] )=( b[icol][l] );( b[icol][l] )=temp;} 
		}
		indxr[i]=irow;
		indxc[i]=icol;
		if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix-2");
		pivinv=1.0/a[icol][icol];
		a[icol][icol]=1.0;
		for (l=1;l<=n;l++) a[icol][l] *= pivinv;
		for (l=1;l<=m;l++) b[icol][l] *= pivinv;
		for (ll=1;ll<=n;ll++)
			if (ll != icol) {
				dum=a[ll][icol];
				a[ll][icol]=0.0;
				for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
				for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;
			}
	}
	for (l=n;l>=1;l--) {
		if (indxr[l] != indxc[l])
			for (k=1;k<=n;k++)
				{temp=( a[k][indxr[l]] );( a[k][indxr[l]] )=( a[k][indxc[l]] );( a[k][indxc[l]] )=temp;} ;
	}
	free_ivector(ipiv,1,n);
	free_ivector(indxr,1,n);
	free_ivector(indxc,1,n);
}


 
# 5 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid0/src/floquet/balanc.c" 1



void balanc(double **a, int n)
{
	int last,j,i;
	double s,r,g,f,c,sqrdx;

	sqrdx= 2.0 * 2.0 ;
	last=0;
	while (last == 0) {
		last=1;
		for (i=1;i<=n;i++) {
			r=c=0.0;
			for (j=1;j<=n;j++)
				if (j != i) {
					c += fabs(a[j][i]);
					r += fabs(a[i][j]);
				}
			if (c && r) {
				g=r/ 2.0 ;
				f=1.0;
				s=c+r;
				while (c<g) {
					f *= 2.0 ;
					c *= sqrdx;
				}
				g=r* 2.0 ;
				while (c>g) {
					f /= 2.0 ;
					c /= sqrdx;
				}
				if ((c+r)/f < 0.95*s) {
					last=0;
					g=1.0/f;
					for (j=1;j<=n;j++) a[i][j] *= g;
					for (j=1;j<=n;j++) a[j][i] *= f;
				}
			}
		}
	}
}

 
# 6 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid0/src/floquet/elmhes.c" 1



void elmhes(double **a, int n)
{
	int m,j,i;
	double y,x;

	for (m=2;m<n;m++) {
		x=0.0;
		i=m;
		for (j=m;j<=n;j++) {
			if (fabs(a[j][m-1]) > fabs(x)) {
				x=a[j][m-1];
				i=j;
			}
		}
		if (i != m) {
			for (j=m-1;j<=n;j++) {y=( a[i][j] );( a[i][j] )=( a[m][j] );( a[m][j] )=y;} 
			for (j=1;j<=n;j++) {y=( a[j][i] );( a[j][i] )=( a[j][m] );( a[j][m] )=y;} 
		}
		if (x) {
			for (i=m+1;i<=n;i++) {
				if ((y=a[i][m-1]) != 0.0) {
					y /= x;
					a[i][m-1]=y;
					for (j=m;j<=n;j++)
						a[i][j] -= y*a[m][j];
					for (j=1;j<=n;j++)
						a[j][m] += y*a[j][i];
				}
			}
		}
	}
}

 
# 7 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid0/src/floquet/hqr.c" 1




void hqr(double **a, int n, double wr[], double wi[])
{
	int nn,m,l,k,j,its,i,mmin;
	double z,y,x,w,v,u,t,s,r,q,p,anorm;

	anorm=fabs(a[1][1]);
	for (i=2;i<=n;i++)
		for (j=(i-1);j<=n;j++)
			anorm += fabs(a[i][j]);
	nn=n;
	t=0.0;
	while (nn >= 1) {
		its=0;
		do {
			for (l=nn;l>=2;l--) {
				s=fabs(a[l-1][l-1])+fabs(a[l][l]);
				if (s == 0.0) s=anorm;
				if ((double)(fabs(a[l][l-1]) + s) == s) break;
			}
			x=a[nn][nn];
			if (l == nn) {
				wr[nn]=x+t;
				wi[nn--]=0.0;
			} else {
				y=a[nn-1][nn-1];
				w=a[nn][nn-1]*a[nn-1][nn];
				if (l == (nn-1)) {
					p=0.5*(y-x);
					q=p*p+w;
					z=sqrt(fabs(q));
					x += t;
					if (q >= 0.0) {
						z=p+ (( p ) >= 0.0 ? fabs( z ) : -fabs( z )) ;
						wr[nn-1]=wr[nn]=x+z;
						if (z) wr[nn]=x-w/z;
						wi[nn-1]=wi[nn]=0.0;
					} else {
						wr[nn-1]=wr[nn]=x+p;
						wi[nn-1]= -(wi[nn]=z);
					}
					nn -= 2;
				} else {
					if (its == 30) nrerror("Too many iterations in hqr");
					if (its == 10 || its == 20) {
						t += x;
						for (i=1;i<=nn;i++) a[i][i] -= x;
						s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]);
						y=x=0.75*s;
						w = -0.4375*s*s;
					}
					++its;
					for (m=(nn-2);m>=l;m--) {
						z=a[m][m];
						r=x-z;
						s=y-z;
						p=(r*s-w)/a[m+1][m]+a[m][m+1];
						q=a[m+1][m+1]-z-r-s;
						r=a[m+2][m+1];
						s=fabs(p)+fabs(q)+fabs(r);
						p /= s;
						q /= s;
						r /= s;
						if (m == l) break;
						u=fabs(a[m][m-1])*(fabs(q)+fabs(r));
						v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));
						if ((double)(u+v) == v) break;
					}
					for (i=m+2;i<=nn;i++) {
						a[i][i-2]=0.0;
						if (i != (m+2)) a[i][i-3]=0.0;
					}
					for (k=m;k<=nn-1;k++) {
						if (k != m) {
							p=a[k][k-1];
							q=a[k+1][k-1];
							r=0.0;
							if (k != (nn-1)) r=a[k+2][k-1];
							if ((x=fabs(p)+fabs(q)+fabs(r)) != 0.0) {
								p /= x;
								q /= x;
								r /= x;
							}
						}
						if ((s= (( p ) >= 0.0 ? fabs( sqrt(p*p+q*q+r*r) ) : -fabs( sqrt(p*p+q*q+r*r) )) ) != 0.0) {
							if (k == m) {
								if (l != m)
								a[k][k-1] = -a[k][k-1];
							} else
								a[k][k-1] = -s*x;
							p += s;
							x=p/s;
							y=q/s;
							z=r/s;
							q /= p;
							r /= p;
							for (j=k;j<=nn;j++) {
								p=a[k][j]+q*a[k+1][j];
								if (k != (nn-1)) {
									p += r*a[k+2][j];
									a[k+2][j] -= p*z;
								}
								a[k+1][j] -= p*y;
								a[k][j] -= p*x;
							}
							mmin = nn<k+3 ? nn : k+3;
							for (i=l;i<=mmin;i++) {
								p=x*a[i][k]+y*a[i][k+1];
								if (k != (nn-1)) {
									p += z*a[i][k+2];
									a[i][k+2] -= p*r;
								}
								a[i][k+1] -= p*q;
								a[i][k] -= p;
							}
						}
					}
				}
			}
		} while (l < nn-1);
	}
}

 
# 8 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid0/src/floquet/rk4.c" 1



void rk4(double y[], double dydx[], int n, double x, double h, double yout[],
	void (*derivs)(double, double [], double []))
{
	int i;
	double xh,hh,h6,*dym,*dyt,*yt;

	dym=dvector(1,n);
	dyt=dvector(1,n);
	yt=dvector(1,n);
	hh=h*0.5;
	h6=h/6.0;
	xh=x+hh;
	for (i=1;i<=n;i++) yt[i]=y[i]+hh*dydx[i];
	(*derivs)(xh,yt,dyt);
	for (i=1;i<=n;i++) yt[i]=y[i]+hh*dyt[i];
	(*derivs)(xh,yt,dym);
	for (i=1;i<=n;i++) {
		yt[i]=y[i]+h*dym[i];
		dym[i] += dyt[i];
	}
	(*derivs)(x+h,yt,dyt);
	for (i=1;i<=n;i++)
		yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);
	free_dvector(yt,1,n);
	free_dvector(dyt,1,n);
	free_dvector(dym,1,n);
}

 
# 9 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid0/src/floquet/rkck.c" 1



void rkck(double y[], double dydx[], int n, double x, double h, double yout[],
	double yerr[], void (*derivs)(double, double [], double []))
{
	int i;
	static double a2=0.2,a3=0.3,a4=0.6,a5=1.0,a6=0.875,b21=0.2,
		b31=3.0/40.0,b32=9.0/40.0,b41=0.3,b42 = -0.9,b43=1.2,
		b51 = -11.0/54.0, b52=2.5,b53 = -70.0/27.0,b54=35.0/27.0,
		b61=1631.0/55296.0,b62=175.0/512.0,b63=575.0/13824.0,
		b64=44275.0/110592.0,b65=253.0/4096.0,c1=37.0/378.0,
		c3=250.0/621.0,c4=125.0/594.0,c6=512.0/1771.0,
		dc5 = -277.00/14336.0;
	double dc1=c1-2825.0/27648.0,dc3=c3-18575.0/48384.0,
		dc4=c4-13525.0/55296.0,dc6=c6-0.25;
	double *ak2,*ak3,*ak4,*ak5,*ak6,*ytemp;

	ak2=dvector(1,n);
	ak3=dvector(1,n);
	ak4=dvector(1,n);
	ak5=dvector(1,n);
	ak6=dvector(1,n);
	ytemp=dvector(1,n);
	for (i=1;i<=n;i++)
		ytemp[i]=y[i]+b21*h*dydx[i];
	(*derivs)(x+a2*h,ytemp,ak2);
	for (i=1;i<=n;i++)
		ytemp[i]=y[i]+h*(b31*dydx[i]+b32*ak2[i]);
	(*derivs)(x+a3*h,ytemp,ak3);
	for (i=1;i<=n;i++)
		ytemp[i]=y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
	(*derivs)(x+a4*h,ytemp,ak4);
	for (i=1;i<=n;i++)
		ytemp[i]=y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);
	(*derivs)(x+a5*h,ytemp,ak5);
	for (i=1;i<=n;i++)
		ytemp[i]=y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);
	(*derivs)(x+a6*h,ytemp,ak6);
	for (i=1;i<=n;i++)
		yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
	for (i=1;i<=n;i++)
		yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
	free_dvector(ytemp,1,n);
	free_dvector(ak6,1,n);
	free_dvector(ak5,1,n);
	free_dvector(ak4,1,n);
	free_dvector(ak3,1,n);
	free_dvector(ak2,1,n);
}

 
# 10 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid0/src/floquet/rkqs.c" 1








void rkqs(double y[], double dydx[], int n, double *x, double htry, double eps,
	double yscal[], double *hdid, double *hnext,
	void (*derivs)(double, double [], double []))
{
	void rkck(double y[], double dydx[], int n, double x, double h,
		double yout[], double yerr[], void (*derivs)(double, double [], double []));
	int i;
	double errmax,h,htemp,xnew,*yerr,*ytemp;

	yerr=dvector(1,n);
	ytemp=dvector(1,n);
	h=htry;
	for (;;) {
		rkck(y,dydx,n,*x,h,ytemp,yerr,derivs);
		errmax=0.0;
		for (i=1;i<=n;i++) errmax= (maxarg1=( errmax ),maxarg2=( fabs(yerr[i]/yscal[i]) ),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2)) ;
		errmax /= eps;
		if (errmax > 1.0) {
			htemp= 0.9 *h*pow(errmax,-0.25 );
			h=(h >= 0.0 ? (maxarg1=( htemp ),maxarg2=( 0.1*h ),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))  : (minarg1=( htemp ),minarg2=( 0.1*h ),(minarg1) < (minarg2) ? (minarg1) : (minarg2)) );
			xnew=(*x)+h;
			if (xnew == *x) nrerror("stepsize underflow in rkqs");
			continue;
		} else {
			if (errmax > 1.89e-4 ) *hnext= 0.9 *h*pow(errmax,-0.2 );
			else *hnext=5.0*h;
			*x += (*hdid=h);
			for (i=1;i<=n;i++) y[i]=ytemp[i];
			break;
		}
	}
	free_dvector(ytemp,1,n);
	free_dvector(yerr,1,n);
}





 
# 11 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid0/src/floquet/odeint.c" 1






extern int kmax,kount;
extern double *xp,**yp,dxsav;

void odeint(double ystart[], int nvar, double x1, double x2, double eps, double h1,
	double hmin, int *nok, int *nbad,
	void (*derivs)(double, double [], double []),
	void (*rkqs)(double [], double [], int, double *, double, double, double [],
	double *, double *, void (*)(double, double [], double [])))
{
	int nstp,i;
	double xsav,x,hnext,hdid,h;
	double *yscal,*y,*dydx;

	yscal=dvector(1,nvar);
	y=dvector(1,nvar);
	dydx=dvector(1,nvar);
	x=x1;
	h= (( x2-x1 ) >= 0.0 ? fabs( h1 ) : -fabs( h1 )) ;
	*nok = (*nbad) = kount = 0;
	for (i=1;i<=nvar;i++) y[i]=ystart[i];
	if (kmax > 0) xsav=x-dxsav*2.0;
	for (nstp=1;nstp<= 10000 ;nstp++) {
		(*derivs)(x,y,dydx);
		for (i=1;i<=nvar;i++)
			yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+ 1.0e-30 ;
		if (kmax > 0 && kount < kmax-1 && fabs(x-xsav) > fabs(dxsav)) {
			xp[++kount]=x;
			for (i=1;i<=nvar;i++) yp[i][kount]=y[i];
			xsav=x;
		}
		if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x;
		(*rkqs)(y,dydx,nvar,&x,h,eps,yscal,&hdid,&hnext,derivs);
		if (hdid == h) ++(*nok); else ++(*nbad);
		if ((x-x2)*(x2-x1) >= 0.0) {
			for (i=1;i<=nvar;i++) ystart[i]=y[i];
			if (kmax) {
				xp[++kount]=x;
				for (i=1;i<=nvar;i++) yp[i][kount]=y[i];
			}
			free_dvector(dydx,1,nvar);
			free_dvector(y,1,nvar);
			free_dvector(yscal,1,nvar);
			return;
		}
		if (fabs(hnext) <= hmin) nrerror("Step size too small in odeint");
		h=hnext;
	}
	nrerror("Too many steps in routine odeint");
}



 
# 12 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid0/src/floquet/indexx.c" 1






void indexx(unsigned long n, double arr[], unsigned long indx[])
{
	unsigned long i,indxt,ir=n,itemp,j,k,l=1;
	int jstack=0,*istack;
	double a;

	istack=ivector(1,50 );
	for (j=1;j<=n;j++) indx[j]=j;
	for (;;) {
		if (ir-l < 7 ) {
			for (j=l+1;j<=ir;j++) {
				indxt=indx[j];
				a=arr[indxt];
				for (i=j-1;i>=1;i--) {
					if (arr[indx[i]] <= a) break;
					indx[i+1]=indx[i];
				}
				indx[i+1]=indxt;
			}
			if (jstack == 0) break;
			ir=istack[jstack--];
			l=istack[jstack--];
		} else {
			k=(l+ir) >> 1;
			itemp=( indx[k] );( indx[k] )=( indx[l+1] );( indx[l+1] )=itemp; ;
			if (arr[indx[l+1]] > arr[indx[ir]]) {
				itemp=( indx[l+1] );( indx[l+1] )=( indx[ir] );( indx[ir] )=itemp; 
			}
			if (arr[indx[l]] > arr[indx[ir]]) {
				itemp=( indx[l] );( indx[l] )=( indx[ir] );( indx[ir] )=itemp; 
			}
			if (arr[indx[l+1]] > arr[indx[l]]) {
				itemp=( indx[l+1] );( indx[l+1] )=( indx[l] );( indx[l] )=itemp; 
			}
			i=l+1;
			j=ir;
			indxt=indx[l];
			a=arr[indxt];
			for (;;) {
				do i++; while (arr[indx[i]] < a);
				do j--; while (arr[indx[j]] > a);
				if (j < i) break;
				itemp=( indx[i] );( indx[i] )=( indx[j] );( indx[j] )=itemp; 
			}
			indx[l]=indx[j];
			indx[j]=indxt;
			jstack += 2;
			if (jstack > 50 ) nrerror("NSTACK too small in indexx.");
			if (ir-i+1 >= j-l) {
				istack[jstack]=ir;
				istack[jstack-1]=i;
				ir=j-1;
			} else {
				istack[jstack]=j-1;
				istack[jstack-1]=l;
				l=i;
			}
		}
	}
	free_ivector(istack,1,50 );
}




 
# 13 "/set/lang/valid0/src/floquet/floquet.c" 2

# 1 "/set/lang/valid0/src/floquet/ludcmp.c" 1





void ludcmp(double **a, int n, int *indx, double *d)
{
	int i,imax,j,k;
	double big,dum,sum,temp;
	double *vv;

	vv=dvector(1,n);
	*d=1.0;
	for (i=1;i<=n;i++) {
		big=0.0;
		for (j=1;j<=n;j++)
			if ((temp=fabs(a[i][j])) > big) big=temp;
		if (big == 0.0) nrerror("Singular matrix in routine ludcmp");
		vv[i]=1.0/big;
	}
	for (j=1;j<=n;j++) {
		for (i=1;i<j;i++) {
			sum=a[i][j];
			for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
		}
		big=0.0;
		for (i=j;i<=n;i++) {
			sum=a[i][j];
			for (k=1;k<j;k++)
				sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
			if ( (dum=vv[i]*fabs(sum)) >= big) {
				big=dum;
				imax=i;
			}
		}
		if (j != imax) {
			for (k=1;k<=n;k++) {
				dum=a[imax][k];
				a[imax][k]=a[j][k];
				a[j][k]=dum;
			}
			*d = -(*d);
			vv[imax]=vv[j];
		}
		indx[j]=imax;
		if (a[j][j] == 0.0) a[j][j]= 1.0e-20; ;
		if (j != n) {
			dum=1.0/(a[j][j]);
			for (i=j+1;i<=n;i++) a[i][j] *= dum;
		}
	}
	free_dvector(vv,1,n);
}


 
# 14 "/set/lang/valid0/src/floquet/floquet.c" 2







 

FILE *p_file1,*p_file2,*p_file3;
double eps,ra,rav,gr,grv,k,pr,sc,dz,dz2,dz4,k2,k4,om,period,precision,dtstart,dtmin,accuracy;
double *inverselaplacian[16 +1],*inutile[16 +1],*monodromy[300 +1],wr[300 +1], wi[300 +1];
int kmax,kount,nok,nbad;
double *xp,**yp,dxsav;

int n,N;

 


 

double w(int i, double y[])
{
if (i>1 && i<n-1) return(y[i-1]);
  else if (i==0 || i==n) return(0);
         else if (i==1) return (y[1]/4.0);
                else if (i==n-1) return (y[n-3]/4.0);

}

double th(int i, double y[])
{
if (i>0 && i<n) return(y[i+n-3]);
  else if (i==0 || i==n) return(0);
}

double xi(int i, double y[])
{
if (i>0 && i<n) return(y[i+2*n-4]);
  else if (i==0) return(y[2*n-3]);
         else if (i==n) return(y[3*n-5]);
}

void second_membre(double t, double y[], double dydt[])
{

int i,j;
double du[300 +1];

for (i=2;i<=n-2;++i)
   du[i-1]=(6*w(i,y) - 4*w(i+1,y) - 4*w(i-1,y) + w(i+2,y) + 
       w(i-2,y))/dz4 - (2*k2*(w(i+1,y) + w(i-1,y) - 2*w(i,y)))/ 
     dz2 + k4*w(i,y) - k2*(gr + grv*sin(om*t))* 
     (th(i,y)*(1 + eps) + xi(i,y));

for (i=1;i<=n-3;++i)
   {
   dydt[i]=0;
   for (j=1;j<=n-3;++j)
      dydt[i]+=inverselaplacian[i][j]*du[j];
   }

for (i=1;i<=n-1;++i)
   dydt[i+n-3]=w(i,y) + ((th(i+1,y) + th(i-1,y) - 2*th(i,y))/dz2 - 
       k2*th(i,y))/pr;

for (i=1;i<=n-1;++i)
   dydt[i+2*n-4]=((xi(i+1,y) + xi(i-1,y) - 2*xi(i,y))/dz2 - k2*xi(i,y))/
     sc - (eps*((th(i+1,y) + th(i-1,y) - 2*th(i,y))/dz2 - 
         k2*th(i,y)))/pr;

}

void read_data()
{
p_file1=fopen("data","r");
fscanf(p_file1,"%i %lf %lf %lf %lf %lf %lf %lf\n%i %lf %lf %lf %lf", &n,&pr,&sc,&ra,&rav,&eps,&k,&om,&kmax,&precision,&dtstart,&dtmin,&accuracy);
fclose(p_file1);
}

void write_data(FILE *p_file)
{
fprintf(p_file,"n=%i pr=%lf sc=%lf ra=%lf (gr=%lf) rav=%lf (grv=%lf) eps=%lf om=%lf k=%lf\n", n,pr,sc,ra,gr,rav,grv,eps,om,k);
}

void update_parameters()
{
N=3*n-5;
dz=2.0/((double) n);
dz2=dz*dz;
dz4=dz2*dz2;
k2=k*k;
k4=k2*k2;
gr=ra/(16.0*pr);
grv=sqrt(om*om*rav/(8*pr));
period=2* 3.14159265358979 /om;
}

void initialisation()
{
int i,j;
update_parameters();
if (kmax!=0)
  {
  dxsav=period/kmax;
  xp=malloc((kmax+1)*sizeof(double));
  for (i=0;i<= 300 ;++i)
     yp[i]=malloc((kmax+1)*sizeof(double));
  }
for (i=0;i<= 16 ;++i)
   {
   inverselaplacian[i]=malloc((16 +1)*sizeof(double));
   inutile[i]=malloc(2*sizeof(double));
   }
for (i=0;i<= 300 ;++i)
   monodromy[i]=malloc((300 +1)*sizeof(double));

for (i=1;i<=n-3;++i)
   for (j=1;j<=n-3;++j)
      {
      inutile[j][1]=1;
      inverselaplacian[i][j]=0;
      }
for (i=2;i<=n-4;++i)
   {
   inverselaplacian[i][i]=-(k2+2.0/dz2);
   inverselaplacian[i][i-1]=1.0/dz2;
   inverselaplacian[i][i+1]=1.0/dz2;
   }
inverselaplacian[1][1]=-(k2+7.0/(4.0*dz2));
inverselaplacian[1][2]=1.0/dz2;
inverselaplacian[n-3][n-3]=-(k2+7.0/(4.0*dz2));
inverselaplacian[n-3][n-4]=1.0/dz2;
gaussj(inverselaplacian,n-3,inutile,1);
for (i=0;i<= 16 ;++i)
   free(inutile[i]);
}

void free_memory()
{
int i;
for (i=0;i<= 16 ;++i)
   free(inverselaplacian[i]);
for (i=0;i<= 300 ;++i)
   free(monodromy[i]);
if (kmax!=0)
  {
  free(xp);
  for (i=1;i<=N;++i)
     free(yp[i]);
  }
}

void write_monodromy_matrix(double **monodromy)
{
int i,j;
double y[300 +1];
update_parameters();
for (i=1;i<=N;++i)
   {  
   for (j=1;j<=N;++j)
      y[j]=0;
   y[i]=1;
   odeint(y,N,0,period,precision,dtstart,dtmin,&nok,&nbad,second_membre,rkqs);
   for (j=1;j<=N;++j)
      monodromy[i][j]=y[j];
    
   }
}

void compute_eigenvalues(double **matrix, int n)
{
balanc(matrix,n);
elmhes(matrix,n);
hqr(matrix,n,wr,wi);
}

void write_matrix(FILE *p_file,double **a,int n)
{
int i,j;
for (i=1;i<=n;++i)
   {
   fprintf(p_file,"\n");
   for (j=1;j<=n;++j)
      fprintf(p_file,"%lf ",a[i][j]);
   }
fprintf(p_file,"\n");
}

double deter(double **a, int n)
{
double d;
int j,indx[300 +1];
ludcmp(a,n,indx,&d);
for (j=1;j<=N;j++)
   d*=a[j][j];
return(d);
}

void determinant(FILE *p_file, double lambda, double *var,double var1, double var2, double varstep)
{
for (*var=var1;*var<=var2;*var+=varstep)
  {
  int i;
  write_monodromy_matrix(monodromy);
  for (i=1;i<=N;++i)
     monodromy[i][i]-=lambda;
  fprintf(p_file,"%e %e\n",*var,deter(monodromy,N));
  }
}



double rtsec(double (*func)(double *, double, double), double lambda, double *v, double x1, double x2, double xacc)
{
	void nrerror(char error_text[]);
	int j;
	double fl,f,dx,swap,xl,rts;

	fl=(*func)(v,lambda,x1);
	f=(*func)(v,lambda,x2);
	if (fabs(fl) < fabs(f)) {
		rts=x1;
		xl=x2;
		swap=fl;
		fl=f;
		f=swap;
	} else {
		xl=x1;
		rts=x2;
	}
	for (j=1;j<= 100 ;j++) {
		dx=(xl-rts)*f/(f-fl);
		xl=rts;
		fl=f;
		rts += dx;
		f=(*func)(v,lambda,rts);
		if (fabs(dx) < xacc || f == 0.0) return rts;
	}
	nrerror("Maximum number of iterations exceeded in rtsec");
	return 0.0;
}



double compatibility(double *p_x, double lambda, double x)
{
int i;
double result;
*p_x=x;
write_monodromy_matrix(monodromy);
  for (i=1;i<=N;++i)
     monodromy[i][i]-=lambda;
 
return (deter(monodromy,N));
}

void compute_locus(FILE *p_file, double lambda, double *ivar,double ivar1, double ivar2, double ivarstep, double *dvar,double dvar0,double ddvar)
{
double dvarst=dvar0;
for (*ivar=ivar1;*ivar<=ivar2;*ivar+=ivarstep)
   {
   dvarst=rtsec(compatibility,lambda,dvar,dvarst-ddvar,dvarst+ddvar,accuracy);




   fprintf(p_file,"%lf %lf\n",*ivar,dvarst);
   fprintf((&__iob[2]) ,"%lf %lf\n",*ivar,dvarst);

   }
}


 


main()
{
int i,j;
double wm[300 +1];
unsigned long index[300 +1];
read_data();
initialisation();
write_data((&__iob[2]) );
 





p_file1=fopen("neutral1","w");
p_file2=fopen("neutral2","w");
 
compute_locus(p_file1,1.0,&k,0.9,5.5,0.1,&ra,-5800,200);
 
fclose(p_file1);
fclose(p_file2);

p_file1=fopen("perioddoubling1","w");
p_file2=fopen("perioddoubling2","w");
 


fclose(p_file1);
fclose(p_file2);

 









free_memory();



exit(0);
}