C++ SparseMatrix::apply方法代码示例

//.........这里部分代码省略.........


   cout << "- writing A to the file stiffmat.m ..." << endl;
   std::ofstream Astream("stiffmat.m");
   Astream << "A=";
   print_matrix(A, Astream);
   Astream << ";" << endl;
   Astream.close();
   cout << "  ... done!" << endl;

  cout << "- set up right-hand side..." << endl;
  Vector<double> b;
  setup_righthand_side(eq, Lambda, b);
  
  
  cout << "- writing b to the file rhs.m ..." << endl;
   std::ofstream bstream("rhs.m");
   bstream << "b=";
   print_vector(b, bstream);
   bstream << ";" << endl;
   bstream.close();
   cout << "  ... done!" << endl;

  
  b.compress(1e-14);
  

  Vector<double> x(Lambda.size()), x2(Lambda.size()), err(Lambda.size()), err2(Lambda.size()); x = 0, x2 = 0;
  unsigned int iterations;
  
  
CG(A, b, x, 1e-8, 5000, iterations);

A.apply(x, err);
  //cout << "Ax " << err << endl << endl;
  
    
  err -= b;
  //cout << "residual " << err << endl;
  cout << " residual (infinity) norm " << linfty_norm(err) << endl;
  
#if 1
  
  // evaluate approximate solution on a grid
    Vector<double> u_values(N+1);
    for (unsigned int i = 0; i <= N; i++) {
      const double point = i*h;
      int id = 0;
      for (set<Index>::const_iterator it(Lambda.begin()); it != Lambda.end(); ++it, ++id) 
	u_values[i] += x[id] * basis.evaluate(0, *it, point)*1./eq.D(*it);
    }
     

    
    
    
    
    
     

    // compute some error-norms
    const double Linfty_error = linfty_norm(u_values-uexact_values);
    cout << "  L_infinity error on a subgrid: " << Linfty_error << endl;
    
    
    const double L2_error = sqrt(l2_norm_sqr(u_values-uexact_values)*h);

int main(int argc, char** argv) {
    
    const int d  = 3;
    const int dT = 3;
    const int jmax = 10;
    const bool normalization = 0;//choose 1 for Laplacian
    
    const unsigned int testcase=5;
    PeriodicTestProblem<testcase> tper;
    Function<1>* uexact = 0;
    switch(testcase) {
        case 1:
            uexact = new Function1();
            break;
        case 2:
            uexact = new Function2();
            break;
        case 3:
            uexact = new Function3();
            break;
        case 4:
            uexact = new Function4();
            break;
        case 5:
            uexact = new Hat();    
        default:
            break;
    }
    
    
    
#ifdef PERIODIC_CDFBASIS
    typedef CDFBasis<d,dT> RBasis;
    typedef PeriodicBasis<RBasis> Basis;
    typedef Basis::Index Index;
    Basis basis;
    basis.set_jmax(jmax);
    
    typedef PeriodicIntervalGramian<RBasis> Problem;
    Problem G(tper, basis);
    const int j0= G.basis().j0();
    
#if 1
    //Plot of wavelet with coefficient mu
    Index mu(3,0,1);
    SampledMapping<1> sm1(basis.evaluate(mu, 8, normalization));
    std::ofstream u_stream1("plot_periodicwavelet.m");
    sm1.matlab_output(u_stream1);
    u_stream1 << "figure;\nplot(x,y);"
            << "title('periodic wavelet/generator with index: " << mu << "');" << endl;
    u_stream1.close();
#endif
    
    set<Index> Lambda;
    Vector<double> x;
    
  for (int j = j0; j <= jmax; j++) {
        Lambda.clear();
        
    cout << "  j=" << j << ":" << endl;
    //Implementation of the index set
    
    for (Index lambda = G.basis().first_generator(j0);; ++lambda) {
        Lambda.insert(lambda);
        if (lambda == G.basis().last_wavelet(j)) break;    
    }
    SparseMatrix<double> A;
    cout << "- set up stiffness matrix..." << endl;
    setup_stiffness_matrix(G, Lambda, A);
//    
    
    
    cout << "- set up right-hand side..." << endl;
    Vector<double> b;
    setup_righthand_side(G, Lambda, b);
//    cout << "- right hand side: " << b << endl << endl;
    
    
    x.resize(Lambda.size()); x = 0;
    Vector<double> residual(Lambda.size()); 
    unsigned int iterations;
    
    CG(A, b, x, 1e-15, 250, iterations);
    cout << "  Galerkin system solved with residual (infinity) norm ";
    A.apply(x, residual);
    residual -= b;
    cout << linfty_norm(residual)
	 << " (" << iterations << " iterations needed)" << endl;
    
    // evaluate approximate solution on a grid
    const unsigned int N = 100;
    const double h = 1./N;
    Vector<double> u_values(N+1);
    for (unsigned int i = 0; i <= N; i++) {
      const double point = i*h;
      int id = 0;
      for (set<Index>::const_iterator it(Lambda.begin()); it != Lambda.end(); ++it, ++id) 
	u_values[i] += x[id] * basis.evaluate(0, *it, point, normalization);
    }
//     cout << "  point values of Galerkin solution: " << u_values << endl;
//.........这里部分代码省略.........