/*
 * GradJCompute.c:
 */

#include "../../../config.h"

#if defined(_PARALLEL_) && defined(_HAVE_PETSC_)

#include "../include/cielo.h"
#include "../modules.h"
#include "./parallel.h"

#undef __FUNCT__ 
#define __FUNCT__ "GradJCompute"
#undef CLEANUP
#define CLEANUP GradJComputeLocalCleanup();

void GradJComputeLocalCleanup(void);

int GradJCompute(ParameterInputs* inputs,FemModel* femmodel){
	
	/*Error management: */
	int noerr=1;
	int i;
	int dummy;
	
	Vec* u_g=NULL;
	Vec* du_g=NULL;
	Vec* du_f=NULL;
	Vec* lambda_f=NULL;
	Vec* lambda_g=NULL;
	double* lambda_g_double=NULL;
	double* u_g_double=NULL;
	Vec* gradj_g=NULL;
	Mat* K_ff0=NULL;
	Mat* K_fs0=NULL;
	
	//Recover uset: */
	uset=femmodel->uset; //external variable

	//Recover solution for this stiffness and right hand side: 
	cielodiagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,inputs,femmodel);

	//Buid Du, difference between observed velocity and model velocity.
	VecToMPISerial(&u_g_double,u_g); VecFree(&u_g);
	
	Dux(&du_g,femmodel->bgpdt, femmodel->bgpdtb, femmodel->est, femmodel->ept, femmodel->mpt,femmodel->lst, u_g_double,femmodel->workspaceparams->u_g_obs,inputs,femmodel->batchparams->analysis_type_enum);


	//Reduce adjoint load from g-set to f-set
	Reducerightside(&du_f, du_g, femmodel->G_mn, K_fs0, femmodel->y_s0, femmodel->flag_y_s0, uset->pv_m, uset->msize, uset->pv_n, uset->nsize, uset->pv_f, uset->fsize, uset->msize,uset->ssize,uset->fsize);
	VecFree(&du_g);MatFree(&K_fs0);

	
	//Solve for adjoint vector: 
	lambda_f=xmalloc(sizeof(Vec));
	Solverx(lambda_f,&dummy,K_ff0,dummy,dummy,du_f,dummy,NULL,0,femmodel->batchparams->solverstring);
	VecFree(&du_f);
	MatFree(&K_ff0);

	
	//Merge back to g set
	Mergesolvec( &lambda_g, lambda_f, femmodel->G_mn, femmodel->y_s0, uset->pv_m, uset->msize, uset->pv_n, uset->nsize, uset->pv_f, uset->fsize, uset->pv_s, uset->ssize, uset->gsize,uset->msize,uset->nsize);
	VecFree(&lambda_f);

	//Compute gradJ for each parameter
	VecToMPISerial(&lambda_g_double,lambda_g);VecFree(&lambda_g);
	
	for(i=0;i<femmodel->workspaceparams->num_control_parameters;i++){
		char* control_type=femmodel->workspaceparams->control_types[i];
	
		Gradjx( &gradj_g, femmodel->bgpdt, femmodel->bgpdtb, femmodel->est, femmodel->ept, femmodel->mpt,femmodel->lst, u_g_double,lambda_g_double,inputs,femmodel->batchparams->analysis_type_enum,control_type);

		WorkspaceParamsSetParameterGradient(femmodel->workspaceparams,gradj_g,control_type);
	}
	
	/*free ressources: */
	xfree((void**)&lambda_g_double);
	xfree((void**)&u_g_double);


	EXIT(noerr);
}

void GradJComputeLocalCleanup(void){
	return;
}
#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_)

