Method MultiMult

linalg/cg.cpp:133–249

  template <class IPTYPE>
  void CGSolver<IPTYPE> :: MultiMult (const BaseVector & f, BaseVector & u, const int dim) const
  {
    try
      {
	// Solve A u = f
        BaseStatusHandler::SetThreadPercentage(0);

	auto d = f.CreateVector();
	auto w = f.CreateVector();
	auto s = f.CreateVector();

	int n = 0;
	Vector<SCAL> al(dim), be(dim), wd(dim), wdn(dim), kss(dim);
	double err;

	if (initialize)
	  {
	    u = 0.0;
	    d = f;
	  }
	else
	  {
	    d = f - (*a) * u;
	  }
	if (c)
	  w = (*c) * d;
	else
	  w = d;

	s = w;
	
	BruteInnerProduct(w,d,wdn);	 

	if (printrates) cout << IM(1) << "0 " << sqrt(L2Norm(wdn)) << endl;
	if (L2Norm(wdn) == 0.0) wdn = 1;	

	if(stop_absolute)
	  err = prec * prec;
	else
	  err = prec * prec * L2Norm (wdn);
	
	double lwstart = log(L2Norm(wdn));
	double lerr = log(err);
	

	while (n++ < maxsteps && L2Norm(wdn) > err && !(BaseStatusHandler::ShouldTerminate()))
	  {
	    w = (*a) * s;

	    wd = wdn;

	    BruteInnerProduct(s,w,kss);
	   
	    //(*testout) << "INNERPROD kss " <<kss << endl;
	    if (L2Norm(kss) == 0.0) break;
	    
	    for(int i = 0; i<dim; i++)
	      al[i] = wd[i] / kss[i];