/*!\file SeqMat.cpp
 * \brief: implementation of the SeqMat object
 */

/*Headers:*/
/*{{{*/
#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <stdio.h>
#include <string.h>
#include "../toolkits.h"
#include "../../shared/shared.h"

/*}}}*/

/*SeqMat constructors and destructor*/
/*FUNCTION SeqMat::SeqMat(){{{*/
SeqMat::SeqMat(){

	this->M=0;
	this->N=0;
	this->matrix=NULL;
}
/*}}}*/
/*FUNCTION SeqMat::SeqMat(int M,int N){{{*/
SeqMat::SeqMat(int pM,int pN){

	this->M=pM;
	this->N=pN;
	this->matrix=NULL;
	if(M*N) this->matrix=xNewZeroInit<IssmDouble>(pM*pN);
}
/*}}}*/
/*FUNCTION SeqMat::SeqMat(int M,int N, IssmDouble sparsity){{{*/
SeqMat::SeqMat(int pM,int pN, IssmDouble sparsity){

	this->M=pM;
	this->N=pN;
	this->matrix=NULL;
	if(M*N) this->matrix=xNewZeroInit<IssmDouble>(pM*pN);
}
/*}}}*/
/*FUNCTION SeqMat::SeqMat(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity){{{*/
SeqMat::SeqMat(IssmDouble* serial_mat,int pM,int pN,IssmDouble sparsity){

	int i,j;

	this->M=pM;
	this->N=pN;
	this->matrix=NULL;
	if(M*N){
		this->matrix=xNewZeroInit<IssmDouble>(pM*pN);
		xMemCpy<IssmDouble>(this->matrix,serial_mat,pM*pN);
	}

}
/*}}}*/
/*FUNCTION SeqMat::SeqMat(int M,int N, int connectivity, int numberofdofspernode){{{*/
SeqMat::SeqMat(int pM,int pN, int connectivity,int numberofdofspernode){

	this->M=pM;
	this->N=pN;
	this->matrix=NULL;
	if(M*N) this->matrix=xNewZeroInit<IssmDouble>(pM*pN);
}
/*}}}*/
/*FUNCTION SeqMat::~SeqMat(){{{*/
SeqMat::~SeqMat(){

	xDelete<IssmDouble>(this->matrix);
	M=0;
	N=0;
}
/*}}}*/

/*SeqMat specific routines: */
/*FUNCTION SeqMat::Echo{{{*/
void SeqMat::Echo(void){

	int i,j;
	_printLine_("SeqMat size " << this->M << "-" << this->N);
	for(i=0;i<M;i++){
		for(j=0;j<N;j++){
			_printString_(this->matrix[N*i+j] << " ");
		}
		_printLine_("");
	}
}
/*}}}*/
/*FUNCTION SeqMat::Assemble{{{*/
void SeqMat::Assemble(void){
		
	/*do nothing*/

}
/*}}}*/
/*FUNCTION SeqMat::Norm{{{*/
IssmDouble SeqMat::Norm(NormMode mode){

	IssmDouble norm;
	IssmDouble absolute;
	int i,j;

	switch(mode){
		case NORM_INF:
			norm=0;
			for(i=0;i<this->M;i++){
				absolute=0;
				for(j=0;j<this->N;j++){
					absolute+=fabs(this->matrix[N*i+j]);
				}
				norm=max(norm,absolute);
			}
			return norm;
			break;
		default:
			_error2_("unknown norm !");
			break;
	}
}
/*}}}*/
/*FUNCTION SeqMat::GetSize{{{*/
void SeqMat::GetSize(int* pM,int* pN){

	*pM=this->M;
	*pN=this->N;

}
/*}}}*/
/*FUNCTION SeqMat::GetLocalSize{{{*/
void SeqMat::GetLocalSize(int* pM,int* pN){
	
	*pM=this->M;
	*pN=this->N;

}
/*}}}*/
/*FUNCTION SeqMat::MatMult{{{*/
void SeqMat::MatMult(SeqVec* X,SeqVec* AX){

	int i,j;
	int XM,AXM;
	IssmDouble dummy;

	X->GetSize(&XM);
	AX->GetSize(&AXM);

	if(M!=AXM)_error2_("A and AX should have the same number of rows!");
	if(N!=XM)_error2_("A and X should have the same number of columns!");

	for(i=0;i<M;i++){
		dummy=0;
		for(j=0;j<N;j++){
			dummy+= this->matrix[N*i+j]*X->vector[j];
		}
		AX->vector[i]=dummy;
	}

}
/*}}}*/
/*FUNCTION SeqMat::Duplicate{{{*/
SeqMat* SeqMat::Duplicate(void){

	IssmDouble dummy=0;

	return new SeqMat(this->matrix,this->M,this->N,dummy);

}
/*}}}*/
/*FUNCTION SeqMat::ToSerial{{{*/
IssmDouble* SeqMat::ToSerial(void){

	IssmDouble* buffer=NULL;

	if(this->M*this->N){
		buffer=xNew<IssmDouble>(this->M*this->N);
		xMemCpy<IssmDouble>(buffer,this->matrix,this->M*this->N);
	}
	return buffer;

}
/*}}}*/
/*FUNCTION SeqMat::SetValues{{{*/
void SeqMat::SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){
	
	int i,j;
	switch(mode){
		case ADD_VAL:
			for(i=0;i<m;i++) for(j=0;j<n;j++) this->matrix[N*idxm[i]+idxn[j]]+=values[n*i+j];
			break;
		case INS_VAL:
			for(i=0;i<m;i++) for(j=0;j<n;j++) this->matrix[N*idxm[i]+idxn[j]]=values[n*i+j];
			break;
		default:
			_error2_("unknown insert mode!");
			break;
	}

}
/*}}}*/
/*FUNCTION SeqMat::Convert{{{*/
void SeqMat::Convert(MatrixType type){

	/*do nothing*/

}
/*}}}*/
