#include "matvecalloc.c"
#include "amino.c"
#include <math.h>

/* definitions */

#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
#define SWAPelmhes(g,h) {y=(g);(g)=(h);(h)=y;}
#define TINY 1.0e-20
#define RADIX 2.0
static double dsqrarg;
#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
static int isqrarg;
#define ISQR(a) ((isqrarg=(a)) == 0.0 ? 0.0 : isqrarg*isqrarg)
static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ? (dmaxarg1):(dmaxarg2))
static double dminarg1,dminarg2;
#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ? (dminarg1):(dminarg2))
static int imaxarg1,imaxarg2;
#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ? (imaxarg1):(imaxarg2))
static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ? (iminarg1):(iminarg2))
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))

#define fabslmy(a) ((a) < 0.0 ? -a : a)
#define SIGNL(a,b) ((b) >= 0.0 ? fabslmy(a) : -fabslmy(a))

/* routines */

double **vec2mat(double *v,long nrl,long nrh);
double  *matcol2vec(long nc,double **m,long nrl, long nrh);
double  *vec_duplicate(double *v,long nrl,long nrh);
int     *ivec_duplicate(int *v,long nrl,long nrh);
double **mat_duplicate(double **m, long nrl, long nrh, long ncl, long nch);
int    **imat_duplicate(int **m, long nrl, long nrh, long ncl, long nch);
double **identity_mat(long n);
int    **identity_imat(long n);
double  *identity_vec(long n);
int     *identity_ivec(long n);
double **vec2diagmat(double *v,long nrl,long nrh);
int    **ivec2idiagmat(int *v,long nrl,long nrh);
double  *vec2double(int *v,long nrl, long nrh);
double **mat2double(int **m,long nrl, long nrh, long ncl, long nch);
int      d_ifsymm(double **m, long nrl, long nrh, long ncl, long nch, double pre);
int      i_ifsymm(int **m, long nrl, long nrh, long ncl, long nch);

double **matT(double **m, long nrl, long nrh, long ncl, long nch);
int    **imatT(int **m, long nrl, long nrh, long ncl, long nch);

double  *dXvec(double a,double *v,long nl,long nh);
int     *iXivec(int a,int *v,long nl,long nh);
double **dXmat(double a,double **m,long,long,long,long);
long double **lddXldmat(long double a,long double **m,long,long,long,long);
int    **iXimat(int a, int **m,long nrl,long nrh,long ncl,long nch);

double  *vecPvec(double *v1,long,long, double *v2,long,long);
int     *ivecPivec(int *v1,long nl1,long nh1,int *v2,long nl2,long nh2);
double **matPmat(double **m1,long,long,long,long,double **m2,long,long,long,long);
int    **imatPimat(int **m1,long,long,long,long,int **m2,long,long,long,long);

double  *vecMvec(double *v1,long,long, double *v2,long,long);
int     *ivecMivec(int *v1,long nl1,long nh1,int *v2,long nl2,long nh2);
double **matMmat(double **m1,long,long,long,long,double **m2,long,long,long,long);
int    **imatMimat(int **m1,long,long,long,long,int **m2,long,long,long,long);

double   vecXvec(double *v1,long,long, double *v2,long,long);
int      ivecXivec(int *v1,long nl1,long nh1,int *v2,long nl2,long nh2);
double **vecTXvecT(double *v1,long nl1,long nh1,double *v2,long nl2,long nh2);
int    **ivecTXivecT(int *v1,long nl1,long nh1,int *v2,long nl2,long nh2);
double   vecXmat_diag(double *v, long nl, long nh, double **m, long nrl,long nrh,long ncl,long nch);
double  *matXvec(double **m,long,long,long,long,double *v,long,long);
int     *imatXivec(int **m,long,long,long,long,int *v,long,long);
double  *imatXvec(int **m,long,long,long,long,double *v,long,long);
double  *vecXmat(double *v,long,long,double **m,long,long,long,long);
int     *ivecXimat(int *v,long,long, int **m,long,long,long,long);
double **matXmat(double **m1,long,long,long,long,double **m2,long,long,long,long);
long double **ldmatXldmat(long double **m1,long,long,long,long,long double **m2,long,long,long,long);
int    **imatXimat(int **m1,long,long,long,long,int **m2,long,long,long,long);
void     UxDxU_1(double **U, double *D, double **U_1, double **m, double  **junk, long n);

double   sqnorm_vec(double *v, long nl,long nh);
double   sqnorm_ivec(int *v, long nl,long nh);
double   sqnorm_mat(double **m, long nrl,long nrh,long ncl,long nch);
double   sqnorm_imat(int **m, long nrl,long nrh,long ncl,long nch);
double	 trace_mat(double **m, long nrl,long nrh,long ncl,long nch);
int	 trace_imat(int **m, long nrl,long nrh,long ncl,long nch);

double **matPOWi(double **m,long nrl,long nrh,long ncl,long nch, int a);
int    **imatPOWi(int **m,long nrl,long nrh,long ncl,long nch, int a);

double **svd_inverse(double **u,double *w,double **vt,long n);
void linsys_gaussjordan(double **a, int n, double **b, int m);
void inv_gaussjordan(double **a, int n);

/*
double **inv_lapack(double **m,long l);
void lse_lapack(double *b,double **A,long m, long n,double *x);
void wls_lapack(double *d,double **A,long m, long n, double *x,double **B, double *y);
void svd_lapack(double **m,long n,double **u,double *w,double **vt);
*/

void     gaussj(double **a, int n, double **b, int m);
void     ludcmp(double **a, int n, int *indx, double *d);
void     lubksb(double **a, int n, int *indx, double b[]);
double  *LUsolve(double **a, long n, double *b);
double **LUinverse(double **a, long n);
double   LUdet(double **a, long n);
void     SVDcmp(double **ai, int m, int n, double **a, double w[], double **v); 
double **SVDinverse(double **a, long m, long n);
double   dpythag(double a, double b);
void     svdcmp(double **a, int m, int n, double w[], double **v);
void     dsvbksb(double **u, double w[], double **v, int m, int n, double b[], double x[]);
double  *SVDsolve(double **a, long m, long n, double *b);
void     choldc(double **a, int n, double p[]);
void     cholsl(double **a, int n, double p[], double b[], double x[]);
double  *CHsolve(double **a, long n, double *b);
void     tred2(double **a, int n, double d[], double e[]);
void     tqli(double d[], double e[], int n, double **z);
void     eigensystem_symm(double **a, long n, double *eval, double **evec, double *helper);
void     balanc(double **a, int n);
void     ldbalanc(long double **a, int n);
void     elmhes(double **a, int n);
void     ldelmhes(long double **a, int n);
void     hqr(double **a, int n, double wr[], double wi[]);
void     ldhqr(long double **a, int n, long double wr[], long double wi[]);
void     eigenvalues(double **a,long n, double *eiR, double *eiI);
void     ldeigenvalues(long double **a,long n, long double *eiR, long double *eiI);
void     mat_prepped(double **a, long n, double ***mp, double *evR, double *evI);
void     ldmat_prepped(long double **a, long n, long double ***mp, long double *evR, long double *evI);

int    **read_aa_imatrix(FILE *mat);
double  *read_aa_vector(FILE *fvec);
int     *read_aa_ivector(FILE *fvec);
void     print_aa_imatrix(int **mat, int n);
void     print_aa_imatrix5(int **mat, int n);
void     fprint_aa_imatrix5(FILE *fo, int **mat, int n);
void     print_aa_imatrix6(int **mat, int n);
void     print_aa_imatrix_dayhoff(int **mat, int n);
void     print_aa_vector(double *vec, int n);
void     print_aa_ivector(int *vec, int n);
void     print_aa_matrix(double **mat, int n);
void     print_aa_ldmatrix(long double **mat, int n);



/* codes */


double **vec2mat(double *v,long nrl,long nrh){
/* vector v[nrl..nrh] is transformed to a matrix m[nrl..nrh][1] */

int i;
double **m;
m=dmatrix(nrl,nrh,1,1);
for(i=nrl;i<=nrh;++i)m[i][1]=v[i];
return m;
}

double  *matcol2vec(long nc, double **m, long nrl, long nrh){
/* transformes the column number nc of the matrix **m from row nrl to row nrh
into a vector: m[nrl.nrh][nc] -> v[nrl.nrh] */

int i;
double *v;
v=dvector(nrl,nrh);
for(i=nrl;i<=nrh;++i)v[i]=m[i][nc];
return v;
}


double  *vec_duplicate(double *v,long nrl,long nrh){
/* vector *v is rewritten to *vnew: v[nrl.nrh] -> vnew[nrl.nrh] */

int i;
double *vnew;
vnew=dvector(nrl,nrh);
for(i=nrl;i<=nrh;++i)vnew[i]=v[i];
return vnew;
}

int  *ivec_duplicate(int *v,long nrl,long nrh){
/* vector *v is rewritten to *vnew: v[nrl.nrh] -> vnew[nrl.nrh] */

int i;
int *vnew;
vnew=ivector(nrl,nrh);
for(i=nrl;i<=nrh;++i)vnew[i]=v[i];
return vnew;
}

double **mat_duplicate(double **m, long nrl, long nrh, long ncl, long nch){
/* matrix **m is rewritten to **mnew: m[nrl.nrh][ncl.nch] -> mnew[nrl.nrh][ncl.nch] */

int i,j;
double **mnew;
mnew=dmatrix(nrl,nrh,ncl,nch);
for(i=nrl;i<=nrh;++i)for(j=ncl;j<=nch;++j)mnew[i][j]=m[i][j];
return mnew;
}

int **imat_duplicate(int **m, long nrl, long nrh, long ncl, long nch){
/* matrix **m is rewritten to **mnew: m[nrl.nrh][ncl.nch] -> mnew[nrl.nrh][ncl.nch] */

int i,j;
int **mnew;
mnew=imatrix(nrl,nrh,ncl,nch);
for(i=nrl;i<=nrh;++i)for(j=ncl;j<=nch;++j)mnew[i][j]=m[i][j];
return mnew;
}

double **identity_mat(long n){
/* allocates square identity matrix of length n+1: m[0.n][0.n] */

int i,j;
double **m;
m=dmatrix(0,n,0,n);
for(i=0;i<=n;++i)for(j=0;j<=n;++j){if(i==j)m[i][j]=1.0;else m[i][j]=0.0;}
return m;
}

int    **identity_imat(long n){
/* allocates square integer identity matrix of length n+1: m[0.n][0.n] */

int i,j;
int **m;
m=imatrix(0,n,0,n);
for(i=0;i<=n;++i)for(j=0;j<=n;++j){if(i==j)m[i][j]=1;else m[i][j]=0;}
return m;
}

double  *identity_vec(long n){
/* allocates identity vector of length n+1: v[0.n] */

int i;
double *v;
v=dvector(0,n);
for(i=0;i<=n;++i)v[i]=1.0;
return v;
}

int     *identity_ivec(long n){
/* allocates identity vector of length n+1: v[0.n] */

int i;
int *v;
v=ivector(0,n);
for(i=0;i<=n;++i)v[i]=1;
return v;
}


double **vec2diagmat(double *v,long nrl,long nrh){
/* writes a vector v as a diagonal of matrix with all other elements equal to 0:
v[nrl.nrh] -> m[nrl.nrh][nrl.nrh]; m[i][i]=v[i], m[i][j]=0 if i!=j */

int i,j;
double **m;
m=dmatrix(nrl,nrh,nrl,nrh);
for(i=nrl;i<=nrh;++i)for(j=nrl;j<=nrh;++j){if(i==j)m[i][j]=v[i];else m[i][j]=0;}
return m;
}


int **ivec2idiagmat(int *v,long nrl,long nrh){
/* writes a vector v as a diagonal of matrix with all other elements equal to 0:
v[nrl.nrh] -> m[nrl.nrh][nrl.nrh]; m[i][i]=v[i], m[i][j]=0 if i!=j */

int i,j;
int **m;
m=imatrix(nrl,nrh,nrl,nrh);
for(i=nrl;i<=nrh;++i)for(j=nrl;j<=nrh;++j){if(i==j)m[i][j]=v[i];else m[i][j]=0;}
return m;
}


double *vec2double(int *v, long nrl, long nrh){
/* converts integer matrix m[nrl.nrh][ncl.nch] to double */

int i;
double *vec;
vec=dvector(nrl,nrh);
for(i=nrl;i<=nrh;++i)vec[i]=1.0*v[i];
return vec;
}

double **mat2double(int **m, long nrl, long nrh, long ncl, long nch){
/* converts integer matrix m[nrl.nrh][ncl.nch] to double */

int i,j;
double **mat;
mat=dmatrix(nrl,nrh,ncl,nch);
for(i=nrl;i<=nrh;++i)for(j=ncl;j<=nch;++j)mat[i][j]=1.0*m[i][j];
return mat;
}

int      d_ifsymm(double **m, long nrl, long nrh, long ncl, long nch, double pre){
/* checks if the matrix m[nrl.nrh][nrl.nrh] is symmetrical with presision pre 
returns 0 if it is and 1 if it is not */

int i,j,flag=0;
for(i=nrl;i<=nrh;++i)for(j=i-nrl+ncl+1;j<=nch;++j)if(fabs(m[i][j]-m[j][i])>=pre)flag=1;
return flag;
}

int      i_ifsymm(int **m, long nrl, long nrh, long ncl, long nch){
/* checks if the matrix m[nrl.nrh][nrl.nrh] is symmetrical
returns 0 if it is and 1 if it is not */

int i,j,flag=0;
for(i=nrl;i<=nrh;++i)for(j=i-nrl+ncl+1;j<=nch;++j)if(m[i][j]!=m[j][i])flag=1;
return flag;
}




double **matT(double **m, long nrl, long nrh, long ncl, long nch){
/* returns the pointer to the transpose matrix mt[ncl.nch][nrl.nrh]*/

int i,j;
double **mt;
mt=dmatrix(ncl,nch,nrl,nrh);
for(i=nrl; i<=nrh; i++)for(j=ncl; j<=nch; j++)mt[j][i]=m[i][j];
return mt;
}

int **imatT(int **m, long nrl, long nrh, long ncl, long nch){
/* returns the pointer to the transpose matrix mt[ncl.nch][nrl.nrh]*/

int i,j;
int **mt;
mt=imatrix(ncl,nch,nrl,nrh);
for(i=nrl; i<=nrh; i++)for(j=ncl; j<=nch; j++)mt[j][i]=m[i][j];
return mt;
}



double  *dXvec(double a,double *v,long nl,long nh){
/* returns a product of double and vector: av[nl.nh]=a.v[nl.nh] */
 
int i;
double *av;
av=dvector(nl,nh);
for(i=nl;i<=nh;i++)av[i]=a*v[i];
return av;
}

int *iXivec(int a, int *v,long nl,long nh){
/* returns a product of integer and integer vector: av[nl.nh]=a.v[nl.nh] */

int i;
int *av;
av=ivector(nl,nh);
for(i=nl;i<=nh;i++)av[i]=a*v[i];
return av;
}

double **dXmat(double a,double **m,long nrl,long nrh,long ncl,long nch){
/* returns a product of double and matrix: mnew[nrl.nrh][ncl.nch]=a.m[nrl.nrh][ncl.nch] */

int i,j;
double **mnew;
mnew=dmatrix(nrl,nrh,ncl,nch);
for(i=nrl;i<=nrh;++i)for(j=ncl;j<=nch;++j)mnew[i][j]=a*m[i][j];
return mnew;
}

long double **lddXldmat(long double a,long double **m,long nrl,long nrh,long ncl,long nch){
/* returns a product of double and matrix: mnew[nrl.nrh][ncl.nch]=a.m[nrl.nrh][ncl.nch] */

int i,j;
long double **mnew;
mnew=ldmatrix(nrl,nrh,ncl,nch);
for(i=nrl;i<=nrh;++i)for(j=ncl;j<=nch;++j)mnew[i][j]=a*m[i][j];
return mnew;
}

int **iXimat(int a, int **m,long nrl,long nrh,long ncl,long nch){
/* returns a product of integer and integer matrix: mnew[nrl.nrh][ncl.nch]=a.m[nrl.nrh][ncl.nch] */

int i,j;
int **mnew;
mnew=imatrix(nrl,nrh,ncl,nch);
for(i=nrl;i<=nrh;++i)for(j=ncl;j<=nch;++j)mnew[i][j]=a*m[i][j];
return mnew;
}



double  *vecPvec(double *v1,long nl1,long nh1, double *v2,long nl2,long nh2){
/* return the vector sum  vPv[nl1.nh1]=v1[nl1.nh1]+v2[nl2.nh2] */

int i;
double *vPv;
if((nh2-nl2)!=(nh1-nl1))nrerror("length of 2 vectors differs in vecPvec");
vPv=dvector(nl1,nh1);
for(i=nl1;i<=nh1;i++)vPv[i]=v1[i]+v2[i-nl1+nl2];
return vPv; 
}

int    *ivecPivec(int *v1,long nl1,long nh1,int *v2,long nl2,long nh2){
/* return the vector sum  vPv[nl1.nh1]=v1[nl1.nh1]+v2[nl2.nh2] */

int i;
int *vPv;
if((nh2-nl2)!=(nh1-nl1))nrerror("length of 2 vectors differs in vecPvec");
vPv=ivector(nl1,nh1);
for(i=nl1;i<=nh1;i++)vPv[i]=v1[i]+v2[i-nl1+nl2];
return vPv;
} 

double **matPmat(double **m1,long nrl1,long nrh1,long ncl1,long nch1,\
double **m2,long nrl2,long nrh2,long ncl2,long nch2){
/*returns the pointer to the matrix mXm[nrl1.nrh1]][ncl1.nch1],
obtained by adding matrix  m1[nrl1.nrh1][ncl1.nch1]
to matrix m2[nrl2.nrh2][ncl2.nch2] */

int i,j;
double **mPm;
if((nch1-ncl1)!=(nch2-ncl2))nrerror("number of columns in the first matrix\
 does not match the number of columns in the second matrix in matPmat");
if((nrh1-nrl1)!=(nrh2-nrl2))nrerror("number of rows in the first matrix\
 does not match the number of rows in the second matrix in matPmat");

mPm=dmatrix(nrl1,nrh1,ncl1,nch1);
for(i=nrl1;i<=nrh1;i++)for(j=ncl1;j<=nch1;j++){
        mPm[i][j]=m1[i][j]+m2[i-nrl1+nrl2][j-ncl1+ncl2];}
return mPm;
}

int **imatPimat(int **m1,long nrl1,long nrh1,long ncl1,long nch1,\
int **m2,long nrl2,long nrh2,long ncl2,long nch2){
/*returns the pointer to the matrix mXm[nrl1.nrh1]][ncl1.nch1], 
obtained by adding matrix  m1[nrl1.nrh1][ncl1.nch1]
to matrix m2[nrl2.nrh2][ncl2.nch2] */

int i,j;
int **mPm;
if((nch1-ncl1)!=(nch2-ncl2))nrerror("number of columns in the first matrix\
 does not match the number of columns in the second matrix in matPmat");
if((nrh1-nrl1)!=(nrh2-nrl2))nrerror("number of rows in the first matrix\
 does not match the number of rows in the second matrix in matPmat");

mPm=imatrix(nrl1,nrh1,ncl1,nch1);
for(i=nrl1;i<=nrh1;i++)for(j=ncl1;j<=nch1;j++){
        mPm[i][j]=m1[i][j]+m2[i-nrl1+nrl2][j-ncl1+ncl2];}
return mPm;
}


double  *vecMvec(double *v1,long nl1,long nh1, double *v2,long nl2,long nh2){
/* return the vector difference vMv[nl1.nh1]=v1[nl1.nh1]-v2[nl2.nh2] */

int i;
double *vMv;
if((nh2-nl2)!=(nh1-nl1))nrerror("length of 2 vectors differs in vecMvec");
vMv=dvector(nl1,nh1);
for(i=nl1;i<=nh1;i++)vMv[i]=v1[i]-v2[i-nl1+nl2];
return vMv; 
}

int    *ivecMivec(int *v1,long nl1,long nh1,int *v2,long nl2,long nh2){
/* return the vector difference  vPv[nl1.nh1]=v1[nl1.nh1]+v2[nl2.nh2] */

int i;
int *vMv;
if((nh2-nl2)!=(nh1-nl1))nrerror("length of 2 vectors differs in vecMvec");
vMv=ivector(nl1,nh1);
for(i=nl1;i<=nh1;i++)vMv[i]=v1[i]-v2[i-nl1+nl2];
return vMv;
} 

double **matMmat(double **m1,long nrl1,long nrh1,long ncl1,long nch1,\
double **m2,long nrl2,long nrh2,long ncl2,long nch2){
/*returns the pointer to the matrix mXm[nrl1.nrh1]][ncl1.nch1],
obtained by subtracting matrix m2[nrl2.nrh2][ncl2.nch2]	 
from the matrix m1[nrl1.nrh1][ncl1.nch1] */

int i,j;
double **mMm;
if((nch1-ncl1)!=(nch2-ncl2))nrerror("number of columns in the first matrix\
 does not match the number of columns in the second matrix in matMmat");
if((nrh1-nrl1)!=(nrh2-nrl2))nrerror("number of rows in the first matrix\
 does not match the number of rows in the second matrix in matMmat");

mMm=dmatrix(nrl1,nrh1,ncl1,nch1);
for(i=nrl1;i<=nrh1;i++)for(j=ncl1;j<=nch1;j++){
        mMm[i][j]=m1[i][j]+m2[i-nrl1+nrl2][j-ncl1+ncl2];}
return mMm;
}

int **imatMimat(int **m1,long nrl1,long nrh1,long ncl1,long nch1,\
int **m2,long nrl2,long nrh2,long ncl2,long nch2){
/*returns the pointer to the matrix mXm[nrl1.nrh1]][ncl1.nch1], 
obtained by subtracting matrix m2[nrl2.nrh2][ncl2.nch2]	 
from the matrix m1[nrl1.nrh1][ncl1.nch1] */

int i,j;
int **mMm;
if((nch1-ncl1)!=(nch2-ncl2))nrerror("number of columns in the first matrix\
 does not match the number of columns in the second matrix in matMmat");
if((nrh1-nrl1)!=(nrh2-nrl2))nrerror("number of rows in the first matrix\
 does not match the number of rows in the second matrix in matMmat");

mMm=imatrix(nrl1,nrh1,ncl1,nch1);
for(i=nrl1;i<=nrh1;i++)for(j=ncl1;j<=nch1;j++){
        mMm[i][j]=m1[i][j]+m2[i-nrl1+nrl2][j-ncl1+ncl2];}
return mMm;
}

double vecXvec(double *v1,long nl1,long nh1,double *v2,long nl2,long nh2){
/* returns the value of a scalar product v1[nl1.nh1].v2[nl2.nh2] */

int i;
double vXv=0;
if((nh2-nl2)!=(nh1-nl1))nrerror("length of 2 vectors differs in vecXvec");
for(i=nl1;i<=nh1;i++)vXv=vXv+v1[i]*v2[i-nl1+nl2];
return vXv;
}

int ivecXivec(int *v1,long nl1,long nh1,int *v2,long nl2,long nh2){
/* returns the value of a scalar product v1[nl1.nh1].v2[nl2.nh2] */

int i;
int vXv=0;
if((nh2-nl2)!=(nh1-nl1))nrerror("length of 2 vectors differs in ivecXivec");
for(i=nl1;i<=nh1;i++)vXv=vXv+v1[i]*v2[i-nl1+nl2];
return vXv;
}

double **vecTXvecT(double *v1,long nl1,long nh1,double *v2,long nl2,long nh2){
/* returns the matrix of a product v1[nl1.nh1].v2[nl2.nh2], 
where v1 is a column vector and v2 is a row vector: vTXvT[nl1.nh1][nl2.nh2] */

int i,j;
double **vTXvT;
vTXvT=dmatrix(nl1,nh1,nl2,nh2);
for(i=nl1;i<=nh1;i++)for(j=nl2;j<=nh2;j++)vTXvT[i][j]=v1[i]*v2[j];
return vTXvT;
}

int **ivecTXivecT(int *v1,long nl1,long nh1,int *v2,long nl2,long nh2){
/* returns the matrix of a  product v1[nl1.nh1].v2[nl2.nh2],
where v1 is a column vector and v2 is a row vector: vTXvT[nl1.nh1][nl2.nh2] */

int i,j;
int **vTXvT;
vTXvT=imatrix(nl1,nh1,nl2,nh2);
for(i=nl1;i<=nh1;i++)for(j=nl2;j<=nh2;j++)vTXvT[i][j]=v1[i]*v2[j];
return vTXvT;
}

double   vecXmat_diag(double *v, long nl, long nh, double **m, long nrl,long nrh,long ncl,long nch){
int i;
double p=0.0;
if((nch-ncl)!=(nrh-nrl))nrerror("non-suare matrix in vecXmat_diag");
if((nh-nl)!=(nch-ncl))nrerror("vector and diagonal lengths are different in vecXmat_diag");
for(i=nl;i<=nh;++i)p=p+v[i]*m[i-nl+nrl][i-nl+ncl];
return p;
}


double *matXvec(double **m,long nrl,long nrh,long ncl,long nch,\
double *v,long nl,long nh){
/*returns the pointer to the vector mXv[nrl.nrh], obtained by multiplying
matrix m[nrl.nrh][ncl.nch] by column vector v[nl.nh] */

int i,j;
double *mXv;
if((nh-nl)!=(nch-ncl))nrerror("vector length does not match the number \
of columns in the matrix in matXvec");
mXv=dvector(nrl,nrh);
for(i=nrl;i<=nrh;i++){
	mXv[i]=0.0;
	for(j=ncl;j<=nch;j++)mXv[i]=mXv[i]+m[i][j]*v[j-ncl+nl];
	}
return mXv;
}

int *imatXivec(int **m,long nrl,long nrh,long ncl,long nch,\
int *v,long nl,long nh){
/*returns the pointer to the vector mXv[nrl.nrh], obtained by multiplying
matrix m[nrl.nrh][ncl.nch] by column vector v[nl.nh] */

int i,j;
int *mXv;
if((nh-nl)!=(nch-ncl))nrerror("vector length does not match the number \
of columns in the matrix in imatXivec");
mXv=ivector(nrl,nrh);
for(i=nrl;i<=nrh;i++){
	mXv[i]=0;
	for(j=ncl;j<=nch;j++)mXv[i]=mXv[i]+m[i][j]*v[j-ncl+nl];
	}
return mXv;
}

double *imatXvec(int **m,long nrl,long nrh,long ncl,long nch,\
double *v,long nl,long nh){
/*returns the pointer to the vector mXv[nrl.nrh], obtained by multiplying
matrix m[nrl.nrh][ncl.nch] by vector v[nl.nh] */

int i,j;
double *mXv;
if((nh-nl)!=(nch-ncl))nrerror("vector length does not match the number \
of columns in the matrix");
mXv=dvector(nrl,nrh);
for(i=nrl;i<=nrh;i++){
	mXv[i]=0.0;
	for(j=ncl;j<=nch;j++)mXv[i]=mXv[i]+1.0*m[i][j]*v[j-ncl+nl];
	}
return mXv;
}



double *vecXmat(double *v,long nl,long nh,double **m,long nrl,\
long nrh,long ncl,long nch){
/*returns the pointer to the vector mXv[nl.nh], obtained by multiplying
row vector v[nl.nh] by matrix m[nrl.nrh][ncl.nch] */

int i,j;
double *mXv;
if((nh-nl)!=(nrh-nrl))nrerror("vector length does not match the number\
 of rows in the matrix");
mXv=dvector(ncl,nch);
for(i=ncl;i<=nch;i++){
	mXv[i]=0.0;
	for(j=nrl;j<=nrh;j++)mXv[i]=mXv[i]+m[j][i]*v[j-nrl+nl];
	}
return mXv;
}

int *ivecXimat(int *v,long nl,long nh,int **m,long nrl,\
long nrh,long ncl,long nch){
/*returns the pointer to the vector mXv[nl.nh], obtained by multiplying
row vector v[nl.nh] by matrix m[nrl.nrh][ncl.nch] */

int i,j;
int *mXv;
if((nh-nl)!=(nrh-nrl))nrerror("vector length does not match the number\
 of rows in the matrix");
mXv=ivector(ncl,nch);
for(i=ncl;i<=nch;i++){
	mXv[i]=0;
	for(j=nrl;j<=nrh;j++)mXv[i]=mXv[i]+m[j][i]*v[j-nrl+nl];
	}
return mXv;
}

double **matXmat(double **m1,long nrl1,long nrh1,long ncl1,long nch1,\
double **m2,long nrl2,long nrh2,long ncl2,long nch2){
/*returns the pointer to the matrix mXm[nrl1.nrh1]][ncl2.nch2], 
obtained by multiplying matrix 	m1[nrl1.nrh1][ncl1.nch1]
by matrix m2[nrl2.nrh2][ncl2.nch2] */

int i,j,k;
double **mXm;

if((nch1-ncl1)!=(nrh2-nrl2))nrerror("number of columns in the first matrix\
 does not match the number of rows in the second matrix");
mXm=dmatrix(nrl1,nrh1,ncl2,nch2);
for(i=nrl1;i<=nrh1;i++)for(j=ncl2;j<=nch2;j++){
	mXm[i][j]=0.0;
	for(k=ncl1;k<=nch1;k++){
	mXm[i][j]=mXm[i][j]+m1[i][k]*m2[k-ncl1+nrl2][j];
		}
	}
return mXm;
}

long double **ldmatXldmat(long double **m1,long nrl1,long nrh1,long ncl1,long nch1,\
long double **m2,long nrl2,long nrh2,long ncl2,long nch2){
/*returns the pointer to the matrix mXm[nrl1.nrh1]][ncl2.nch2], 
obtained by multiplying matrix 	m1[nrl1.nrh1][ncl1.nch1]
by matrix m2[nrl2.nrh2][ncl2.nch2] */

int i,j,k;
long double **mXm;

if((nch1-ncl1)!=(nrh2-nrl2))nrerror("number of columns in the first matrix\
 does not match the number of rows in the second matrix");
mXm=ldmatrix(nrl1,nrh1,ncl2,nch2);
for(i=nrl1;i<=nrh1;i++)for(j=ncl2;j<=nch2;j++){
	mXm[i][j]=0.0;
	for(k=ncl1;k<=nch1;k++){
	mXm[i][j]=mXm[i][j]+m1[i][k]*m2[k-ncl1+nrl2][j];
		}
	}
return mXm;
}


int **imatXimat(int **m1,long nrl1,long nrh1,long ncl1,long nch1,\
int **m2,long nrl2,long nrh2,long ncl2,long nch2){
/*returns the pointer to the matrix mXm[nrl1.nrh1]][ncl2.nch2], 
obtained by multiplying matrix 	m1[nrl1.nrh1][ncl1.nch1]
by matrix m2[nrl2.nrh2][ncl2.nch2] */

int i,j,k;
int **mXm;
if((nch1-ncl1)!=(nrh2-nrl2))nrerror("number of columns in the first matrix\
 does not match the number of rows in the second matrix");
mXm=imatrix(nrl1,nrh1,ncl2,nch2);
for(i=nrl1;i<=nrh1;i++)for(j=ncl2;j<=nch2;j++){
	mXm[i][j]=0;
	for(k=ncl1;k<=nch1;k++){
	mXm[i][j]=mXm[i][j]+m1[i][k]*m2[k-ncl1+nrl2][j];
		}
	}
return mXm;
}



void     UxDxU_1(double **U, double *D, double **U_1, double **m, double **junk,\
 long n){
/* multiplies matrix U[1..n][1..n] by diagonal matrix from vector D[1..n]
and by matrix U_1[1..n][1..n] and puts the result into m[1..n][1..n]
using auxilary matrix junk[1..n][1..n] */
int i,j,k;
for(i=1;i<=n;++i)for(j=1;j<=n;++j)junk[i][j]=D[i]*U_1[i][j];
for(i=1;i<=n;i++)for(j=1;j<=n;j++){m[i][j]=0;
        for(k=1;k<=n;k++){
        m[i][j]=m[i][j]+U[i][k]*junk[k][j];
		}}
}




double   sqnorm_vec(double *v, long nl,long nh){
/* returns the norm (length) of vector v: norm=(Sum_{i=nl}^{nh} v[i]*v[i])^0.5 */

int i;
double norm=0.0;
for(i=nl;i<=nh;++i)norm=norm+v[i]*v[i];
norm=sqrt(norm);
return norm;
}

double   sqnorm_ivec(int *v, long nl,long nh){
/* returns the norm (length) of integer vector v: norm=(Sum_{i=nl}^{nh} v[i]*v[i])^0.5 */

int i;
int norm=0;
double norm1=0;
for(i=nl;i<=nh;++i)norm=norm+v[i]*v[i];
norm1=sqrt(norm);
return norm1;
}
 

double   sqnorm_mat(double **m, long nrl,long nrh,long ncl,long nch){
/* returns the norm of matrix m: 
norm=(Sum_{i=nrl}^{nrh} Sum_{j=ncl}^{nch} m[i][j]*m[i][j])^0.5 */

int i,j;
double norm=0.0;
for(i=nrl;i<=nrh;++i)for(j=ncl;j<=nch;++j)norm=norm+m[i][j]*m[i][j];
norm=sqrt(norm);
return norm; 
}

double   sqnorm_imat(int **m, long nrl,long nrh,long ncl,long nch){
/* returns the norm of integer matrix m: 
norm=(Sum_{i=nrl}^{nrh} Sum_{j=ncl}^{nch} m[i][j]*m[i][j])^0.5 */

int i,j;
int norm=0;
double norm1=0.0;
for(i=nrl;i<=nrh;++i)for(j=ncl;j<=nch;++j)norm=norm+m[i][j]*m[i][j];
norm1=sqrt(norm);
return norm1;
}

double   trace_mat(double **m, long nrl,long nrh,long ncl,long nch){
/* returns the trace of matrix Sum_{i=nrl}^{j=nrh}m[i][i-nrl+ncl] */

int i;
double trace=0;
if((nrh-nrl)!=(nch-ncl))nrerror("Non-diagonal matrix in trace_mat");
for(i=nrl;i<=nrh;++i)trace=trace+m[nrl][i-nrl+ncl];
return trace;
}

int   trace_imat(int **m, long nrl,long nrh,long ncl,long nch){
/* returns the trace of integer matrix Sum_{i=nrl}^{j=nrh}m[i][i-nrl+ncl] */

int i;
int trace=0;
if((nrh-nrl)!=(nch-ncl))nrerror("Non-diagonal matrix in trace_imat");
for(i=nrl;i<=nrh;++i)trace=trace+m[nrl][i-nrl+ncl];
return trace;  
}

double **matPOWi(double **m,long nrl,long nrh,long ncl,long nch, int a){
/* returns the pointer to the matrix mPOWi[nrl.nrh][ncl.nch]=
m[nrl.nrh][ncl.nch]^a for positive integer a */

int i,j;
double **mPOWi;
if(a<=0)nrerror("non-positive power in matPOWi");
if((nch-ncl)!=(nrh-nrl))nrerror("non-square matrix in matPOWi");
mPOWi=dmatrix(nrl,nrh,ncl,nch);
for(i=nrl;i<=nrh;++i)for(j=ncl;j<=nch;++j){if((i-nrl+ncl)==j)mPOWi[i][j]=1.0;
else mPOWi[i][j]=0.0;}
for(i=1;i<=a;++i){
mPOWi=matXmat(mPOWi,nrl,nrh,ncl,nch,m,nrl,nrh,ncl,nch);
		 }
return mPOWi;
}

int    **imatPOWi(int **m,long nrl,long nrh,long ncl,long nch, int a){
/* returns the pointer to the matrix mPOWi[nrl.nrh][ncl.nch]=
m[nrl.nrh][ncl.nch]^a for positive integer a */  

int i,j;
int **mPOWi;
if(a<=0)nrerror("non-positive power in matPOWi");
if((nch-ncl)!=(nrh-nrl))nrerror("non-square matrix in matPOWi");
mPOWi=imatrix(nrl,nrh,ncl,nch);
for(i=nrl;i<=nrh;++i)for(j=ncl;j<=nch;++j){if((i-nrl+ncl)==j)mPOWi[i][j]=1;
else mPOWi[i][j]=0;}
for(i=1;i<=a;++i){
mPOWi=imatXimat(mPOWi,nrl,nrh,ncl,nch,m,nrl,nrh,ncl,nch);
                 }
return mPOWi;
} 


double **svd_inverse(double **u,double *w,double **vt,long l){

/* inversion of decomposed square matrix **u *w *vt */

int i,j;
double ww,w1;
double **m, **ut, **v, **c, **temp, **cc;
m=dmatrix(1,l,1,l);
ut=dmatrix(1,l,1,l);
v=dmatrix(1,l,1,l);
c=dmatrix(1,l,1,l);
cc=dmatrix(1,l,1,l);
temp=dmatrix(1,l,1,l);

for(i=1;i<=l;++i){
	if(fabs(w[i])<0.0000000000001){
		printf("\n one of the eigenvalues is %25.20f\n",w[i]);
		nrerror("Possibly singular matrix");
		}
	c[i][i]=1/w[i];
	cc[i][i]=w[i];
	}
ut=matT(u,1,l,1,l);
v=matT(vt,1,l,1,l);
m=matXmat(c,1,l,1,l,ut,1,l,1,l);
m=matXmat(v,1,l,1,l,m,1,l,1,l);

temp=matXmat(cc,1,l,1,l,vt,1,l,1,l);
temp=matXmat(u,1,l,1,l,temp,1,l,1,l);


temp=matXmat(m,1,l,1,l,temp,1,l,1,l);

ww=0.0;
for(i=1;i<=4;++i)for(j=1;j<=4;++j){
	if(i==j)w1=temp[i][j]-1.0;
	else w1=temp[i][j];
	ww=ww+w1*w1;
	}
if(ww<=0)ww=0;
else ww=sqrt(ww);
if(ww>0.0000000001)nrerror("the precision of inversion is unsatisfactory");

return m;
}


void linsys_gaussjordan(double **a, int n, double **b, int m){

/* inverse matrix will be in **a solutions will be in **b */

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++)SWAP(a[irow][l],a[icol][l])
for(l=1;l<=m;l++)SWAP(b[irow][l],b[icol][l])
}
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++)
			SWAP(a[k][indxr[l]],a[k][indxc[l]]);
}
free_ivector(ipiv,1,n);
free_ivector(indxr,1,n);
free_ivector(indxc,1,n);
}



void inv_gaussjordan(double **a, int n){
/* inverse matrix will be in **a */

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++)SWAP(a[irow][l],a[icol][l])
}
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(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=n;l>=1;l--){
	if(indxr[l] != indxc[l])
		for(k=1;k<=n;k++)
			SWAP(a[k][indxr[l]],a[k][indxc[l]]);
}
free_ivector(ipiv,1,n);
free_ivector(indxr,1,n);
free_ivector(indxc,1,n);
}



/* LAPACK-using routines */

/*
double **inv_lapack(double **m,long l){
int i;
double **m_1,**u,**vt,*w;
u=dmatrix(1,l,1,l);
vt=dmatrix(1,l,1,l);
m_1=dmatrix(1,l,1,l);
w=dvector(1,l);
svd_lapack(m,l,u,w,vt);
m_1=svd_inverse(u,w,vt,l);
return m_1;
}

void lse_lapack(double *b,double **A,long m, long n,double *x){
int i,j;
double *a,*x1,*work,*s,rcond;
int nrhs,lda,ldb,rank,lwork,info;
lwork=10*m;
nrhs=1;
lda=m;
ldb=m;
rcond=0.00000000001;
s=dvector(0,n-1);
a=dvector(0,m*n-1);
x1=dvector(0,m-1);
work=dvector(0,lwork-1);
a=mat2forvect(A,1,m,1,n);
for(i=1;i<=m;++i)x1[i-1]=b[i];
dgelss_(&m,&n,&nrhs,a,&lda,x1,&ldb,s,&rcond,&rank,work,&lwork,&info);
if(info!=0)nrerror("error occured in DGELSS in LAPACK somewhere");
for(i=1;i<=m;++i)x[i]=x1[i-1];
}


void wls_lapack(double *d,double **A,long m, long n, double *x,double \
**B,double *y){ int i;
double *d1,*A1,*x1,*y1,*B1,*work;
int lda,ldb,lwork,info;
lwork=20*m;
lda=m;
ldb=m;
x1=dvector(0,n-1);
y1=dvector(0,m-1);
A1=dvector(0,m*n-1);
d1=dvector(0,m-1);
B1=dvector(0,m*m-1);
work=dvector(0,lwork-1);

A1=mat2forvect(A,1,m,1,n);

B1=mat2forvect(B,1,m,1,m);
for(i=1;i<=m;++i)d1[i-1]=d[i]; 

dggglm_(&m,&n,&m,A1,&lda,B1,&ldb,d1,x1,y1,work,&lwork,&info);

if(info!=0)nrerror("error occured in DGGGLM in LAPACK somewhere");
for(i=1;i<=n;++i)x[i]=x1[i-1];
for(i=1;i<=m;++i)y[i]=y1[i-1]; 
}


void svd_lapack(double **m,long l,double **u,double *w,double **vt){

*//* svd decomposition of **m[1..l][1..l] into **u,*w,and *vt *//*

int i,j,ll;
double ww,w1;
char jobu='A';
char jobvt='A';
int mm,n,lda,ldu,ldvt,info,lwork;
double *d,*u1,*vt1,*work,*c,**temp;

mm=l;n=l;lda=l;ldu=l;ldvt=l;lwork=100*l;ll=l*l;

d=dvector(1,ll);
u1=dvector(1,ll);
vt1=dvector(1,ll);
work=dvector(1,lwork);
c=dvector(1,l);
temp=dmatrix(1,l,1,l);

d=mat2forvect(m,1,l,1,l);
dgesvd_(&jobu,&jobvt,&mm,&n,d,&lda,c,u1,&ldu,vt1,&ldvt,work,&lwork,&info);
if(info!=0)nrerror("error occured in LAPACK somewhere");

for(i=1;i<=l;++i)w[i]=c[i-1];
forvect2mat(u1,1,l,1,l,u);
forvect2mat(vt1,1,l,1,l,vt);


temp=foreigen2diagmat(c,1,l);


temp=matXmat(temp,1,l,1,l,vt,1,l,1,l);
temp=matXmat(u,1,l,1,l,temp,1,l,1,l);

ww=0;
for(i=1;i<=4;++i)for(j=1;j<=4;++j){
	w1=(temp[i][j]-m[i][j]);
	ww=ww+w1*w1;
	}
if(ww<=0)ww==0;
else ww=sqrt(ww);
if(ww>0.00000001)nrerror("the precision of SVD decomposition is unsatisfactory");

}

*/

/* end of LAPACK-using routines */





void gaussj(double **a, int n, double **b, int m)

/* 

A * x=B_k, k=1,..,m 
Linear system solved by Gauss-Jordan elimination;
a[1..n][1..n] input matrix A
b[1..n][1..m] input of m right-hand side vectors
a is replaced by its inverse and 
b is replaced by corresponding set of solutions

*/

{
	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++) SWAP(a[irow][l],a[icol][l])
			for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
		}
		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++)
				SWAP(a[k][indxr[l]],a[k][indxc[l]]);
	}
	free_ivector(ipiv,1,n);
	free_ivector(indxr,1,n);
	free_ivector(indxc,1,n);
}





void ludcmp(double **a, int n, int *indx, double *d)

/*
internal routine for LU decomposition used somewhere else
*/

{
	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]=TINY;nrerror("POSSIBLY SINGULAR MATRIX in ludcmp: FATAL");exit(0);}
		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);
}




void lubksb(double **a, int n, int *indx, double b[])

/*
internal routine for LU decomposition used somewhere else
*/

{
	int i,ii=0,ip,j;
	double sum;

	for (i=1;i<=n;i++) {
		ip=indx[i];
		sum=b[ip];
		b[ip]=b[i];
		if (ii)
			for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
		else if (sum) ii=i;
		b[i]=sum;
	}
	for (i=n;i>=1;i--) {
		sum=b[i];
		for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
		b[i]=sum/a[i][i];
	}
}


double *LUsolve(double **a, long n, double *b){

/*
Solves system a*x=b, a[1..n][1..n], b[1..n]
returns x, x[1..n]
does not destroy originals
Good way to solve non-singular systems
Uses LU decomposition
*/

int *indx;
double **a1,*x,d=0.0;

indx=ivector(1,n);
a1=mat_duplicate(a,1,n,1,n);
x=vec_duplicate(b,1,n);

ludcmp(a1,n,indx,&d);
lubksb(a1,n,indx,x);

return x;
}


double **LUinverse(double **a, long n){
/*
inverses a[1..n][1..n] using LU decomposition
does not destroy originals
*/

double **y,d=0,*col,**a1;
int i,j,*indx;

indx=ivector(1,n);
col=dvector(1,n);
y=dmatrix(1,n,1,n);
a1=mat_duplicate(a,1,n,1,n);
ludcmp(a1,n,indx,&d);
for(j=1;j<=n;j++){
	for(i=1;i<=n;i++)col[i]=0.0;
	col[j]=1.0;
	lubksb(a1,n,indx,col);
	for(i=1;i<=n;i++)y[i][j]=col[i];
		 }
return y;
}


double LUdet(double **a, long n){

/*
calculates the determinant of a[1..n][1..n] using LU decomposition
does not destroy original matrix
*/

double **a1,d=0.0;
int j, *indx;

indx=ivector(1,n);
a1=mat_duplicate(a,1,n,1,n);
ludcmp(a1,n,indx,&d);
for(j=1;j<=n;j++)d *= a1[j][j];
return d;
}


void SVDcmp(double **ai, int m, int n, double **a, double w[], double **v)

/*
ai = a * diag (w) * vT
Given matrix A as ai[1..m][1..n] computes SVD as A=U*W*VT
matrix U is in a[1..m][1..n]
singular values are in vector w[1..n]
matrix V (not its transpose) is in v[1..n][1..n]
*/

{
	int flag,i,its,j,jj,k,l,nm;
	double anorm,c,f,g,h,s,scale,x,y,z,*rv1;
	rv1=dvector(1,n);
	for(i=1;i<=m;++i)for(j=1;j<=n;++j)a[i][j]=ai[i][j];
	g=scale=anorm=0.0;
	for (i=1;i<=n;i++) {
		l=i+1;
		rv1[i]=scale*g;
		g=s=scale=0.0;
		if (i <= m) {
			for (k=i;k<=m;k++) scale += fabs(a[k][i]);
			if (scale) {
				for (k=i;k<=m;k++) {
					a[k][i] /= scale;
					s += a[k][i]*a[k][i];
				}
				f=a[i][i];
				g = -SIGN(sqrt(s),f);
				h=f*g-s;
				a[i][i]=f-g;
				for (j=l;j<=n;j++) {
					for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
					f=s/h;
					for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
				}
				for (k=i;k<=m;k++) a[k][i] *= scale;
			}
		}
		w[i]=scale *g;
		g=s=scale=0.0;
		if (i <= m && i != n) {
			for (k=l;k<=n;k++) scale += fabs(a[i][k]);
			if (scale) {
				for (k=l;k<=n;k++) {
					a[i][k] /= scale;
					s += a[i][k]*a[i][k];
				}
				f=a[i][l];
				g = -SIGN(sqrt(s),f);
				h=f*g-s;
				a[i][l]=f-g;
				for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
				for (j=l;j<=m;j++) {
					for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
					for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
				}
				for (k=l;k<=n;k++) a[i][k] *= scale;
			}
		}
		anorm=DMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
	}
	for (i=n;i>=1;i--) {
		if (i < n) {
			if (g) {
				for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g;
				for (j=l;j<=n;j++) {
					for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
					for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
				}
			}
			for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
		}
		v[i][i]=1.0;
		g=rv1[i];
		l=i;
	}
	for (i=IMIN(m,n);i>=1;i--) {
		l=i+1;
		g=w[i];
		for (j=l;j<=n;j++) a[i][j]=0.0;
		if (g) {
			g=1.0/g;
			for (j=l;j<=n;j++) {
				for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
				f=(s/a[i][i])*g;
				for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
			}
			for (j=i;j<=m;j++) a[j][i] *= g;
		} else for (j=i;j<=m;j++) a[j][i]=0.0;
		++a[i][i];
	}
	for (k=n;k>=1;k--) {
		for (its=1;its<=30;its++) {
			flag=1;
			for (l=k;l>=1;l--) {
				nm=l-1;
				if ((double)(fabs(rv1[l])+anorm) == anorm) {
					flag=0;
					break;
				}
				if ((double)(fabs(w[nm])+anorm) == anorm) break;
			}
			if (flag) {
				c=0.0;
				s=1.0;
				for (i=l;i<=k;i++) {
					f=s*rv1[i];
					rv1[i]=c*rv1[i];
					if ((double)(fabs(f)+anorm) == anorm) break;
					g=w[i];
					h=dpythag(f,g);
					w[i]=h;
					h=1.0/h;
					c=g*h;
					s = -f*h;
					for (j=1;j<=m;j++) {
						y=a[j][nm];
						z=a[j][i];
						a[j][nm]=y*c+z*s;
						a[j][i]=z*c-y*s;
					}
				}
			}
			z=w[k];
			if (l == k) {
				if (z < 0.0) {
					w[k] = -z;
					for (j=1;j<=n;j++) v[j][k] = -v[j][k];
				}
				break;
			}
			if (its == 30) nrerror("no convergence in 30 dsvdcmp iterations");
			x=w[l];
			nm=k-1;
			y=w[nm];
			g=rv1[nm];
			h=rv1[k];
			f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
			g=dpythag(f,1.0);
			f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
			c=s=1.0;
			for (j=l;j<=nm;j++) {
				i=j+1;
				g=rv1[i];
				y=w[i];
				h=s*g;
				g=c*g;
				z=dpythag(f,h);
				rv1[j]=z;
				c=f/z;
				s=h/z;
				f=x*c+g*s;
				g = g*c-x*s;
				h=y*s;
				y *= c;
				for (jj=1;jj<=n;jj++) {
					x=v[jj][j];
					z=v[jj][i];
					v[jj][j]=x*c+z*s;
					v[jj][i]=z*c-x*s;
				}
				z=dpythag(f,h);
				w[j]=z;
				if (z) {
					z=1.0/z;
					c=f*z;
					s=h*z;
				}
				f=c*g+s*y;
				x=c*y-s*g;
				for (jj=1;jj<=m;jj++) {
					y=a[jj][j];
					z=a[jj][i];
					a[jj][j]=y*c+z*s;
					a[jj][i]=z*c-y*s;
				}
			}
			rv1[l]=0.0;
			rv1[k]=f;
			w[k]=x;
		}
	}
	free_dvector(rv1,1,n);
}


double dpythag(double a, double b)
/* computes (a^2+b^2)^0.5 without underflow or overflow */

{
	double absa,absb;
	absa=fabs(a);
	absb=fabs(b);
	if (absa > absb) return absa*sqrt(1.0+DSQR(absb/absa));
	else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+DSQR(absa/absb)));
}


double **SVDinverse(double **a, long m, long n){
/* inverts a matrix a[1.m][1.n]
by SVD to produce inv[1.n][1.m]
matrix not destoyed */

int i,j,k;
double wmax,wmin;
double **u,*w,**v,**t,**inv;
inv=dmatrix(1,n,1,m);
t=dmatrix(1,n,1,m);
u=dmatrix(1,m,1,n);
w=dvector(1,n);
v=dmatrix(1,n,1,n);
SVDcmp(a,m,n,u,w,v);

/* scan for eigenvalues relatively close to 0 */
wmax=0.0;
for(j=1;j<=n;j++)if(w[j]>wmax)wmax=w[j];
wmin=wmax*1.0e-12;
for(j=1;j<=n;j++)if(w[j]<wmin)
{printf("The eigenvalue number %d is %.5E, which is essentially zero\n",j,w[j]);
printf("Inverse does not exist.\n");
printf("If you object, study and change SVDinverse routine in matrixalgebra.c\n\n");
exit(1);}

for(i=1;i<=n;++i)for(j=1;j<=m;++j)t[i][j]=u[j][i]/w[i];
for(i=1;i<=n;++i)for(j=1;j<=m;++j){
	inv[i][j]=0.0;
	for(k=1;k<=m;++k)inv[i][j]=inv[i][j]+v[i][k]*t[k][j];
				  }
return inv;
}





void svdcmp(double **a, int m, int n, double w[], double **v)
{

/*
a = a * diag (w) * vT
Given matrix A as a[1..m][1..n] computes SVD as A=U*W*VT
matrix U is in a[1..m][1..n], so original is destroyed
singular values are in vector w[1..n]
matrix V (not its transpose) is in v[1..n][1..n]
*/



	int flag,i,its,j,jj,k,l,nm;
	double anorm,c,f,g,h,s,scale,x,y,z,*rv1;

	rv1=dvector(1,n);
	g=scale=anorm=0.0;
	for (i=1;i<=n;i++) {
		l=i+1;
		rv1[i]=scale*g;
		g=s=scale=0.0;
		if (i <= m) {
			for (k=i;k<=m;k++) scale += fabs(a[k][i]);
			if (scale) {
				for (k=i;k<=m;k++) {
					a[k][i] /= scale;
					s += a[k][i]*a[k][i];
				}
				f=a[i][i];
				g = -SIGN(sqrt(s),f);
				h=f*g-s;
				a[i][i]=f-g;
				for (j=l;j<=n;j++) {
					for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
					f=s/h;
					for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
				}
				for (k=i;k<=m;k++) a[k][i] *= scale;
			}
		}
		w[i]=scale *g;
		g=s=scale=0.0;
		if (i <= m && i != n) {
			for (k=l;k<=n;k++) scale += fabs(a[i][k]);
			if (scale) {
				for (k=l;k<=n;k++) {
					a[i][k] /= scale;
					s += a[i][k]*a[i][k];
				}
				f=a[i][l];
				g = -SIGN(sqrt(s),f);
				h=f*g-s;
				a[i][l]=f-g;
				for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
				for (j=l;j<=m;j++) {
					for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
					for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
				}
				for (k=l;k<=n;k++) a[i][k] *= scale;
			}
		}
		anorm=DMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
	}
	for (i=n;i>=1;i--) {
		if (i < n) {
			if (g) {
				for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g;
				for (j=l;j<=n;j++) {
					for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
					for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
				}
			}
			for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
		}
		v[i][i]=1.0;
		g=rv1[i];
		l=i;
	}
	for (i=IMIN(m,n);i>=1;i--) {
		l=i+1;
		g=w[i];
		for (j=l;j<=n;j++) a[i][j]=0.0;
		if (g) {
			g=1.0/g;
			for (j=l;j<=n;j++) {
				for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
				f=(s/a[i][i])*g;
				for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
			}
			for (j=i;j<=m;j++) a[j][i] *= g;
		} else for (j=i;j<=m;j++) a[j][i]=0.0;
		++a[i][i];
	}
	for (k=n;k>=1;k--) {
		for (its=1;its<=30;its++) {
			flag=1;
			for (l=k;l>=1;l--) {
				nm=l-1;
				if ((double)(fabs(rv1[l])+anorm) == anorm) {
					flag=0;
					break;
				}
				if ((double)(fabs(w[nm])+anorm) == anorm) break;
			}
			if (flag) {
				c=0.0;
				s=1.0;
				for (i=l;i<=k;i++) {
					f=s*rv1[i];
					rv1[i]=c*rv1[i];
					if ((double)(fabs(f)+anorm) == anorm) break;
					g=w[i];
					h=dpythag(f,g);
					w[i]=h;
					h=1.0/h;
					c=g*h;
					s = -f*h;
					for (j=1;j<=m;j++) {
						y=a[j][nm];
						z=a[j][i];
						a[j][nm]=y*c+z*s;
						a[j][i]=z*c-y*s;
					}
				}
			}
			z=w[k];
			if (l == k) {
				if (z < 0.0) {
					w[k] = -z;
					for (j=1;j<=n;j++) v[j][k] = -v[j][k];
				}
				break;
			}
			if (its == 30) nrerror("no convergence in 30 dsvdcmp iterations");
			x=w[l];
			nm=k-1;
			y=w[nm];
			g=rv1[nm];
			h=rv1[k];
			f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
			g=dpythag(f,1.0);
			f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
			c=s=1.0;
			for (j=l;j<=nm;j++) {
				i=j+1;
				g=rv1[i];
				y=w[i];
				h=s*g;
				g=c*g;
				z=dpythag(f,h);
				rv1[j]=z;
				c=f/z;
				s=h/z;
				f=x*c+g*s;
				g = g*c-x*s;
				h=y*s;
				y *= c;
				for (jj=1;jj<=n;jj++) {
					x=v[jj][j];
					z=v[jj][i];
					v[jj][j]=x*c+z*s;
					v[jj][i]=z*c-x*s;
				}
				z=dpythag(f,h);
				w[j]=z;
				if (z) {
					z=1.0/z;
					c=f*z;
					s=h*z;
				}
				f=c*g+s*y;
				x=c*y-s*g;
				for (jj=1;jj<=m;jj++) {
					y=a[jj][j];
					z=a[jj][i];
					a[jj][j]=y*c+z*s;
					a[jj][i]=z*c-y*s;
				}
			}
			rv1[l]=0.0;
			rv1[k]=f;
			w[k]=x;
		}
	}
	free_dvector(rv1,1,n);
}


void dsvbksb(double **u, double w[], double **v, int m, int n, double b[], double x[])
/*
Solves Ax=B, where A is specified by arrays

u[1.m][1.n]    w[1.n]    v[1.n][1.n]
b[1.m] is the right-hand side vector

results are in x[1.n]

no input is destroyed
*/
{
	int jj,j,i;
	double s,*tmp;

	tmp=dvector(1,n);
	for (j=1;j<=n;j++) {
		s=0.0;
		if (w[j]) {
			for (i=1;i<=m;i++) s += u[i][j]*b[i];
			s /= w[j];
		}
		tmp[j]=s;
	}
	for (j=1;j<=n;j++) {
		s=0.0;
		for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj];
		x[j]=s;
	}
	free_dvector(tmp,1,n);
}



double *SVDsolve(double **a, long m, long n, double *b){

/*
Solves system a*x=b, a[1..m][1..n], b[1..m]
returns x, x[1..n]

number of linear equations is m, number of unknowns is n

does not destroy originals
Good way to solve systems even singular
Uses SVD decomposition
*/

int j;
double *x,**u,*w,**v;
double wmax,wmin;

u=dmatrix(1,m,1,n);
w=dvector(1,n);
v=dmatrix(1,n,1,n);
x=dvector(1,n);

SVDcmp(a,m,n,u,w,v);

/* scan for eigenvalues relatively close to 0 */

wmax=0.0;
for(j=1;j<=n;j++)if(w[j]>wmax)wmax=w[j];
wmin=wmax*1.0e-12;
for(j=1;j<=n;j++)if(w[j]<wmin){
/*
printf("The eigenvalue number %d is %.5E, which is close to zero\n",j,w[j]);
printf("It is forced to be zero to improve performance.\n");
printf("If you object, study and change  SVDsolve routine in matrixalgebra.c\n\n");
*/
w[j]=0.0;}

dsvbksb(u,w,v,m,n,b,x);

return x;
}



void choldc(double **a, int n, double p[])

/*
Cholesky decomposition
APPLIED ONLY FOR positive definine (all eigenvalues positive, or v.A.v > 0 for all vectors v)
		 symmetrical (a_ij=a_ji) matrix a[1.n][1.n]
A=L.LT
upper triangle of A is needed and not modified. L is returned in lower triangle
diagonal elements of L are returned to p[1.n]
THIS FAILS if matrix is not positive definite
*/
{
	int i,j,k;
	double sum;

	for (i=1;i<=n;i++) {
		for (j=i;j<=n;j++) {
			for (sum=a[i][j],k=i-1;k>=1;k--) sum -= a[i][k]*a[j][k];
			if (i == j) {
				if (sum <= 0.0)
					nrerror("choldc failed");
				p[i]=sqrt(sum);
			} else a[j][i]=sum/p[i];
		}
	}
}


void cholsl(double **a, int n, double p[], double b[], double x[])

/*
Solves the set of equations A.x=b, where A is a positive-definite symmetric matrix.
Takes input from the choldc routine: a and p are exatly as they were output by it
b is right-hand side vector. Solutions are in x.
No input is destroyed
*/
{
	int i,k;
	double sum;

	for (i=1;i<=n;i++) {
		for (sum=b[i],k=i-1;k>=1;k--) sum -= a[i][k]*x[k];
		x[i]=sum/p[i];
	}
	for (i=n;i>=1;i--) {
		for (sum=x[i],k=i+1;k<=n;k++) sum -= a[k][i]*x[k];
		x[i]=sum/p[i];
	}
}

double  *CHsolve(double **a, long n, double *b){
/* 
solves A.x=b for x[1.n], where A is symmetric positive-definite a[1.n][1.n]
b[1.n] is a right-hand side vector 
input is not destroyed
*/

double *x,*p,**a1;
x=dvector(1,n);
p=dvector(1,n);
a1=mat_duplicate(a,1,n,1,n);

print_aa_matrix(a1,n);

choldc(a1,n,p);
cholsl(a1,n,p,b,x);
return x;
}



void tred2(double **a, int n, double d[], double e[])
/* reduces a real symmetric matrix a[1.n][1.n] to tridiagonal form
a is replaced by orthoginal matrix, necessary further
d[1..n] - diagonal elements of the tridiagonal matrix;
e[1..n] - off-diagonal elements with e[1]=0
*/

{
	int l,k,j,i;
	double scale,hh,h,g,f;

	for (i=n;i>=2;i--) {
		l=i-1;
		h=scale=0.0;
		if (l > 1) {
			for (k=1;k<=l;k++)
				scale += fabs(a[i][k]);
			if (scale == 0.0)
				e[i]=a[i][l];
			else {
				for (k=1;k<=l;k++) {
					a[i][k] /= scale;
					h += a[i][k]*a[i][k];
				}
				f=a[i][l];
				g=(f >= 0.0 ? -sqrt(h) : sqrt(h));
				e[i]=scale*g;
				h -= f*g;
				a[i][l]=f-g;
				f=0.0;
				for (j=1;j<=l;j++) {
					a[j][i]=a[i][j]/h;
					g=0.0;
					for (k=1;k<=j;k++)
						g += a[j][k]*a[i][k];
					for (k=j+1;k<=l;k++)
						g += a[k][j]*a[i][k];
					e[j]=g/h;
					f += e[j]*a[i][j];
				}
				hh=f/(h+h);
				for (j=1;j<=l;j++) {
					f=a[i][j];
					e[j]=g=e[j]-hh*f;
					for (k=1;k<=j;k++)
						a[j][k] -= (f*e[k]+g*a[i][k]);
				}
			}
		} else
			e[i]=a[i][l];
		d[i]=h;
	}
	d[1]=0.0;
	e[1]=0.0;
	/* Contents of this loop can be omitted if eigenvectors not
			wanted except for statement d[i]=a[i][i]; */
	for (i=1;i<=n;i++) {
		l=i-1;
		if (d[i]) {
			for (j=1;j<=l;j++) {
				g=0.0;
				for (k=1;k<=l;k++)
					g += a[i][k]*a[k][j];
				for (k=1;k<=l;k++)
					a[k][j] -= g*a[k][i];
			}
		}
		d[i]=a[i][i];
		a[i][i]=1.0;
		for (j=1;j<=l;j++) a[j][i]=a[i][j]=0.0;
	}
}

void tqli(double d[], double e[], int n, double **z)
/*
determines eigenvalues and eigenvectors of real symmetric tridiagonal matrix:
on input d[1..n] dagonal elements of tridiagonal matrix; on output - eigenvalues;
e[1..n]
the k-th column of z returns eigenvector corresponding to d[k];
on input z is matrix, outputed by tred2 
*/

{
	int m,l,iter,i,k;
	double  s,r,p,g,f,dd,c,b;

	for (i=2;i<=n;i++) e[i-1]=e[i];
	e[n]=0.0;
	for (l=1;l<=n;l++) {
		iter=0;
		do {
			for (m=l;m<=n-1;m++) {
				dd=fabs(d[m])+fabs(d[m+1]);
				if ((double)(fabs(e[m])+dd) == dd) break; 
			/*	if (fabs(e[m]) <1.0e-14) break;  */

			}
			if (m != l) {
				if (iter++ == 30) nrerror("Too many iterations in tqli");
				g=(d[l+1]-d[l])/(2.0*e[l]);
				r=dpythag(g,1.0);
				g=d[m]-d[l]+e[l]/(g+SIGN(r,g));
				s=c=1.0;
				p=0.0;
				for (i=m-1;i>=l;i--) {
					f=s*e[i];
					b=c*e[i];
					e[i+1]=(r=dpythag(f,g));
					if (r == 0.0) {
						d[i+1] -= p;
						e[m]=0.0;
						break;
					}
					s=f/r;
					c=g/r;
					g=d[i+1]-p;
					r=(d[i]-g)*s+2.0*c*b;
					d[i+1]=g+(p=s*r);
					g=c*r-b;
					for (k=1;k<=n;k++) {
						f=z[k][i+1];
						z[k][i+1]=s*z[k][i]+c*f;
						z[k][i]=c*z[k][i]-s*f;
					}
				}
				if (r == 0.0 && i >= l) continue;
				d[l] -= p;
				e[l]=g;
				e[m]=0.0;
			}
		} while (m != l);
	}
}


void eigensystem_symm(double **a, long n, double *eval, double **evec, double *helper)
/* solves eigensystem for symmetric real matrix a[1.n][1.n]
stores eigenvalues in eval[1..n] and correspoding vectors as columns in evec[1..n][1..n]
helper[1..n] is a storage vector to save time for its allocation
*/
{
int i,j;
for(i=1;i<=n;++i)for(j=1;j<=n;++j)evec[i][j]=a[i][j];
tred2(evec,n,eval,helper);
tqli(eval,helper,n,evec);
}



void balanc(double **a, int n)
/* balances a[1.n][1.n]
needed for eigenvalue calculation 
a destryed on output */
{
	int last,j,i;
	double s,r,g,f,c,sqrdx;

	sqrdx=RADIX*RADIX;
	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/RADIX;
				f=1.0;
				s=c+r;
				while (c<g) {
					f *= RADIX;
					c *= sqrdx;
				}
				g=r*RADIX;
				while (c>g) {
					f /= RADIX;
					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;
				}
			}
		}
	}
}




void elmhes(double **a, int n)
/* reduces balanced matrix a[1.n][1.n] to Hessenberg form */
{
	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++) SWAPelmhes(a[i][j],a[m][j])
			for (j=1;j<=n;j++) SWAPelmhes(a[j][i],a[j][m])
		}
		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];
				}
			}
		}
	}
}


void hqr(double **a, int n, double wr[], double wi[])
/* 
calculates eigenvalues of an upper Hessenberg matrix a[1.n][1.n] 
real and imaginary parts of eigenvalues are in wr and 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+SIGN(z,p);
						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=SIGN(sqrt(p*p+q*q+r*r),p)) != 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);
	}
}



void ldbalanc(long double **a, int n)
/* balances a[1.n][1.n]
needed for eigenvalue calculation 
a destryed on output */
{
	int last,j,i;
	long double s,r,g,f,c,sqrdx;

	sqrdx=RADIX*RADIX;
	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 += fabslmy(a[j][i]);
					r += fabslmy(a[i][j]);
				}
			if (c && r) {
				g=r/RADIX;
				f=1.0;
				s=c+r;
				while (c<g) {
					f *= RADIX;
					c *= sqrdx;
				}
				g=r*RADIX;
				while (c>g) {
					f /= RADIX;
					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;
				}
			}
		}
	}
}




void ldelmhes(long double **a, int n)
/* reduces balanced matrix a[1.n][1.n] to Hessenberg form */
{
	int m,j,i;
	long double y,x;

	for (m=2;m<n;m++) {
		x=0.0;
		i=m;
		for (j=m;j<=n;j++) {
			if (fabslmy(a[j][m-1]) > fabslmy(x)) {
				x=a[j][m-1];
				i=j;
			}
		}
		if (i != m) {
			for (j=m-1;j<=n;j++) SWAPelmhes(a[i][j],a[m][j])
			for (j=1;j<=n;j++) SWAPelmhes(a[j][i],a[j][m])
		}
		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];
				}
			}
		}
	}
}



void ldhqr(long double **a, int n, long double wr[], long double wi[])


/* 
calculates eigenvalues of an upper Hessenberg matrix a[1.n][1.n] 
real and imaginary parts of eigenvalues are in wr and wi
*/

{
	int nn,m,l,k,j,its,i,mmin;
	long double z,y,x,w,v,u,t,s,r,q,p,anorm;

	anorm=fabslmy(a[1][1]);
	for (i=2;i<=n;i++)
		for (j=(i-1);j<=n;j++)
			anorm += fabslmy(a[i][j]);
	nn=n;
	t=0.0;
	while (nn >= 1) {
		its=0;
		do {
			for (l=nn;l>=2;l--) {
				s=fabslmy(a[l-1][l-1])+fabslmy(a[l][l]);
				if (s == 0.0) s=anorm;
				if ((long double)(fabslmy(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;
/* should be sqrtl, but not compiles */	z=sqrt(fabslmy(q));
					x += t;
					if (q >= 0.0) {
						z=p+SIGNL(z,p);
						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 ldhqr");
					if (its == 10 || its == 20) {
						t += x;
						for (i=1;i<=nn;i++) a[i][i] -= x;
						s=fabslmy(a[nn][nn-1])+fabslmy(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=fabslmy(p)+fabslmy(q)+fabslmy(r);
						p /= s;
						q /= s;
						r /= s;
						if (m == l) break;
						u=fabslmy(a[m][m-1])*(fabslmy(q)+fabslmy(r));
						v=fabslmy(p)*(fabslmy(a[m-1][m-1])+fabslmy(z)+fabslmy(a[m+1][m+1]));
						if ((long 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 ((long double)(x=fabslmy(p)+fabslmy(q)+fabslmy(r)) != 0.0) {
								p /= x;
								q /= x;
								r /= x;
							}
						}
/* should be sqrtl, but not compiles */	if ((long double)(s=SIGNL(sqrt(p*p+q*q+r*r),p)) != 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);
	}
}




void eigenvalues(double **a,long n, double *eiR, double *eiI){
/* calculates eigenvalues of a square matrix a[1.n][1.n] */

double **m;
m=mat_duplicate(a,1,n,1,n);
balanc(m,n);
elmhes(m,n);
hqr(m,n,eiR,eiI);
}


void ldeigenvalues(long double **a,long n, long double *eiR, long double *eiI){
/* calculates eigenvalues of a square matrix a[1.n][1.n] */
int i,j;
long double **m;
m=ldmatrix(1,n,1,n);
for(i=1;i<=n;++i)for(j=1;j<=n;++j)m[i][j]=a[i][j];
ldbalanc(m,n);
ldelmhes(m,n);
ldhqr(m,n,eiR,eiI);
}


void mat_prepped(double **a, long n, double ***mp, double *evR, double *evI){
/* prepares a tensor ***mp[1.n][1.n][1.n] for taking any function of a 
symmetric real matrix a[1.n][1.n] 
eigenvalues are in evR[1.n]  (real part) and evI[1.n] (imaginary) and 
matrix, corresponding to ev?[k] is in mp[k][1.n][1.n] */

int i,i1,j,k;
double **m,diff,f=0.9999999;
m=dmatrix(1,n,1,n);

eigenvalues(a,n,evR,evI);

for(k=1;k<=n;++k){
	for(i=1;i<=n;++i)for(j=1;j<=n;++j){
		if(i==j)mp[k][i][j]=1.0;
		else mp[k][i][j]=0.0;	  }

	for(i=1;i<=n;++i){
		if(i==k)continue;
		diff=evR[k]-evR[i];
	if(fabs(diff)<=1.0e-10)nrerror("not all eigenvalues are different in mat_prepped");

	for(i1=1;i1<=n;++i1)for(j=1;j<=n;++j){
		if(i1==j)m[i1][j]=f*(a[i1][j]-evR[i])/diff;
		else m[i1][j]=f*a[i1][j]/diff;}
	mp[k]=matXmat(mp[k],1,n,1,n,m,1,n,1,n);
			 }
		 }
for(i=1;i<=20;++i)for(k=1;k<=n;++k)mp[k]=dXmat(1.0/f,mp[k],1,20,1,20);
}


void ldmat_prepped(long double **a, long n, long double ***mp, long double *evR, long double *evI){
/* prepares a tensor ***mp[1.n][1.n][1.n] for taking any function of a 
symmetric real matrix a[1.n][1.n] 
eigenvalues are in evR[1.n]  (real part) and evI[1.n] (imaginary) and 
matrix, corresponding to ev?[k] is in mp[k][1.n][1.n] */

int i,i1,j,k;
long double **m,diff;
double **ad;
/* double *evId,*evRd; */
m=ldmatrix(1,n,1,n);
ad=dmatrix(1,n,1,n);
/*
evRd=dvector(1,n);
evId=dvector(1,n);
*/
for(i=1;i<=n;++i)for(j=1;j<=n;++j)ad[i][j]=a[i][j];

/* eigenvalues(ad,n,evRd,evId); */

ldeigenvalues(a,n,evR,evI);

/*
for(i=1;i<=n;++i){evR[i]=evRd[i];evI[i]=evId[i];}
*/
printf("done eigenvalues\n");
for(k=1;k<=n;++k){
	for(i=1;i<=n;++i)for(j=1;j<=n;++j){
		if(i==j)mp[k][i][j]=1.0;
		else mp[k][i][j]=0.0;	  }

	for(i=1;i<=n;++i){
		if(i==k)continue;
		diff=evR[k]-evR[i];
	if(fabslmy(diff)<=1.0e-10)nrerror("not all eigenvalues are different in mat_prepped");

	for(i1=1;i1<=n;++i1)for(j=1;j<=n;++j){
		if(i1==j)m[i1][j]=(a[i1][j]-evR[i])/diff;
		else m[i1][j]=a[i1][j]/diff;}
	mp[k]=ldmatXldmat(mp[k],1,n,1,n,m,1,n,1,n);
			 }
		 }
/* for(i=1;i<=20;++i)for(k=1;k<=n;++k)mp[k]=lddXldmat(1.0/f,mp[k],1,20,1,20); */
}


void mat_EXP_prepped(double ***mp, long n, double *eiR, double **m, \
double p){ 
/* calculates m[1.n][1.n]=exp(A[1.n][1.n]*p) from
prepped matrix with real eignevalues eiR[1.k.n] and n corresponding
matrices mp[1.k.n][1.n][1.n]*/
int i,j,k;
double d;
d=exp(p*eiR[1]);
for(i=1;i<=n;++i)for(j=1;j<=n;++j)m[i][j]=d*mp[1][i][j];
for(k=2;k<=n;++k){
	d=exp(p*eiR[k]);
	for(i=1;i<=n;++i)for(j=1;j<=n;++j)m[i][j]=m[i][j]+d*mp[k][i][j];
		 }
}
void mat_EXPDI_prepped(double ***mp, long n, double *eiR, double **m, \
double p){
/* calculates m[1.n][1.n]=?exp(A[1.n][1.n]*p) from
prepped matrix with real eignevalues eiR[1.k.n] and n corresponding
matrices mp[1.k.n][1.n][1.n]*/
int i,j,k;
double d;
d=1.0-p*eiR[1];d=1.0/d;
for(i=1;i<=n;++i)for(j=1;j<=n;++j)m[i][j]=d*mp[1][i][j];
for(k=2;k<=n;++k){
        d=1.0-p*eiR[k];d=1.0/d;
        for(i=1;i<=n;++i)for(j=1;j<=n;++j)m[i][j]=m[i][j]+d*mp[k][i][j];
                 }
}


/*---------------------------------------------------*/

/* 
work with NCBI-style formatted amino-acid integer matrices that look like this

#  Matrix made by matblas from blosum62.iij
#  * column uses minimum score
#  BLOSUM Clustered Scoring Matrix in 1/2 Bit Units
#  Blocks Database = /data/blocks_5.0/blocks.dat
#  Cluster Percentage: >= 62
#  Entropy =   0.6979, Expected =  -0.5209
   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  Z  X  *
A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0 -2 -1  0 -4 
R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3 -1  0 -1 -4 
N -2  0  6  1 -3  0  0  0  1 -3 -3  0 -2 -3 -2  1  0 -4 -2 -3  3  0 -1 -4 
D -2 -2  1  6 -3  0  2 -1 -1 -3 -4 -1 -3 -3 -1  0 -1 -4 -3 -3  4  1 -1 -4 
C  0 -3 -3 -3  9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -4 
Q -1  1  0  0 -3  5  2 -2  0 -3 -2  1  0 -3 -1  0 -1 -2 -1 -2  0  3 -1 -4 
E -1  0  0  2 -4  2  5 -2  0 -3 -3  1 -2 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4 
G  0 -2  0 -1 -3 -2 -2  6 -2 -4 -4 -2 -3 -3 -2  0 -2 -2 -3 -3 -1 -2 -1 -4 
H -2  0  1 -1 -3  0  0 -2  8 -3 -3 -1 -2 -1 -2 -1 -2 -2  2 -3  0  0 -1 -4 
I -1 -3 -3 -3 -1 -3 -3 -4 -3  4  2 -3  1  0 -3 -2 -1 -3 -1  3 -3 -3 -1 -4 
L -1 -2 -3 -4 -1 -2 -3 -4 -3  2  4 -2  2  0 -3 -2 -1 -2 -1  1 -4 -3 -1 -4 
K -1  2  0 -1 -3  1  1 -2 -1 -3 -2  5 -1 -3 -1  0 -1 -3 -2 -2  0  1 -1 -4 
M -1 -1 -2 -3 -1  0 -2 -3 -2  1  2 -1  5  0 -2 -1 -1 -1 -1  1 -3 -1 -1 -4 
F -2 -3 -3 -3 -2 -3 -3 -3 -1  0  0 -3  0  6 -4 -2 -2  1  3 -1 -3 -3 -1 -4 
P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4  7 -1 -1 -4 -3 -2 -2 -1 -2 -4 
S  1 -1  1  0 -1  0  0  0 -1 -2 -2  0 -1 -2 -1  4  1 -3 -2 -2  0  0  0 -4 
T  0 -1  0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1  1  5 -2 -2  0 -1 -1  0 -4 
W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1  1 -4 -3 -2 11  2 -3 -4 -3 -2 -4 
Y -2 -2 -2 -3 -2 -1 -2 -3  2 -1 -1 -2 -1  3 -3 -2 -2  2  7 -1 -3 -2 -1 -4 
V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4 -3 -2 -1 -4 
B -2 -1  3  4 -3  0  1 -1  0 -3 -4  0 -3 -3 -2  0 -1 -4 -3 -3  4  1 -1 -4 
Z -1  0  0  1 -3  3  4 -2  0 -3 -3  1 -1 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4 
X  0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2  0  0 -2 -1 -1 -1 -1 -1 -4 
* -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4  1 

everyting in the line after #  at the beginning of the file is ignored
*/


int  **read_aa_imatrix(FILE *fmat){

/* read matrix from file *fmat */

int i,ri,rj,c,flag,j,t;
int col[31],row[31];
char stri[11];
int **mat;

mat=imatrix(0,25,0,25);
for(i=0;i<=25;++i)for(j=0;j<=25;++j)mat[i][j]=0;

i=0;
ri=0;
rj=0;

flag=0;


while( (c=getc(fmat)) != EOF){

if(flag==0 && c=='#'){flag=-1;continue;}
else if(flag==-1 && c=='\n'){flag=0;continue;}
else if(flag==-1){continue;}
else if(flag==0 && c==' '){flag=1;continue;}
else if(flag==1 && c=='\n'){flag=0;continue;}
else if(flag==1 && c==' '){continue;}
else if(flag==1){
		++i;
		if(i>=25){nrerror("matrix has more than 25 columns: FATAL");exit(0);}
		col[i]=am2numBZX(c);
		continue;
		}
else if(flag==0 && c!=' ' && c!='#'){
		ri=0;
		++rj;
		if(rj>=25){nrerror("matrix has more than 25 rows: FATAL");exit(0);}
		row[rj]=am2numBZX(c);
		flag=2;
		continue;
		}
else if (flag==2 && c==' '){for(i=0;i<=10;++i){stri[i]=' ';}j=0;continue;}
else if (flag==2 && c=='\n'){flag=0;continue;}
else if (flag==2){flag=3;stri[++j]=c;if(j>10){nrerror("string too long:FATAL");exit(0);}continue;} 
else if (flag==3 && c==' ' || flag==3 && c=='\n'){
		/*	printf("%s|\n",stri); */
			j=0;
			++ri;
			t=atoi(stri);
			mat[row[rj]][col[ri]]=t;
/*			printf("%2d %2d %2d\n",t,rj,ri); */
			if (c=='\n')flag=0;else flag=2;
			continue;
			}
else if (flag==3){stri[++j]=c;continue;}

}

return mat;
}


double  *read_aa_vector(FILE *fvec){

/* read vector of floats from file *fvec in a format:

# Dayhoff amino acid equilibrium frequencies
W   0.010493678893
F   0.039771658566
Y   0.029916218247
M   0.014752761601
L   0.085357443759
I   0.036885613389
V   0.064717601578
A   0.087126757719
C   0.033473507356
G   0.088612083614
P   0.050679914837
T   0.058541664191
S   0.069577076005
N   0.040432285538
Q   0.038255436979
D   0.046871932808
E   0.049530476922
H   0.033618414315
R   0.040903890308
K   0.080481583376

*/

int i;
char *s;
double *vec;

s=cvector(0,100);
vec=dvector(0,25);
for(i=0;i<=25;++i)vec[i]=0.0;

for(;;) {
        s=fgets(s,100,fvec);
	if(s==NULL) break;
	if(s[0]=='#')continue;
	if(s[0]=='\n')continue;
	if(s[0]==' ')continue;
	i=am2numBZX(s[0]);
	s[0]=' ';
/*	i = sscanf(s,"%c %f",c,d); */
	vec[i]=atof(s);
	}
return vec;
}

int  *read_aa_ivector(FILE *fvec){

/* read vector of floats from file *fvec 
*/

int i;
char *s;
int *vec;

s=cvector(0,100);
vec=ivector(0,25);
for(i=0;i<=25;++i)vec[i]=0;

for(;;) {
        s=fgets(s,100,fvec);
	if(s==NULL) break;
	if(s[0]=='#')continue;
	if(s[0]=='\n')continue;
	if(s[0]==' ')continue;
	i=am2numBZX(s[0]);
	s[0]=' ';
/*	i = sscanf(s,"%c %f",c,d); */
	vec[i]=atoi(s);
	}
return vec;
}


void print_aa_imatrix(int **mat, int n){
/* prints inetger amino acid matrix mat[1.n][1.n], indices mean amino acids */

int i,j;
for(i=0;i<=n;++i){
if (i==0){
	printf(" ");
	for(j=1;j<=n;++j)printf("%3c",am[j]);
	printf("\n");
	continue;
	}
for(j=0;j<=n;++j){
if(j==0){printf("%1c",am[i]);continue;}
printf("%3d",mat[i][j]);}
printf("\n");
	          }
}

void print_aa_imatrix5(int **mat, int n){
/* prints inetger amino acid matrix mat[1.n][1.n], indices mean amino acids */

int i,j;
for(i=0;i<=n;++i){
if (i==0){
	printf(" ");
	for(j=1;j<=n;++j)printf("%5c",am[j]);
	printf("\n");
	continue;
	}
for(j=0;j<=n;++j){
if(j==0){printf("%1c",am[i]);continue;}
printf("%5d",mat[i][j]);}
printf("\n");
	          }
}
void fprint_aa_imatrix5(FILE *fo, int **mat, int n){
/* prints inetger amino acid matrix mat[1.n][1.n], indices mean amino acids */

int i,j;
for(i=0;i<=n;++i){
if (i==0){
	fprintf(fo," ");
	for(j=1;j<=n;++j)fprintf(fo,"%5c",am[j]);
	fprintf(fo,"\n");
	continue;
	}
for(j=0;j<=n;++j){
if(j==0){fprintf(fo,"%1c",am[i]);continue;}
fprintf(fo,"%5d",mat[i][j]);}
fprintf(fo,"\n");
	          }
}
void print_aa_imatrix6(int **mat, int n){
/* prints inetger amino acid matrix mat[1.n][1.n], indices mean amino acids */

int i,j;
for(i=0;i<=n;++i){
if (i==0){
	printf(" ");
	for(j=1;j<=n;++j)printf("%6c",am[j]);
	printf("\n");
	continue;
	}
for(j=0;j<=n;++j){
if(j==0){printf("%1c",am[i]);continue;}
printf("%6d",mat[i][j]);}
printf("\n");
	          }
}

void print_aa_imatrix_dayhoff(int **mat, int n){
/* prints inetger amino acid matrix mat[1.n][1.n], indices mean amino acids */

int i,j;
for(i=0;i<=n;++i){
if (i==0){
	printf(" ");
	for(j=1;j<=n;++j)printf("%4c",am_dayhoff[j]);
	printf("\n");
	continue;
	}
for(j=0;j<=n;++j){
if(j==0){printf("%1c",am_dayhoff[i]);continue;}
printf("%4d",mat[i][j]);}
printf("\n");
	          }
}

void print_aa_vector(double *vec, int n){
/* printfs amino acid vector vec[1.n][1.n] */

int i;
for(i=1;i<=n;++i)printf("%1c %14.12f\n",am[i],vec[i]);
}

void print_aa_ivector(int *vec, int n){
/* printfs amino acid vector vec[1.n][1.n] */

int i;
for(i=1;i<=n;++i)printf("%1c %10d\n",am[i],vec[i]);
}


void print_aa_matrix(double **mat, int n){
/* prints amino acid matrix mat[1.n][1.n], indices mean amino acids */

int i,j;
for(i=0;i<=n;++i){
if (i==0){	
	printf(" ");
	for(j=1;j<=n;++j)printf("%7c",am[j]);
	printf("\n");
	continue;
	}
for(j=0;j<=n;++j){
if(j==0){printf("%1c",am[i]);continue;}
printf("%7.4f",mat[i][j]);}
printf("\n");
	          }
}
void print_aa_ldmatrix(long double **mat, int n){
/* prints amino acid matrix mat[1.n][1.n], indices mean amino acids */

int i,j;
for(i=0;i<=n;++i){
if (i==0){	
	printf(" ");
	for(j=1;j<=n;++j)printf("%7c",am[j]);
	printf("\n");
	continue;
	}
for(j=0;j<=n;++j){
if(j==0){printf("%1c",am[i]);continue;}
printf("%7.4f",(double)mat[i][j]);}
printf("\n");
	          }
}


void print_matrix(double **m,long n1,long n2,char *txt){
int i,j;

printf("%s\n",txt);
for(i=1;i<=n1;i++){
  for(j=1;j<=n2;j++)printf("%10.4f",m[i][j]);
  printf("\n");
		  }
}
