SuperLU.edp

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 


problem laplace(uh,vh,solver=sparsesolver,tgv=1e5) =                    //  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,ps="Laplace.eps",value=true); 
} 

// 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=sparsesolver,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)defaulttoGMRES(); 
  if(i==1)defaultsolver(); 
} 
cout << " fin.. \n";