Method ComputeMatrices

linalg/commutingAMG.cpp:233–274

  

  void AMG_H1 :: ComputeMatrices (const BaseSparseMatrix & mat)
  {
    pmat = &mat;

    if (0) // if (prol)
      {
	// define smoothing blocks
	int nv = prol->Width();
	Array<int> cnt(nv);
	cnt = 0;
	for (int i = 0; i < prol->Height(); i++)
	  for (int j = 0; j < prol->GetRowIndices(i).Size(); j++)
	    cnt[prol->GetRowIndices(i)[j]]++;
	
	Table<int> smblocks(cnt);
	cnt = 0;
	for (int i = 0; i < prol->Height(); i++)
	  for (int j = 0; j < prol->GetRowIndices(i).Size(); j++)
	    {
	      int jj = prol->GetRowIndices(i)[j];
	      smblocks[jj][cnt[jj]] = i;
	      cnt[jj]++;
	    }
	// *testout << "prol = " << endl << *prol << endl;
	// *testout << "smoothing blocks = " << endl << *smblocks << endl;
	
	bjacobi = mat.CreateBlockJacobiPrecond (make_shared<Table<int>> (smblocks));
      }
    jacobi = mat.CreateJacobiPrecond ();

    if (recAMG)
      {
	coarsemat = shared_ptr<BaseSparseMatrix>(mat.Restrict (*prol));
	recAMG -> ComputeMatrices (*coarsemat);
	inv = 0;
      }
    else 
      {
	mat.SetInverseType ("sparsecholesky");
        inv = mat.InverseMatrix();
      }
  }

  size_t AMG_H1 :: NZE() const
  {