load "SuperLu"

verbosity=0;

{

cout << "laplace solving with SuperLu" << endl;

mesh Th=square(10,10);
fespace Vh(Th,P1);     // P1 FE space
 Vh uh,vh;              // unkown and test function. 
 func f=1;                 //  right hand side function 
 func g=0;                 //  boundary condition function
 
 verbosity= 3;
problem laplace(uh,vh,solver="SUPERLU",tgv=1e5,sym=0) =                    //  definion of  the problem 
    int2d(Th)( dx(uh)*dx(vh) + dy(uh)*dy(vh) ) //  bilinear form
  - int2d(Th)( f*vh )                        	     //  linear form
  + on(1,2,3,4,uh=g) ;                      //  bou  ndary condition form

  laplace; // solve the problem plot(uh); // to see the result
 
  plot(uh,value=true); 
   verbosity= 1;
   assert(abs(uh[].max - 0.0730987)< 0.01);
}

// FFCS: need to declare it globally to print out its value for regression tests
complex[int] lastx(4);

for(int i=0;i<3;++i)
{
  if(i==0)  cout << "resolution SuperLU" <<endl; 
  if(i==1)  cout << "resolution GMRES" <<endl;
  if(i==2)  cout << "resolution UMFPACK" <<endl; 
  {
    matrix A = 
      [[ 0,  1,  0, 10],
       [ 0,  0,  2,  0],
       [ 0,  0,  0,  3],
       [ 4,  0,  0,  0]];
    real[int] xx = [ 4,1,2,3], x(4), b(4); // xb(4),bbb(4);
    b=A*xx;
    cout << "b="  << b  << endl;
    cout << "xx=" << xx << endl;
    set(A,solver="SUPERLU",sparams="DiagPivotThresh=0.05,ColPerm=MMD_AT_PLUS_A,Equil=NO"); 
    x = A^-1*b;	
    cout << "x=" << endl; cout << x << endl;	
  }
  
  {
    matrix<complex> A = 
      [[  0, 1i,  0, 10],
       [  0,  0, 2i,  0],
       [  0,  0,  0, 3i],
       [ 4i,  0,  0,  0]];
    complex[int] xx = [ 4i,1i,2i,3i], x(4), b(4);
    b = A*xx;
    cout << "b="  << b << endl;
    cout << "xx=" << xx << endl;
    set(A,solver=sparsesolver);
    x = A^-1*b;
    cout << "x=" << endl; cout << x << endl;
    lastx=x;
  }
  if(i==0) DefaultSolver="GMRES";
  if(i==1) DefaultSolver="UMFPACK";
}
cout << " fin.. \n";




