#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define NRANSI
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <stdarg.h>

#define NR_END 1
#define FREE_ARG char*

/* nrutil.h **/
static double sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0f ? 0.0f : sqrarg*sqrarg)

#define SIGN(a,b) ((b) >= 0.0f ? (double)fabs(a) : -(double)fabs(a))

static double maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
        (maxarg1) : (maxarg2))

static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
        (iminarg1) : (iminarg2))

/* end nrutil.h */

/** nr.h **/
static void  svdcmp(double **a, int m, int n, double w[], double **v);
static void  svbksb(double **u, double w[], double **v, int m, int n, double b[], double x[]);
static void ludcmp(double **a, int n, int *indx, double *d);
static void lubksb(double **a, int n, int *indx, double b[]);


static double **convert_matrix(double *a, long nrl, long nrh, long ncl, long nch);
static double *vector(long nl, long nh);
static void  free_vector(double *v, long nl, long nh);
static double **matrix(long nrl, long nrh, long ncl, long nch);
static void  free_matrix(double **m, long nrl, long nrh, long ncl, long nch);
static void  free_convert_matrix(double **b, long nrl, long nrh, long ncl, long nch);
static int *ivector(long nl, long nh);
static void free_ivector(int *v, long nl, long nh);

/** end nr **/

/* Funcoes utilitarias */
static void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
	fprintf(stderr,"Numerical Recipes run-time error...\n");
	fprintf(stderr,"%s\n",error_text);
	fprintf(stderr,"...now exiting to system...\n");
	exit(1);
}

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

static double **convert_matrix(double *a, long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix m[nrl..nrh][ncl..nch] that points to the matrix
declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
and ncol=nch-ncl+1. The routine should be called with the address
&a[0][0] as the first argument. */
{
	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
	double **m;

	/* allocate pointers to rows */
	m=(double **) malloc((size_t) ((nrow+NR_END)*sizeof(double*)));
	if (!m) nrerror("allocation failure in convert_matrix()");
	m += NR_END;
	m -= nrl;

	/* set pointers to rows */
	m[nrl]=a-ncl;
	for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
	/* return pointer to array of pointers to rows */
	return m;
}


static double *vector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
	double *v;

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


static int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
	int *v;

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

static void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

static void free_vector(double *v, long nl, long nh)
/* free a double vector allocated with vector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

static double **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
	double **m;

	/* allocate pointers to rows */
	m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m += NR_END;
	m -= nrl;

	/* allocate rows and set pointers to them */
	m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
	m[nrl] += NR_END;
	m[nrl] -= ncl;

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

	/* return pointer to array of pointers to rows */
	return m;
}

static void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by matrix() */
{
	free((FREE_ARG) (m[nrl]+ncl-NR_END));
	free((FREE_ARG) (m+nrl-NR_END));
}


static void free_convert_matrix(double **b, long nrl, long nrh, long ncl, long nch)
/* free a matrix allocated by convert_matrix() */
{
	free((FREE_ARG) (b+nrl-NR_END));
}



static void svdcmp(double **a, int m, int n, double w[], double **v)
{
	double pythag(double a, double b);
	int flag,i,its,j,jj,k,l,nm;
	double anorm,c,f,g,h,s,scale,x,y,z,*rv1;

	rv1=vector(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 += (double) 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 += (double)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=FMAX(anorm,((double)fabs(w[i])+(double)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=(double)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=pythag(f,g);
					w[i]=h;
					h=(double)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 svdcmp 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.0f*h*y);
			g=pythag(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=pythag(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=pythag(f,h);
				w[j]=z;
				if (z) {
					z=1.0f/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_vector(rv1,1,n);
}


static void svbksb(double **u, double w[], double **v, int m, int n, double b[], double x[])
{
	int jj,j,i;
	double s,*tmp;

	tmp=vector(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_vector(tmp,1,n);
}


#define TINY (double)1.0e-20;

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

	vv=vector(1,n);
	*d=1.0;
	for (i=1;i<=n;i++) {
		big=0.0;
		for (j=1;j<=n;j++)
			if ((temp=(double)fabs(a[i][j])) > big) big=temp;
		if (big == 0.0) nrerror("Singular matrix in routine LUDCMP");
		vv[i]=(double)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=(double)(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.0f) a[j][j]=TINY;
		if (j != n) {
			dum=(double)1.0/(a[j][j]);
			for (i=j+1;i<=n;i++) a[i][j] *= dum;
		}
	}
	free_vector(vv,1,n);
}

#undef TINY

static void lubksb(double **a,int n,int *indx,double *b)
{
	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];
	}
}


#undef NRANSI
/* (C) Copr. 1986-92 Numerical Recipes Software "!15L1. */

/*******************************************************/
/* * Interface para as funcoes do Numerical Recipes * *
 * * Marcelo Gattass - 28 de setembro de 2003        * */
/*******************************************************/


/*
 * Acha uma solucao nao trivial x para o Sistema Ax=0
 * A mxn, com m>n e rank(A)< n (preferencialmente = n-1).
 *
 * O valor de retorno indica erro. Ele e' a razao do menor 
 * para o maior valor singular de A. Se o rank(A) for menor 
 * que n, este valor deveria ser zero (~=1.0e-6 para esta
 * implemtacao com precisao simples).
 *  
 * Adaptada de Numerical Recipes Software 
 */ 
double nrSVDAx0(double **a, int m, int n, double *x)
{
	double **u=matrix(1,m,1,n);
	double *w=vector(1,n);
	double **v=matrix(1,n,1,n);

	double wmax,wmin;
	int i,j,jmin,nzeros=0,error=0;
		
	/* copia a matrix original em u */
	for (i=1;i<=m;i++)
		for (j=1;j<=n;j++) {
				u[i][j]=a[i-1][j-1];
			}

    /* perform decomposition */
	svdcmp(u,m,n,w,v);

	/* testa se a tem solucao nao trivial */
	wmax=w[1]; wmin = w[1];
	for (j=2;j<=n;j++){
		    if (w[j] < wmin)
			{
				wmin=w[j];
				jmin =j;
			}
			if (w[j] > wmax) wmax=w[j];
	}

	/* copia a coluna de v correspondente ao menor valor singular para x */
	for (j=1;j<=n;j++)
        x[j-1] = v[j][jmin];

	free_matrix(v,1,n,1,n);
	free_matrix(u,1,m,1,n);
	free_vector(w,1,n);

	return (wmin/wmax);
}

/*
 * Acha uma solucao x para o Sistema Ax=b
 * A mxn, com m>n.
 *
 * Adaptada de Numerical Recipes Software 
 */ 
void nrSVDAxb(double **a, int m, int n, double *b, double *x)
{
	double **u=matrix(1,m,1,n);
	double  *w=vector(1,n);
	double **v=matrix(1,n,1,n);
	
	double *x1=x-1;      /* mudando a origem para compatibilizar com a notacao do Livro */
	double *b1=b-1;
    
	double wmax,wmin;
	int i,j;
	
	/* copia a matrix original em u */
	for (i=1;i<=m;i++){
		for (j=1;j<=n;j++) {
				u[i][j]=a[i-1][j-1];
			}
	}

    /* perform decomposition */
	svdcmp(u,m,n,w,v);

	/* Calcula o maior valor singular */
	wmax=w[1];
	for (j=2;j<=n;j++){
		if (w[j] > wmax) wmax=w[j];
	}
	
	/* calcula o valor minimo que um valor singular pode ter */
	wmin=(double)1.e-6*wmax;

	/* zera os valores singulares que estao abaixo do minimo */
	for (j=1;j<=n;j++){
		if (w[j]<wmin) w[j]=0.0f;
	}

	svbksb(u,w,v,m,n,b1,x1);

    free_matrix(v,1,n,1,n);
	free_matrix(u,1,m,1,n);
	free_vector(w,1,n);
	return;
}


/*
 * Computa as matrizes [U],[W] e [V] da decomposicao singular da matriz 
 * [A] = [U][W][V]T
 * [A] mxn, [U]mxm (colunas ortogonais), [W]nxn (diagonal) e [V] nxn (ortogonal).
 *
 * Adaptada de Numerical Recipes Software 
 */ 
void nrSVD(double **a, int m, int n, double **u, double *w, double **v)
{
	double **u1=matrix(1,m,1,n);
	double **v1=matrix(1,n,1,n);
	double  *w1=w-1;

	int i,j;
		
	/* copia a matrix original em u */
	for (i=1;i<=m;i++) {
		for (j=1;j<=n;j++) {
				u1[i][j]=a[i-1][j-1];
			}
	}

    /* perform decomposition */
	svdcmp(u1,m,n,w1,v1);

	/* copia a matrix u1 em u */
	for (i=1;i<=m;i++) {
		for (j=1;j<=n;j++) {
				u[i-1][j-1]=u1[i][j];
			}
	}

	/* copia a matrix v1 em v */
	for (i=1;i<=n;i++) {
		for (j=1;j<=n;j++) {
				v[i-1][j-1]=v1[i][j];
			}
	}

	free_matrix(v1,1,n,1,n);
	free_matrix(u1,1,m,1,n);
}

/*
 * Acha uma solucao x para o Sistema Ax=b
 * A nxn, via decomposicao LU.
 *
 * Adaptada de Numerical Recipes Software 
 */ 
void nrLUAxb(double **a, int n, double *b, double *x)
{
	double **lu=matrix(1,n,1,n);
	int    *indx=ivector(1,n);
	double  d;
	
	double *x1=x-1;      /* mudando a origem para compatibilizar com a notacao do Livro */
    
	int i,j;
	
	/* copia a matrix original em lu */
	for (i=1;i<=n;i++){
		for (j=1;j<=n;j++) {
				lu[i][j]=a[i-1][j-1];
			}
	}

    /* decompoe */
    ludcmp(lu,n,indx,&d);

	/* copia o vetor b em x */
	for (j=0;j<n;j++)
		x[j]=b[j];


	/* retrosubstituicao */
	lubksb(lu, n, indx, x1);

    free_matrix(lu,1,n,1,n);
	free_ivector(indx,1,n);
	return;
}

/*
 * Computes the nxn matrix [X] such that [A][X]=[B]
 * If [B] is set to be the identity matrix, then [X]
 * is the inverse of [A].  (uses LU decomposition).
 *
 * Adaptada de Numerical Recipes Software 
 */ 
void nrLUAXB(double **a, int n, double **b, double **x)
{
	double **lu=matrix(1,n,1,n); 
	double *tmp=vector(1,n);
	int    *indx=ivector(1,n);
	double  d;
	
	int i,j;
	
	/* copia a matrix original em lu */
	for (i=1;i<=n;i++){
		for (j=1;j<=n;j++) {
				lu[i][j]=a[i-1][j-1];
			}
	}

    /* LU decomposition */
    ludcmp(lu,n,indx,&d);

	/* copy each column of [b] into {tmp}, solve and copy into [x] */
	for (j=0;j<n;j++) {
		for (i=0;i<n;i++) tmp[i+1]=b[i][j];
        lubksb(lu, n, indx, tmp);
		for (i=0;i<n-1;i++) x[i][j]=tmp[i+1];
	}


    free_matrix(lu,1,n,1,n);
	free_vector(tmp,1,n);
	free_ivector(indx,1,n);
	return;
}




/*
 * Simple algebra
 */

void nrAb(double **a, double *b, int n, int m, double *x)
{
	int i,j;
	for (i=0;i<n;i++)
	{
		x[i]=0.0;
		for (j=0;j<m;j++)
			x[i]+=a[i][j]*b[j];
	}
}

/* Dimensions: [A]=mxn, {b}=m and {x}=n. */
void nrAtb(double **a, double *b, int m, int n, double *x)
{
	int i,j;
	for (i=0;i<m;i++)
	{
		x[i]=0.0;
		for (j=0;j<n;j++)
			x[i]+=a[j][i]*b[j];
	}
}


void nrAB(double **a, double **b, int n, int p, int m, double **x)
{
	int i,j,k;
	for (i=0;i<n;i++) {
		for (j=0;j<m;j++)
		{
			x[i][j]=0.0;
			for (k=0;k<p;k++)
                x[i][j]+=a[i][k]*b[k][j];
		}
	}
}

/* Dimensions: [A]=nxp, [B]=mxp and [X]=nxm. */
void nrABt(double **a, double **b, int n, int p, int m, double **x)
{
	int i,j,k;
	for (i=0;i<n;i++) {
		for (j=0;j<m;j++)
		{
			x[i][j]=0.0;
			for (k=0;k<p;k++)
                x[i][j]+=a[i][k]*b[j][k];
		}
	}
}

/* Dimensions: [A]=pxn, [B]=pxm and [X]=nxm */
void nrAtB(double **a, double **b, int p, int n, int m, double **x)
{
	int i,j,k;
	for (i=0;i<n;i++) {
		for (j=0;j<m;j++)
		{
			x[i][j]=0.0;
			for (k=0;k<p;k++)
                x[i][j]+=a[k][i]*b[k][j];
		}
	}
}

void nrAddMat(double **a, double **b, int n, int m, double **x)
{
	int i,j;
	for (i=0;i<n;i++)
		for (j=0;j<m;j++)
			x[i][j]=a[i][j]+b[i][j];
}

void nrSubMat(double **a, double **b, int n, int m, double **x)
{
	int i,j;
	for (i=0;i<n;i++)
		for (j=0;j<m;j++)
			x[i][j]=a[i][j]-b[i][j];
}


void nrScaMat(double **a, int n, int m, double s, double **x)
{
	int i,j;
	for (i=0;i<n;i++)
		for (j=0;j<m;j++)
			x[i][j]=s*a[i][j];
}



void nrAt(double **a, int n, int m, double **x)
{
	int i,j;
	for (i=0;i<n;i++)
		for (j=0;j<m;j++)
			x[j][i]=a[i][j];
}

void nrMatCopy(double **a, int n, int m, double **x)
{
	int i,j;
	for (i=0;i<n;i++) 
		for (j=0;j<m;j++) 
			x[i][j]=a[i][j];
}

void nrAddVec(double *u, double *v, int n, double *x)
{
	int i;
	for (i=0;i<n;i++) x[i]=u[i]+v[i];
}

void nrSubVec(double *u, double *v, int n, double *x)
{
	int i;
	for (i=0;i<n;i++) x[i]=u[i]-v[i];
}

void nrScaVec(double *u, int n,  double s, double *x)
{
	int i;
	for (i=0;i<n;i++) x[i]=s*u[i];
}

double nrDotProd(double *u, double *v, int n)
{
	double tmp=0;
	int i;
	for (i=0;i<n;i++) 
		tmp+=u[i]*v[i];
	return tmp;
}


/*
 * Utilities functions
 */

double **nrNewMat(int n, int m)
{
	int i;
	double **a = (double **) malloc(n*sizeof(double*));
    for (i=0;i<n;i++) 
	  a[i] = (double *) malloc(m*sizeof(double));
	return a;
}

double**nrNewMatCopy(double **a, int n, int m)
{
	double **x = nrNewMat(n,m);
	int i,j;

	for (i=0;i<n;i++) 
		for (j=0;j<m;j++) 
			x[i][j]=a[i][j];

	return x;
}


double ** nrNewMatInit( int n, int m, ... ) 
{
	va_list arglist;
	double **a=nrNewMat(n,m);
	int i,j;
	double val;

	va_start( arglist, m );
	for(i=0;i<n;i++) 
	{
		for(j=0;j<m;j++)
    {
      val = va_arg( arglist, double );
			a[i][j] = (double)val;
		}
	}
	va_end( arglist );
	return a;
}


double **nrNewIdentMat(int n)
{
	int i,j;
	double **a = (double **) malloc(n*sizeof(double*));
    for (i=0;i<n;i++) 
	  a[i] = (double *) malloc(n*sizeof(double));

    for (i=0;i<n;i++) 
	{
		for (j=0;j<n;j++) a[i][j]=0;
		a[i][i]=1;
	}

	return a;
}


void nrZeroMat(double **a, int n, int m)
{
	int i,j;
    for (i=0;i<n;i++) 
		for (j=0;j<m;j++) 
			a[i][j]=0;
}



void nrFreeMat(double **a, int n)
{
	int i;
    for (i=0;i<n;i++)
          if (a[i]) free(a[i]);
	if (a) free(a);
}


double *nrNewVec(int n)
{
	double * x  = (double *) malloc(n*sizeof(double));
	return x;
}

double *nrNewVecCopy(double *v, int n)
{
	int i;
	double * x  = (double *) malloc(n*sizeof(double));
    for (i=0;i<n;i++) x[i]=v[i]; 
	return x;
}

double *nrNewVecInit(int n, ...)
{
	va_list arglist;
	int i;
	double val;
	double *a=nrNewVec(n);

	va_start( arglist, n );
	for(i=0;i<n;i++) 
	{
    val = va_arg( arglist, double );
   	a[i] = (double)val;
	}
	va_end( arglist );
	return a;
}

void nrZeroVec(double *v, int n)
{
	int i;
    for (i=0;i<n;i++) v[i]=0;
}

void nrFreeVec(double *x)
{
	if (x) free(x);
}




void nrShowMat(char * title, double **a, int n, int m)
{
	int i,j;
	printf("\n%s",title);
	for (i=0;i<n;i++){
		for (j=0;j<m;j++){
            if ((j%5)==0) printf("\n");
			printf("%12.6f ", a[i][j]);
		}
	}
}
void nrShowVec(char * title, double *v, int n)
{
	int j;
	printf("\n%s",title);
	for (j=0;j<n;j++){
            if ((j%5)==0) printf("\n");
			printf("%12.6f ", v[j]);
		}
}

void nrAddMatTensor(double** a, double *v, int n)
{
	int i,j;

   for (i=0;i<n;i++) {
      for (j=0;j<n;j++) { 
			a[i][j]+=(v[i]*v[j]);
      }
   }

}

