/*This is a collection of functions about the basic matrix operations
 * for solving linear equations and related problems. In this package,
 * matrices are assumed to be N-by-N square matrix, and the vectors
 * to be length of N.
 * Most of the functions in this package are independent of others.
 * 
 * This program was written by Xiaoke Yang @ School of Automation Science and
 * Electrical Engineering, Beihang University. You can copy, modify or 
 * redistribute it freely with this header kept intact. Welcome reports of 
 * bugs, suggestions, etc. yxkmlstar@gmail.com
 *
 * Last modified by Xiaoke Yang, Oct.31,2007
 */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N		3			/*The matrix dimension or vector length*/
#define eps		1e-10		/*The error tolerance or accuracy of computation*/
#define DEBUG	1			/*The flag of debugging mode, 1 for on, 0 for off*/
#define IMAX	200			/*The maximum iteration count, used to terminate */
							/*possible infinite iteration*/

/*Generate a N-by-N random matrix*/
/*not implemented yet*/
void randM(double M[N][N]){
}

/*Extract the diagonal entries from a N-by-N matrix*/
/*not implemented yet*/
void diag(double M[N][N],double D[N]){

}

/*Generate a N-by-N diagonal matrix using the given vector*/
/*not implemented yet*/
void diagg(double v[N],double M[N][N]){
}

/*Generate a N-by-N diagonal matrix of all ones*/
/*not implemented yet*/
void ones(){
}

/*Generate an N-dimensional identity matrix In*/
/*not implemented yet*/
void eye(){
}

/*Get the upper triangular matrix of a square matrix*/
/*not implemented yet*/
void triu(){
}

/*Get the lower triangular matrix of a square matrix*/
/*not implemented yet*/
void tril(){
}

/* Description: copy the values of N-by-N matrix src to Matrix dest
 * Calling syntax: matcopy(src,dest);
 * */
void matcopy(double src[N][N], double dest[N][N]){
	int i,j;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			dest[i][j]=src[i][j];
}

/* Description: copy the values of vector src[N] to vector dest[N]
 * Calling syntax: veccopy(src,dest)
 * */
void veccopy(double src[N],double dest[N]){
	int i;
	for(i=0;i<N;i++)
		dest[i]=src[i];
}

/* Description: print the values of a N-by-N matrix on the screen
 * Calling syntax: matdisp(mat);
 * */
void matdisp(double mat[N][N]){
	int i,j;
	for(i=0;i<N;i++){
		for(j=0;j<N;j++){
			printf("%.7e\t",mat[i][j]);
		}
		printf("\b\n");
	}
}

/* Description: print the values of a vector of length N on the screen
 * Calling syntax: vecdisp(vec);
 * */
void vecdisp(double vec[N]){
	int i;
	for(i=0;i<N;i++){
		printf("%.7e\t",vec[i]);
	}
	printf("\b\n");
}

/* Description: check if a matrix is row diagonally dominant
 * Calling syntax: val=diagdomr(A);
 * */
int diagdomr(double A[N][N]){
	int isdom=1;	/*the return value, 0 for non-dominant, 1 for dominant*/
	int i,j;		/*the iterator*/
	double s;
	for(i=0;i<N;i++){
		s=0;						/*initialize s as 0 for each row*/
		for(j=0;j<N;j++){
			if(j!=i)
				s+=fabs(A[i][j]);	/*the sum of the magnitude of all*/
		}							/*non-diagnoal entries in each row*/
		/*print the debugging information*/
		if(DEBUG==1){
			printf("the sum of magnitude of all non-diagonal entries in row %d is %f.\n",i+1,s);
			printf("the magnitude of diagonal entries in row %d is %f.\n",i+1,fabs(A[i][i]));
		}
		if((fabs(A[i][i])-s)<eps){	/*A(i,i)<s or A(i,i)=s*/
			isdom=0;
			break;
		}
	}
	return isdom;
}

/* Description: check if a matrix is column diagonally dominant
 * Calling syntax: val=diagdomc(A);
 * */
int diagdomc(double A[N][N]){
	int isdom=1;	/*the return value, 0 for non-dominant, 1 for dominant*/
	int i,j;		/*the iterator*/
	double s;
	for(j=0;j<N;j++){
		s=0;						/*initialize s as 0 for each row*/
		for(i=0;i<N;i++){
			if(i!=j)
				s+=fabs(A[i][j]);	/*the sum of the magnitude of all*/
		}							/*non-diagnoal entries in each row*/
		if(fabs(A[j][j])-s<eps){	/*A(j,j)<s or A(j,j)=s*/
			isdom=0;
			break;
		}
	}
	return isdom;
}

/* Description: matrix multiplication, res=A*B
 * Calling syntax: matmult(A,B,res);
 * */
void matmult(double A[N][N],double B[N][N],double res[N][N]){
	int i,j,k;
	double s;
	for(i=0;i<N;i++){			/*i iterates through the row of res[][]*/
		for(j=0;j<N;j++){		/*j iterates through the column of res[][]*/
			s=0;
			for(k=0;k<N;k++)	/*k iterates through row of A and column of B*/
				s+=A[i][k]*B[k][j];
			res[i][j]=s;
		}
	}
}

/* Function: void matmultvec(double A[N][N],double b[N],double res[N])
 * Description: The multiplication of a matrix and a vector, res=A*b
 * Input: A[N][N]-a square matrix; 
 *		  b[N]-a vector of N; 
 *		  res[N]-the resultant vector;
 * Output: none
 * Calling Syntax: matmultvec(A,b,res);
 */
void matmultvec(double A[N][N],double b[N],double res[N]){
	int i,k;
	double s;
	for(i=0;i<N;i++){			/*i iterates through of row of res[]*/
		s=0;
		for(k=0;k<N;k++)		/*k iterates through row of A and column of b*/
			s+=A[i][k]*b[k];
		res[i]=s;
	}
}

/* Description: the qutient of a vector b[N] divided by a constant c, q=b/c.
 * Calling Syntax: vecmultcon(b,c,q);
 * */
void vecmultcon(double b[N], double c, double q[N]){
	int i;
	for(i=0;i<N;i++)
		q[i]=b[i]*c;
}		


/* Description: read data from a file*/
/*not implemented yet*/
void matread(char *fname, double A){
}

/* Description: sum norm or l1 norm of vectors on Rn
 * Calling syntax: val=vecnorm1(vec);
 * */
double vecnorm1(double vec[N]){
	int i;
	double s=0;
	for(i=0;i<N;i++)
		s+=fabs(vec[i]);
	return s;
}

/* Description: Euclidean norm or l2 norm of vectors on Rn
 * Calling syntax: val=vecnorm2(vec);
 * */
double vecnorm2(double vec[N]){
	int i;
	double s=0;
	for(i=0;i<N;i++)
		s+=pow(vec[i],2);
	return sqrt(s);
}

/* Description: max norm or l_inf norm of vectors on Rn
 * Calling syntax: val=vecnorminf(vec);
 * */
double vecnorminf(double vec[N]){
	int i;
	double m;
	m=vec[0];
	for(i=1;i<N;i++)
		if(vec[i]>m)
			m=vec[i];
}

/* Description: lp norm of vectors on Rn, p>=1.
 * Calling syntax: val=vecnormp(vec,p)
 * */
double vecnormp(double vec[N], int p){
	int i;
	double s=0;
	for(i=0;i<N;i++)
		s+=pow(vec[i],(double)p);
	return pow(s,1.0/p);
}

/* Description: minus operation on vectors on Rn, res[N]=x[N]-y[N] subtracts 
 * vector x from y
 * Calling syntax: vecminus(x,y,res);
 * */
void vecminus(double x[N],double y[N],double res[N]){
	int i;
	for(i=0;i<N;i++)
		res[i]=x[i]-y[i];
}

/* Description: plus operation on vectors on Rn, res[N]=x[N]+y[N] adds
 * vectors x and y
 * Calling syntax: vecplus(x,y,res)
 * */
void vecplus(double x[N],double y[N],double res[N]){
	int i;
	for(i=0;i<N;i++)
		res[i]=x[i]+y[i];
}

/* Description: Euclidean distance between two vectors x[N] and y[N] on Rn,
 * this function depends on function vecminus() and vecnorm2().
 * Calling syntax: val=vecdist(x,y);
 * */
double vecdist(double x[N],double y[N]){
	double d[N];
	vecminus(x,y,d);
	return vecnorm2(d);
}

/* Description: inner product or element-by-element multiplication of 
 * two vectors a[N] and b[N] on Rn
 * Calling syntax: val=innerprod(a,b);
 * */
double innerprod(double a[N], double b[N]){
	int i;
	double s;
	s=0;
	for(i=0;i<N;i++)
		s+=a[i]*b[i];
	return s;
}