/*!\file: Reducematrixfromgton
 * \brief
 */ 

#ifdef HAVE_CONFIG_H
	#include "config.h"
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#undef __FUNCT__ 
#define __FUNCT__ "Reducematrixfromgton"

#include "../toolkits/toolkits.h"


void  Reducematrixfromgton(Mat* pKnn,Mat Kgg,Mat Gmn,double* pv_m,int msize, double* pv_n,int nsize,int flag){


	/*output: */
	Mat Knn=NULL;

	/*intermediary: */
	Mat Kmm=NULL;
	Mat Kmn=NULL;
	Mat Knm=NULL;
	Mat KnmGmn=NULL;
	Mat KmmGmn=NULL;
	Mat tGmn=NULL;
	Mat tGmnKmn=NULL;
	PetscScalar a=1;

	if (msize){
		
		/*Partition K_gg*/
		MatPartition(&Kmm, Kgg, pv_m,msize,pv_m,msize); 
		MatPartition(&Kmn, Kgg, pv_m,msize,pv_n,nsize); 
		MatPartition(&Knm, Kgg, pv_n,nsize,pv_m,msize); 
		MatPartition(&Knn, Kgg, pv_n,nsize,pv_n,nsize);
	
		/*Reduce K_gg to K_nn*/
		//K_nn = K_nn + K_nm * G_mn;
		MatMatMult(Knm,Gmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&KnmGmn);
		MatAXPY(Knn,a,KnmGmn,DIFFERENT_NONZERO_PATTERN);

		if (flag!=1){
			//K_mn = K_mn + K_mm * G_mn;
			MatMatMult(Kmm,Gmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&KmmGmn);
			MatAXPY(Kmn,a,KmmGmn,DIFFERENT_NONZERO_PATTERN);
			
			//K_nn = K_nn + G_mn' * K_mn;
			#if _PETSC_VERSION_ == 2
			MatTranspose(Gmn,&tGmn);
			#else
			MatTranspose(Gmn,MAT_INITIAL_MATRIX,&tGmn);
			#endif
			MatMatMult(tGmn,Kmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&tGmnKmn);
			MatAXPY(Knn,a,tGmnKmn,DIFFERENT_NONZERO_PATTERN);
		}
	}
	else{
		/*msize=0 just return input*/
		MatDuplicate(Kgg,MAT_COPY_VALUES,&Knn);
	}	

	/*Free ressources:*/
	MatFree(&Kmm);
	MatFree(&Kmn);
	MatFree(&Knm);
	MatFree(&KnmGmn);
	MatFree(&KmmGmn);
	MatFree(&tGmn);
	MatFree(&tGmnKmn);
	
	/*Assign output pointers:*/
	*pKnn=Knn;
}