load "medit"
//
// SET UP THE PROBLEM AND READ DATA
//
include "settingANDproblemOnline.idp";
//
// FOR GRAPHICAL OUTPUT
//
border OX(t=2,3){x=t;y=0;z=0;}
border OY(t=-1,1){x=2;y=t;z=0;}

//
// LOOP ON PARAMETER SPACE TO CHECK RESULTS
//
int NN=k0max-k0min;
real ds=1;
real[int] errBR(NN+1),errOUTPUT(NN+1),k0vec(NN+1),Tiempo(NN+1);

k0vec=k0min:ds:k0max;
for(int jj=0;jj<=k0max-k0min;jj++)
{

	k0=k0min+jj*ds;

real cpu=clock();

   include "parameters.idp";

//
// COMPUTE GLOBAL MATRIX PARAMETER (DEPENDENT NOW)
//
for(int i=0;i<samples;i++)
{

for(int j=0;j<samples;j++)
         {        
 Ro(i,j)=k0*A0(i,j)+kleft*A1(i,j)+kright*A2(i,j)+Bi*A3(i,j)+aa*A4(i,j);  

         }
}

//
// COEFFICIENTS OF SOLUTION
//
real[int] alphao(samples);
//
// HERE WE KEEP SOLUTION
//
Vh RBufinal=0;

RBufinal=0;
matrix So;
So=Ro;
set(So,solver=sparsesolver);
alphao=So^-1*Rrhsso;

//
// COMPUTATION OF RBu USING alpha ARRAY
//
for(int i=0;i<samples;i++)
{
RBufinal=RBufinal+alphao[i]*usample[i];
}
//
// GET TIME OF RB COMPUTATION
//
real RBt=clock()-cpu;
// GET OUTPUT
real aux1=int2d(Th3,700)(RBufinal);

cout<<" CPU TIME WITH REDUCED BASIS = "<<RBt<<endl;
cout<<" OUTPUT WITH REDUCED BASIS = "<<aux1<<endl;
int mm=40;
real[int] viso(mm+1);
for (int i=0;i<viso.n;i++)
viso[i]=RBufinal[].min+i*(RBufinal[].max-RBufinal[].min)/mm;

//plot(cmm="REDUCED BASIS COMPUTATION FOR k0 ="+k0,RBufinal,nbiso=40,value=1,fill=1
//,viso=viso(0:viso.n-1),hsv=colorhsv);

//
// COMPUTATION FEM SOLUTION
// 
      cpu=clock();
		  include "ParametersAndSolve.idp";
//	medit("Solution ",Th3,u);
   real FEMt=clock()-cpu; // TIME
	 real aux2=int2d(Th3,700)(u); // OUTPUT WITH FEM
	  cout << " CPU TIME NO REDUCED BASIS =  " << FEMt << endl;
	  cout<<" OUTPUT NO REDUCED BASIS = "<<aux2<<endl;

//
// TRANSLATE SOLUTION TO PLOT BOTH AT THE SAME TIME 
//
Vhc uc;
uc=u(x,y-5,z);
//
// TAKE AVERAGE VALUES FOR THE ISOLINES
//
for (int i=0;i<viso.n;i++)
viso[i]=(RBufinal[].min+u[].min)
    +i*((RBufinal[].max+u[].max)-(RBufinal[].min+u[].min))/mm;
viso=0.5*viso;  
//  	
plot(cmm="REDUCED BASIS vs FEM SOLUTION for k0 ="+k0,RBufinal,uc,OX(3),OY(1),value=1,fill=1,viso=viso(0:viso.n-1),hsv=colorhsv);

//
// ABSOLUTE ERROR EN OUTPUT 
//
 errOUTPUT[jj]=abs(aux1- aux2);
 //
// GAIN FACTOR
//
 Tiempo[jj]=FEMt/RBt;
//
// ABSOLUTE ERROR IN SOLUTIONS 
// 
Vh error=RBufinal-u;
//
//  NORMA  DE ERROR 
//
real part1= int3d(Th3)( dx(error)^2+dy(error)^2+dz(error)^2); 
real part2=int3d(Th3)(error^2);                   
real normdiff=sqrt(part1+part2*invvol);  
  errBR[jj]= normdiff;            
cout <<"ABSOLUTE ERROR BETWEEN SOLUTIONS "<<errBR[jj]<<" FOR h= "<<1.0/nn<<endl;
cout <<"ABSOLUTE ERROR BETWEEN OUTPUTS "<<errOUTPUT[jj]<<" FOR h= "<<1.0/nn<<endl;
cout <<"TIME GAIN FACTOR "<<Tiempo[jj]<<" FOR h= "<<1.0/nn<<endl;
cout <<"REDUCED BASIS DIM =  "<<samples<<" FOR h= "<<1.0/nn<<endl;
}


real[int] OX1(2),OX2(2),OY1(2),OY2(2);
OX1[0]=k0vec.min;
OX1[1]=k0vec.max;
OX2[0]=0;
OX2[1]=0;

OY1[0]=0;
OY1[1]=0;
OY2[0]=0;
OY2[1]=errBR.max;

plot([k0vec,errBR],cmm="SOLUTION: ABSOLUTE ERRORS FOR h = "
+1.0/nn+ ". Max "+errBR.max+". Min "+errBR.min,value=1,ps="errRBnn"+nn+".eps");

//plot([OX1,OX2],[OY1,OY2],[k0vec,errBR],cmm="SOLUTION: ABSOLUTE ERRORS FOR h = "
//+1.0/nn+ ". Max "+errBR.max+". Min "+errBR.min,value=1,ps="errRBnn"+nn+".eps");

OY2[1]=errOUTPUT.max;
plot([k0vec,errOUTPUT],cmm="OUTPUT: ABSOLUTE ERRORS FOR h = "
+1.0/nn+ ". Max "+errOUTPUT.max+". Min "+errOUTPUT.min,value=1,ps="errOUTPUTnn"+nn+".eps");
//plot([OX1,OX2],[OY1,OY2],[k0vec,errOUTPUT],cmm="OUTPUT: ABSOLUTE ERRORS FOR h = "
//+1.0/nn+ ". Max "+errOUTPUT.max+". Min "+errOUTPUT.min,value=1,ps="errOUTPUTnn"+nn+".eps");

plot([k0vec,Tiempo],cmm="TIME GAIN FACTOR :  FOR h = "
+1.0/nn+ ". Max "+Tiempo.max+". Min "+Tiempo.min,value=1,ps="Tiemponn"+nn+".eps");
//plot([OX1,OX2],[OY1,OY2],[k0vec,Tiempo],cmm="TIME GAIN FACTOR :  FOR h = "
//+1.0/nn+ ". Max "+Tiempo.max+". Min "+Tiempo.min,value=1,ps="Tiemponn"+nn+".eps");