Method BlockJacobiPrecond

linalg/blockjacobi.cpp:317–582

  ///
  template <class TM, class TV_ROW, class TV_COL>
  BlockJacobiPrecond<TM, TV_ROW, TV_COL> ::
  BlockJacobiPrecond (shared_ptr<const SparseMatrix<TM,TV_ROW,TV_COL>> amat, 
		      shared_ptr<Table<int>> ablocktable, bool cumulate_block_diags)
    : BaseBlockJacobiPrecond(ablocktable), mat(amat), 
      invdiag(ablocktable->Size())
  {
    static Timer t("BlockJacobiPrecond ctor"); RegionTimer reg(t);
    static Timer tinv("BlockJacobiPrecond ctor inv");
    static Timer tget("BlockJacobiPrecond ctor get");
    static Timer tprep("BlockJacobiPrecond ctor prep");
    static Timer tpar("BlockJacobiPrecond ctor par");
    cout << IM(3) << "BlockJacobi Preconditioner, constructor called, #blocks = " << blocktable->Size() << endl;


    // double prevtime = WallTime();

    
    // find nze element in all blocks together 
    /*
    nze = 0;
    for (auto block : *blocktable)
      for (auto row : block)
	nze += amat.GetRowIndices(row).Size();
    */
    nze = 
      ParallelReduce (blocktable->Size(),
                      [&] (size_t i)
                      {
                        size_t nze = 0;
                        for (auto row : (*blocktable)[i])
                          nze += amat->GetRowIndices(row).Size();
                        return nze;
                      },
                      [] (size_t a, size_t b) { return a+b; },
                      size_t(0));

    /*
    size_t totmem = 0;
    for (auto i : Range (*blocktable))
      totmem += sqr ((*blocktable)[i].Size());
    */
    size_t totmem = 
      ParallelReduce (blocktable->Size(),
                      [&] (size_t i) { return sqr ((*blocktable)[i].Size()); },
                      [] (size_t a, size_t b) { return a+b; },
                      size_t(0));

    bigmem.SetSize(totmem);

    cout << IM(5) << "avg entrysize:   " << blocktable->AsArray().Size()/blocktable->Size() << endl;
    cout << IM(5) << "avg entrysize^2: " << totmem/blocktable->Size() << endl;    
    
    totmem = 0;
    for (auto i : Range (*blocktable))
      {
        size_t bs = (*blocktable)[i].Size();
        new ( & invdiag[i] ) FlatMatrix<TM> (bs, bs, bigmem.Addr(totmem));
        totmem += sqr (bs);