/*!\file Reduceloadfromgtofx
 * \brief reduce loads from g set to f set 
 */

#include "./Reduceloadfromgtofx.h"

#undef __FUNCT__ 
#define __FUNCT__ "Reduceloadfromgtofx"

void	Reduceloadfromgtofx( Vec* ppf, Vec pg, Mat G_mn, Mat Kfs, Vec y_s, NodeSets* nodesets){

	/*output: */
	Vec pf=NULL;

	/*intermediary*/
	Vec pn=NULL;
	Vec pm=NULL;
	Mat tG_mn=NULL;
	int tG_mnm,tG_mnn;
	Vec tG_mnpm=NULL;
	Vec Kfsy_s=NULL;
	int Kfsm,Kfsn;
	PetscScalar a;
	

	if(!pg){
		pf=NULL;
	}
	else{

		/* Reduce pg to pn:*/
		if (nodesets->GetMSize()){

			VecPartition(&pn, pg, nodesets->GetPV_N(), nodesets->GetNSize());
			VecPartition(&pm, pg, nodesets->GetPV_M(), nodesets->GetMSize());
		
			/*Create G_mn'*/
			MatTranspose(G_mn,&tG_mn);  

			/*Create G_mn' * pm*/
			MatGetLocalSize(tG_mn,&tG_mnm,&tG_mnn);
			tG_mnpm=NewVecFromLocalSize(tG_mnm);
			MatMultPatch(tG_mn,pm,tG_mnpm);

			/*add pn to G_mn' * pm*/ 
			a=1;
			VecAXPY(pn,a,tG_mnpm);  
		}
		else{
			VecDuplicate(pg,&pn);  
			VecCopy(pg,pn);  
		}
		
		/*% Reduce pn to pf:*/

		if (nodesets->GetSSize()){

			VecPartition(&pf, pn, nodesets->GetPV_F(),nodesets->GetFSize());

			/*pf = pf - Kfs * y_s;*/
			MatGetLocalSize(Kfs,&Kfsm,&Kfsn);
			Kfsy_s=NewVecFromLocalSize(Kfsm);
			MatMultPatch(Kfs,y_s,Kfsy_s);

			a=-1;
			VecAXPY(pf,a,Kfsy_s);  
		}
		else{
			/*pf=pn*/
			VecDuplicate(pn,&pf);  
			VecCopy(pn,pf);  
		}
	}
	
	
	/*Assign correct pointer*/
	*ppf=pf;

	/*Free ressources and return*/
	VecFree(&pn);
	VecFree(&pm);
	MatFree(&tG_mn);
	VecFree(&tG_mnpm);
	VecFree(&Kfsy_s);

}
