sparse-cmatrix.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

//  sparse matrix test  ---
// example of the new matrix feature in version 1.40
// -------------------------------------------------
mesh  TH = square(3,4);
mesh  th = square(2,3);
mesh  Th = square(4,4);

complex ccc;
ccc= 1;
cout << ccc << endl;
fespace VH(TH,P1);
fespace Vh(th,P1);
fespace Wh(Th,P1);

matrix RB= interpolate(VH,Vh);  // build interpolation matrix Vh->Vh 
matrix RBB= interpolate(Wh,Vh);  // build interpolation matrix
matrix<complex> B=RB;
B = B*(1+2i);

matrix<complex> BB=RBB;

varf vA(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))+ int1d(Th)(u*v); 
matrix<complex> A=vA(Wh,Wh);
Vh<complex> ml=0;
cout << " ml " << ml[] << endl;
varf vML(u,v) = int2d(th)(1.*v);
ml[]=vML(0,Vh); // build the P1 mass lump of P1
cout << ml[] << endl;
matrix<complex> ML(ml[]); // matrix diagonal
cout << "ML="<<ML << endl;
cout << "B="<<B << endl;
matrix<complex> BML=B*ML; // a faire 
matrix<complex> tB=B';        // transpose and conjugate 
cout << "tB=" << tB << endl;
matrix<complex> MLtB=ML'*B'; // 

cout << "BML="<<BML << endl;
cout << "MLtB=" << MLtB << endl;

// WARNING if UMFPACK is not install
// the UMFPACK solver is replace by LU 
//  but LU need skyline matrix 
if(HaveUMFPACK)
  set(A,solver=UMFPACK); // set a solver 
else 
  set(A,solver=GMRES); // set a solver 
  VH<complex> uH=0;
  Vh<complex> uh=x+y+1i*(x-y);
  uH[]= B*uh[];
  Vh uHr = imag(uH);
  plot(uHr,wait=1);
  matrix<complex> BtA = BB'*A;
  matrix<complex> BtAB = BtA*BB;
if(HaveUMFPACK)  
  set(BtAB,solver=UMFPACK);  
else 
  set(BtAB,solver=GMRES);  
  Vh<complex> ff=1+1i;
  Vh<complex> xx;
  Vh xxr;
  cout << " ------ " << endl;
  
  xx[]=BtAB^-1*ff[];
  cout << " ------ " << endl;
  xx[]=BtAB^-1*ff[];
  cout << " ------ " << endl;
  xxr=imag(xx);
  plot(xxr, wait=1);

{
  int N=10;
  complex [int,int] A(N,N);
  complex [int] a(N),b(N),bb(N);
  A =0;
  for (int i=0;i<N;i++)
    {
      A(i,i)=1.+i;
      if(i+1 < N)    A(i,i+1)=-i-1i*i;
      a[i]=i*(1.+2i);
    }
  b=A*a;

  cout << " b =" << b << endl ;
  cout << " a =" << a << endl ;
  cout << " b'*b (hermissian product) = " << b'*b << endl;
  cout << " a'*a = " << a'*a << endl;
  assert( abs(imag(b'*b)) <1e-5);
  cout << "xxxx\n"; 
  matrix<complex> sparseA=A;

  cout << sparseA << endl;
  sparseA = 2*sparseA+sparseA;
  sparseA = 4*sparseA+sparseA*(5+1i); //  * 27
  matrix<complex> sparseB=sparseA;//+sparseA+sparseA; ;
  cout << sparseA << endl;
  cout << sparseB << endl; // *81 
  cout << "sparseB = " << sparseB(0,0) << endl;
  // ajoute version  2.0-2
  sparseA=A;
  verbosity=4;
  if(HaveUMFPACK)		
    set(sparseA,solver=UMFPACK,tolpivot=1e-10,tolpivotsym=1e-9);  
  else 
    set(sparseA,solver=GMRES);  
  bb=sparseA^-1*a;
  verbosity=1;
  b = sparseA*bb;
  b -= a;
  cout << " nb coef sparseA " << sparseA.nbcoef << endl; 
  cout << " ||b.||_1  " << b.l1 << endl;
  cout << " ||b.||_2  " << b.l2 << endl;
  cout << " ||b.||_infty  " << b.linfty << endl;
  assert(b.l1 < 1e-10);
}