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

/*
 *	Finds a non-trivial solution x to the system Ax=0
 *  A mxn, with m>n and rank(a)< n (preferably = n-1).
 *
 *  Return value indicates error. It's the quotient of the
 *  smallest by the biggest singular value of A.  If rank(a)
 *  is less than n, value should be zero.  (~=1.0e-6 to this
 *  implementation with simple precision).
 */
double nrSVDAx0(double **a, int m, int n, double *x);


/*
 *	Finds a solution to the system Ax=b, A mxn with m>n
 */
void  nrSVDAxb(double **a, int m, int n, double *b, double *x);

/*
 * Computes matrices [U],[W] e [V] of a single value decomposition
 * [A] = [U][W][V]T
 * [A] mxn, [U]mxn (orthogonal columns), [W]nxn (diagonal) e [V] nxn (orthogonal).
 *
 */ 
void nrSVD(double **a, int m, int n, double **u, double *w, double **v);

/*
 * Finds solution x to the system Ax=b
 * A nxn, throuhg LU decomposition.
 *
 * Adapted from Numerical Recipes. 
 */ 
void nrLUAxb(double **a, int n, double *b, double *x);


/*
 * 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);


/* Convenience Functions*/

/*
 * Computes the product of matrix [A] by the vector {b} 
 * {x} is the result.
 * Dimensions: [A]=nxm, {b}=m and {x}=n.
 * @warning  {x} CANNOT use the same memory area as {b}
 */
void nrAb(double **a, double *b, int n, int m, double *x);

/*
 * Computes the product of the transpose of matrix [A] 
 * by the vector {b}. {x} is the result.
 * Dimensions: [A]=mxn, {b}=m and {x}=n.
 * @warning  {x} CANNOT use the same memory area as {b}
 */
void nrAtb(double **a, double *b, int n, int m, double *x);


/*
 * Computes the product of matrix [A] by the matrix [B] 
 * [X] is the result.
 * Dimensions: [A]=nxp, [B]=pxm and [X]=nxm.
 * @warning  [X] CANNOT use the same memory area as [B]
 */
void nrAB(double **a, double **b, int n, int p, int m, double **x);


/*
 * Computes the product of matrix [A] by the transpose of 
 * matrix [B]. [X] is the result.
 * Dimensions: [A]=nxp, [B]=mxp and [X]=nxm.
 * @warning  [X] CANNOT use the same memory area as [B]
 */
void nrABt(double **a, double **b, int n, int p, int m, double **x);

/*
 * Computes the product of the transpose of matrix [A] by the  
 * matrix [B]. [X] is the result.
 * Dimensions: [A]=pxn, [B]=pxm and [X]=nxm.
 * @warning  [X] CANNOT use the same memory area as [B]
 */
void nrAtB(double **a, double **b, int p, int n, int m, double **x);

/*
 * Computes the transpose of matrix [A]. 
 * [X] is the result.
 * Dimensions: [A]=nxm and  [X]=mxn.
 * @warning  [X] CANNOT use the same memory area as [A]
 */
void nrAt(double **a, int n, int m, double **x);


/*
 * Computes [X] = [A] + [B]
 * Dimensions: [A],[B]  and  [X]=nxm.
 */
void nrAddMat(double **a, double **b, int n, int m, double **x);

/*
 * Computes [X] = [A] - [B]
 * Dimensions: [A],[B]  and  [X]=nxm.
 */
void nrSubMat(double **a, double **b, int n, int m, double **x);

/*
 * Computes [X] = [A]
 * Dimensions: [A]  and  [X]=nxm.
 */
void nrMatCopy(double **a, int n, int m, double **x);

/*
 * Computes [X] = s [A] (s is a scalar value).
 * Dimensions: [A] and  [X]=nxm. 
 */
void nrScaMat(double **a, int n, int m, double s, double **x);

/*
 * Computes {x} = {u} + {v}
 * Dimensions: {x},{u} and {v} = n.
 */
void nrAddVec(double *u, double *v, int n, double *x);

/*
 * Computes {x} = {u} - {v}
 * Dimensions: {x},{u} and {v} = n.
 */
void nrSubVec(double *u, double *v, int n, double *x);

/*
 * Computes {x} = x {v} (s is a scalar value).
 * Dimensions: {x},{u}= n.
 */
void nrScaVec(double *u, int n,  double s, double *x);



/*
 * Computes the dot product of two vectors of dimension n.
 */
double nrDotProd(double *u, double *v, int n);


/** 
 * Memory allocation utilities
 * Note that matrices are allocated in memory using 
 * common C techniques for dinamically allocated matrices.
 * That is:
 *   1) the (i,j) element of a matrix "a" can be addressed as a[i][j]
 *   2) indices of a nxm matrix range from [0,n-1] and [0,m-1].
 *   3) matrices are allocades in a row by row scheme.
 */


/*
 *  Creates and returns a standard C dynamically allocated matrix
 *  of dimensions nxm (n rows and m columns).
 *  
 */
double **nrNewMat(int n, int m);

/*
 * Create and returns a new standard C dynamically allocated matrix
 * that is a copy matrix [A]. [A] nxm. 
 * @warning  You can change or free [A] without affecting its copy.
 */
double **nrNewMatCopy(double **a, int n, int m);

/*
 * Create and returns a new standard C dynamically allocated matrix [I].
 * [I]nxn is the identity matrix.
 */
double **nrNewIdentMat(int n);

/**
 *   Creates and returns a new standard C dynamically allocated matrix
 *   with the elements given in the parameter list.
 *
 *   @parm n number of rows
 *   @parm m number of columns
 *   @parm ... matrix elements given by row.
 *
 *   Exemple: to obtain the matrix
 *           1 -2 5
 *           2  4 3
 *   use:
 *          double **a=nrCreateMAtrixFromData(2,3, 1.,-2.,5.,  2.,4.,3.);
 * 
 *   @warning allways express values as reals, i.e., use "1.0" instead of "1".
 */
double ** nrNewMatInit( int n, int m, ... );


/*
 * Clears (assings zero) all components of the matrix [A].
 * Dimensions of [A]=nxm.
 */
void nrZeroMat(double **a, int n, int m);


void nrFreeMat(double **a, int n);


/*
 *  Returns a standard C dynamic allocated array.
 */
double *nrNewVec(int n);

/*
 * Returns a new standard C dynamically allocated vector
 * that is a copy of vector {v}. Dimension of {v}=n.
 * @warning  You can change or free {v} without affecting its copy.
 */
double *nrNewVecCopy(double *v, int n);

/*
 * Clears (assings zero) all components of the vector {v}.
 * Dimension of {v}=n.
 */
void nrZeroVec(double *v, int n);


/*
 * Returns a new standard C dynamically allocated vector
 * that is initialized with the elements given in the parameter list.
 * Dimension of {v}=n.
 */
double *nrNewVecInit(int n, ...);


void nrFreeVec(double *x);


/***********************
 * Convenient functions*
 ***********************/

/*
 * Prints a matrix on the stdio.
 */
void nrShowMat(char * title, double **a, int n, int m);

/*
 * Prints a vector on the stdio.
 */
void nrShowVec(char * title, double *v, int n);

/*
 * Add to matrix a[n][n] the the tensor product of a vector
 * a[i][j]+=v[i]*v[j]
 */
void nrAddMatTensor(double** a, double *v, int n);



