/*!\file NormalizeConstraintsx
 * \brief normalize and reduce constraints matrix Rmg to Gmn:
 */

#include "./NormalizeConstraintsx.h"

#include "../shared/shared.h"

#undef __FUNCT__ 
#define __FUNCT__ "NormalizeConstraintsx"

void NormalizeConstraintsx( Mat* pGmn, Mat Rmg, NodeSets* nodesets){

	/*output: */
	Mat Gmn=NULL;

	/*intermediary: */
	int i;
	Mat Rmm=NULL;
	Mat Rmn=NULL;
	Mat InvRmm=NULL;
	int msize,nsize;
	double* pv_m=NULL;
	double* row_m=NULL;
	double* pv_n=NULL;
	int     null=0;

	/*Recover data: */
	msize=nodesets->GetMSize();
	nsize=nodesets->GetNSize();
	pv_m=nodesets->GetPV_M();
	pv_n=nodesets->GetPV_N();


	if (msize){

		/*Build row partitioning vector for Rmm and Rmn: */
		row_m=(double*)xmalloc(msize*sizeof(double));
		for(i=0;i<msize;i++)row_m[i]=i+1; //matlab indexing

		/*Partition Rmg into Rmm and Rmn: */
		MatPartition(&Rmm,Rmg,row_m,msize,pv_m,msize); //equivalent of Rmm=Rmg(:,mset);
		MatPartition(&Rmn,Rmg,row_m,msize,pv_n,nsize);

		/*Invert Rmm: */
		MatInvert(&InvRmm,Rmm);

		/*Do Gmn=-(Rmm-1)*Rmn :*/
		MatScale(InvRmm,-1.0);
		MatMatMult(InvRmm,Rmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Gmn);

	}
	else{
		Gmn=NULL;
	}

	/*Free ressources:*/
	xfree((void**)&row_m);
	MatFree(&Rmm);
	MatFree(&Rmn);
	MatFree(&InvRmm);
	
	/*Assign output pointers:*/
	*pGmn=Gmn;
}
