MGPoissonSolver.cpp 19.2 KB
Newer Older
1
#ifdef HAVE_SAAMG_SOLVER
2
//#define DBG_STENCIL
3

kraus's avatar
kraus committed
4
#include "Solvers/MGPoissonSolver.h"
5 6 7 8 9 10 11 12

#include "Structure/BoundaryGeometry.h"
#include "ArbitraryDomain.h"
#include "EllipticDomain.h"
#include "BoxCornerDomain.h"
//#include "RectangularDomain.h"

#include "Algorithms/PartBunch.h"
13
#include "Track/Track.h"
gsell's avatar
gsell committed
14 15 16 17 18 19 20 21
#include "Physics/Physics.h"
#include "Attributes/Attributes.h"
#include "ValueDefinitions/RealVariable.h"
#include "AbstractObjects/OpalData.h"
#include "Utilities/OpalException.h"
#include "Attributes/Attributes.h"
#include "ValueDefinitions/RealVariable.h"
#include "AbstractObjects/OpalData.h"
kraus's avatar
kraus committed
22
#include "Utilities/OpalOptions.h"
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-local-typedefs"
#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_Operator.h"
#include "EpetraExt_RowMatrixOut.h"
#include "Epetra_Import.h"

#include "Teuchos_CommandLineProcessor.hpp"

#include "BelosLinearProblem.hpp"
#include "BelosRCGSolMgr.hpp"
#include "BelosEpetraAdapter.hpp"
#include "BelosBlockCGSolMgr.hpp"

#include "ml_MultiLevelPreconditioner.h"
#include "ml_MultiLevelOperator.h"
#include "ml_epetra_utils.h"

#include "Isorropia_Exception.hpp"
#include "Isorropia_Epetra.hpp"
#include "Isorropia_EpetraRedistributor.hpp"
#include "Isorropia_EpetraPartitioner.hpp"

#pragma GCC diagnostic pop

gsell's avatar
gsell committed
51
#include <algorithm>
52
#include <omp.h>
gsell's avatar
gsell committed
53

54 55
using Teuchos::RCP;
using Teuchos::rcp;
gsell's avatar
gsell committed
56 57
using Physics::c;

58
MGPoissonSolver::MGPoissonSolver ( PartBunch &beam,
59 60 61 62 63 64 65
                                   Mesh_t *mesh,
                                   FieldLayout_t *fl,
                                   std::vector<BoundaryGeometry *> geometries,
                                   std::string itsolver,
                                   std::string interpl,
                                   double tol,
                                   int maxiters,
66
                                   std::string precmode):
67
    itsBunch_m(&beam),
68 69 70 71 72
    mesh_m(mesh),
    layout_m(fl),
    geometries_m(geometries),
    tol_m(tol),
    maxiters_m(maxiters),
73
    Comm(Ippl::getComm()) {
74

gsell's avatar
gsell committed
75 76
    domain_m = layout_m->getDomain();

77
    for (int i = 0; i < 3; i++) {
gsell's avatar
gsell committed
78 79 80 81 82 83 84 85 86 87 88 89 90
        hr_m[i] = mesh_m->get_meshSpacing(i);
        orig_nr_m[i] = domain_m[i].length();
    }

    if(itsolver == "CG") itsolver_m = AZ_cg;
    else if(itsolver == "BiCGStab") itsolver_m = AZ_bicgstab;
    else if(itsolver == "GMRES") itsolver_m = AZ_gmres;
    else throw OpalException("MGPoissonSolver", "No valid iterative solver selected!");

    precmode_m = STD_PREC;
    if(precmode == "STD") precmode_m = STD_PREC;
    else if(precmode == "HIERARCHY") precmode_m = REUSE_HIERARCHY;
    else if(precmode == "REUSE") precmode_m = REUSE_PREC;
91
    else if(precmode == "NO") precmode_m = NO;
gsell's avatar
gsell committed
92 93 94 95 96 97

    tolerableIterationsCount_m = 2;
    numIter_m = -1;
    forcePreconditionerRecomputation_m = false;

    hasParallelDecompositionChanged_m = true;
98
    repartFreq_m = 1000; 
gsell's avatar
gsell committed
99
    useRCB_m = false;
100
    if(Ippl::Info->getOutputLevel() > 1)
gsell's avatar
gsell committed
101 102 103 104
        verbose_m = true;
    else
        verbose_m = false;

105
    // Find CURRENT geometry
gsell's avatar
gsell committed
106 107
    currentGeometry = geometries_m[0];
    if(currentGeometry->getFilename() == "") {
108
        if(currentGeometry->getTopology() == "ELLIPTIC"){
109
            bp = new EllipticDomain(currentGeometry, orig_nr_m, hr_m, interpl);
110
	} else if (currentGeometry->getTopology() == "BOXCORNER") {
111
            bp = new BoxCornerDomain(currentGeometry->getA(), currentGeometry->getB(), currentGeometry->getC(), currentGeometry->getLength(),currentGeometry->getL1(), currentGeometry->getL2(), orig_nr_m, hr_m, interpl);
112
            bp->Compute(itsBunch_m->get_hr());
gsell's avatar
gsell committed
113
        } else {
114 115
            throw OpalException("MGPoissonSolver::MGPoissonSolver",
                                "Geometry not known");
gsell's avatar
gsell committed
116
        }
117
    } else 
118 119
	bp = new ArbitraryDomain(currentGeometry, orig_nr_m, hr_m, interpl);

gsell's avatar
gsell committed
120 121 122 123
    Map = 0;
    A = Teuchos::null;
    LHS = Teuchos::null;
    RHS = Teuchos::null;
124
    MLPrec = 0;
125
    prec_m = Teuchos::null;
126

gsell's avatar
gsell committed
127 128
    numBlocks_m = Options::numBlocks;
    recycleBlocks_m = Options::recycleBlocks;
129
    nLHS_m = Options::nLHS;
gsell's avatar
gsell committed
130 131
    SetupMLList();
    SetupBelosList();
132
    problem_ptr = rcp(new Belos::LinearProblem<ST, MV, OP>);
133
    // setup Belos solver
gsell's avatar
gsell committed
134
    if(numBlocks_m == 0 || recycleBlocks_m == 0)
135
        solver_ptr = rcp(new Belos::BlockCGSolMgr<double, MV, OP>());
gsell's avatar
gsell committed
136
    else
137 138
        solver_ptr = rcp(new Belos::RCGSolMgr<double, MV, OP>());
    solver_ptr->setParameters(rcp(&belosList, false));
gsell's avatar
gsell committed
139 140 141 142
    convStatusTest = rcp(new Belos::StatusTestGenResNorm<ST, MV, OP> (tol));
    convStatusTest->defineScaleForm(Belos::NormOfRHS, Belos::TwoNorm);

    //all timers used here
143 144 145 146 147 148 149 150
    FunctionTimer1_m = IpplTimings::getTimer("BGF-IndexCoordMap");
    FunctionTimer2_m = IpplTimings::getTimer("computeMap");
    FunctionTimer3_m = IpplTimings::getTimer("IPPL to RHS");
    FunctionTimer4_m = IpplTimings::getTimer("ComputeStencil");
    FunctionTimer5_m = IpplTimings::getTimer("ML");
    FunctionTimer6_m = IpplTimings::getTimer("Setup");
    FunctionTimer7_m = IpplTimings::getTimer("CG");
    FunctionTimer8_m = IpplTimings::getTimer("LHS to IPPL");
gsell's avatar
gsell committed
151 152 153 154 155 156 157 158 159 160
}


MGPoissonSolver::MGPoissonSolver(PartBunch &beam):
    layout_m(&beam.getFieldLayout()), mesh_m(&beam.getMesh()), itsBunch_m(&beam), LHS(0), Comm(Ippl::getComm()) {

    throw OpalException("MGPoissonSolver", "Copy Constructor not implemented yet");
}

MGPoissonSolver::~MGPoissonSolver() {
161 162 163 164 165 166 167 168 169 170 171 172 173 174
    if (Map) delete Map; Map = 0;
    if (MLPrec) delete MLPrec; MLPrec = 0;
    A = Teuchos::null;
    LHS = Teuchos::null;
    RHS = Teuchos::null;
    prec_m = Teuchos::null;
    solver_ptr = Teuchos::null;
    problem_ptr = Teuchos::null;
}

void MGPoissonSolver::deletePtr() {
    A = Teuchos::null;
    LHS = Teuchos::null;
    RHS = Teuchos::null;
175 176
    if(Map) delete Map;
    if(MLPrec) delete MLPrec; MLPrec=0;
177
    prec_m = Teuchos::null;
gsell's avatar
gsell committed
178 179 180 181 182 183
}

void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr, double zshift) {
    throw OpalException("MGPoissonSolver", "not implemented yet");
}

184
void MGPoissonSolver::computeMap(NDIndex<3> localId) {
185 186
    if (itsBunch_m->getLocalTrackStep()%repartFreq_m == 0){
        deletePtr(); 
gsell's avatar
gsell committed
187
        if(useRCB_m)
188
            redistributeWithRCB(localId);
189
        else
190
            IPPLToMap3D(localId);
191 192

        extrapolateLHS();
193 194
    }
}
gsell's avatar
gsell committed
195

196 197 198 199 200
void MGPoissonSolver::extrapolateLHS() {
// Aitken-Neville
// Pi0 (x) := yi , i = 0 : n
// Pik (x) := (x − xi ) Pi+1,k−1(x) − (x − xi+k ) Pi,k−1(x) /(xi+k − xi )
// k = 1, . . . , n, i = 0, . . . , n − k.
201
    //we also have to redistribute LHS
202

203 204 205 206 207 208 209 210 211
    if (Teuchos::is_null(LHS)){
        LHS = rcp(new Epetra_Vector(*Map));
        LHS->PutScalar(1.0);
    } else {
        RCP<Epetra_Vector> tmplhs = rcp(new Epetra_Vector(*Map));
        Epetra_Import importer(*Map, LHS->Map());
        tmplhs->Import(*LHS, importer, Add);
        LHS = tmplhs;
    }
gsell's avatar
gsell committed
212

213
    //...and all previously saved LHS
214
    std::deque< Epetra_Vector >::iterator it = OldLHS.begin();
215 216 217 218 219 220 221 222 223 224
    if(OldLHS.size() > 0) {
        int n = OldLHS.size();
        for(int i = 0; i < n; ++i) {
            Epetra_Vector tmplhs = Epetra_Vector(*Map);
            Epetra_Import importer(*Map, it->Map());
            tmplhs.Import(*it, importer, Add);
            *it = tmplhs;
            ++it;
        }
    }
gsell's avatar
gsell committed
225

226 227 228
    // extrapolate last OldLHS.size LHS to get a new start vector
    it = OldLHS.begin();
    if(nLHS_m == 0 || OldLHS.size()==0)
229 230 231
        LHS->PutScalar(1.0);
    else if(OldLHS.size() == 1)
        *LHS = *it;
232 233 234
    else if(OldLHS.size() == 2)
        LHS->Update(2.0, *it++, -1.0, *it, 0.0);
    else if(OldLHS.size() > 2) {
235
        int n = OldLHS.size();
236 237 238
        P = rcp(new Epetra_MultiVector(*Map, nLHS_m, false));
        for(int i = 0; i < n; ++i) {
           *(*P)(i) = *it++;
239
        }
240 241 242 243
        for(int k = 1; k < n; ++k) {
           for(int i = 0; i < n - k; ++i) {
              (*P)(i)->Update(-(i + 1) / (float)k, *(*P)(i + 1), (i + k + 1) / (float)k);
           }
244 245
        }
        *LHS = *(*P)(0);
246 247
     } else
        throw OpalException("MGPoissonSolver", "Invalid number of old LHS: " + OldLHS.size());
gsell's avatar
gsell committed
248 249
}

250

gsell's avatar
gsell committed
251 252 253 254 255
// given a charge-density field rho and a set of mesh spacings hr,
// compute the electric field and put in eg by solving the Poisson's equation
// XXX: use matrix stencil in computation directly (no Epetra, define operators
// on IPPL GRID)
void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
256

gsell's avatar
gsell committed
257 258 259
    nr_m[0] = orig_nr_m[0];
    nr_m[1] = orig_nr_m[1];
    nr_m[2] = orig_nr_m[2];
260

261 262
    bp->setGlobalMeanR(itsBunch_m->getGlobalMeanR());
    bp->setGlobalToLocalQuaternion(itsBunch_m->getGlobalToLocalQuaternion());
gsell's avatar
gsell committed
263
    bp->setNr(nr_m);
264 265
    NDIndex<3> localId = layout_m->getLocalNDIndex();

266
    IpplTimings::startTimer(FunctionTimer1_m);
267
    if(!itsBunch_m->getLocalTrackStep())
268
        bp->Compute(hr, localId);
269
    IpplTimings::stopTimer(FunctionTimer1_m);
270 271

    // Define the Map
272
    IpplTimings::startTimer(FunctionTimer2_m);
273
    computeMap(localId);
274
    IpplTimings::stopTimer(FunctionTimer2_m);
gsell's avatar
gsell committed
275

276
    // Allocate the RHS with the new Epetra Map
277 278
    if (Teuchos::is_null(RHS)) 
        RHS = rcp(new Epetra_Vector(*Map));
279
    RHS->PutScalar(0.0);
280 281 282
    // get charge densities from IPPL field and store in Epetra vector (RHS)
    IpplTimings::startTimer(FunctionTimer3_m);
    int id = 0;
283
    float scaleFactor = itsBunch_m->getdT();
284
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
285
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
286
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
287
                  if (bp->isInside(idx, idy, idz))
288
                    RHS->Values()[id++] = 4*M_PI*rho[idx][idy][idz].get()/scaleFactor;
gsell's avatar
gsell committed
289 290 291
            }
        }
    }
292
    IpplTimings::stopTimer(FunctionTimer3_m);
293

gsell's avatar
gsell committed
294
    // build discretization matrix
295
    IpplTimings::startTimer(FunctionTimer4_m);
296
    if(Teuchos::is_null(A))
297
        A = rcp(new Epetra_CrsMatrix(Copy, *Map,  7, true));
298
    ComputeStencil(hr, RHS);
299
    IpplTimings::stopTimer(FunctionTimer4_m);
300

gsell's avatar
gsell committed
301 302 303
#ifdef DBG_STENCIL
    EpetraExt::RowMatrixToMatlabFile("DiscrStencil.dat", *A);
#endif
304

305
    IpplTimings::startTimer(FunctionTimer5_m);
306
    if(!MLPrec) {
307 308 309 310
        MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList_m);
    } else if (precmode_m == REUSE_HIERARCHY) {
        MLPrec->ReComputePreconditioner();
    } else if (precmode_m == REUSE_PREC){
311
    }
312
    IpplTimings::stopTimer(FunctionTimer5_m);
313

gsell's avatar
gsell committed
314 315
    // setup preconditioned iterative solver
    // use old LHS solution as initial guess
316
    IpplTimings::startTimer(FunctionTimer6_m);
317 318 319 320

    problem_ptr->setOperator(A);
    problem_ptr->setLHS(LHS);
    problem_ptr->setRHS(RHS);
321 322
    if(Teuchos::is_null(prec_m))
        prec_m = Teuchos::rcp ( new Belos::EpetraPrecOp ( rcp(MLPrec,false)));
323 324 325 326
    problem_ptr->setLeftPrec(prec_m);
    solver_ptr->setProblem( problem_ptr);
    if (!problem_ptr->isProblemSet()) {
        if (problem_ptr->setProblem() == false) {
327
            ERRORMSG("Belos::LinearProblem failed to set up correctly!" << endl);
328
        }
gsell's avatar
gsell committed
329
    }
330
    IpplTimings::stopTimer(FunctionTimer6_m);
gsell's avatar
gsell committed
331 332
    std::ofstream timings;
    char filename[50];
333
    sprintf(filename, "timing_MX%d_MY%d_MZ%d_nProc%d_recB%d_numB%d_nLHS%d", orig_nr_m[0], orig_nr_m[1], orig_nr_m[2], Comm.NumProc(), recycleBlocks_m, numBlocks_m, nLHS_m);
gsell's avatar
gsell committed
334 335 336
    if(Comm.MyPID() == 0) timings.open(filename, std::ios::app);
    double time = MPI_Wtime();

337
    IpplTimings::startTimer(FunctionTimer7_m);
338
    	solver_ptr->solve();
339
    IpplTimings::stopTimer(FunctionTimer7_m);
gsell's avatar
gsell committed
340 341 342 343 344 345 346 347

    time = MPI_Wtime() - time;
    double minTime, maxTime, avgTime;
    Comm.MinAll(&time, &minTime, 1);
    Comm.MaxAll(&time, &maxTime, 1);
    Comm.SumAll(&time, &avgTime, 1);
    avgTime /= Comm.NumProc();
    if(Comm.MyPID() == 0) timings <<
348
                                      solver_ptr->getNumIters() << "\t" <<
gsell's avatar
gsell committed
349 350 351 352 353 354
                                      //time <<" "<<
                                      minTime << "\t" <<
                                      maxTime << "\t" <<
                                      avgTime << "\t" <<
                                      numBlocks_m << "\t" <<
                                      recycleBlocks_m << "\t" <<
355
                                      nLHS_m << "\t" <<
gsell's avatar
gsell committed
356 357 358 359 360
                                      //OldLHS.size() <<"\t"<<
                                      std::endl;
    if(Comm.MyPID() == 0) timings.close();

    // Store new LHS in OldLHS
361 362
    if(nLHS_m > 1) OldLHS.push_front(*(LHS.get()));
    if(OldLHS.size() > nLHS_m) OldLHS.pop_back();
gsell's avatar
gsell committed
363 364

    //now transfer solution back to IPPL grid
365 366
    IpplTimings::startTimer(FunctionTimer8_m);
    id = 0;
gsell's avatar
gsell committed
367
    rho = 0.0;
368
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
369
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
370
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
371 372
                  NDIndex<3> l(Index(idx, idx), Index(idy, idy), Index(idz, idz));
                  if (bp->isInside(idx, idy, idz))
373
                     rho.localElement(l) = LHS->Values()[id++]*scaleFactor;
gsell's avatar
gsell committed
374 375 376
            }
        }
    }
377
    IpplTimings::stopTimer(FunctionTimer8_m);
378 379 380 381 382 383 384

    if(itsBunch_m->getLocalTrackStep()+1 == (long long)Track::block->localTimeSteps.front()) {
        A = Teuchos::null;
        LHS = Teuchos::null;
        RHS = Teuchos::null;
        prec_m = Teuchos::null;
    }	
gsell's avatar
gsell committed
385 386 387
}


388
void MGPoissonSolver::redistributeWithRCB(NDIndex<3> localId) {
gsell's avatar
gsell committed
389 390 391

    int numMyGridPoints = 0;

392
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
393
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
394
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
395
                 if (bp->isInside(idx, idy, idz))
gsell's avatar
gsell committed
396
                    numMyGridPoints++;
397
             }
gsell's avatar
gsell committed
398
        }
399
     }
gsell's avatar
gsell committed
400 401 402 403 404 405 406 407 408 409

    Epetra_BlockMap bmap(-1, numMyGridPoints, 1, 0, Comm);
    Teuchos::RCP<const Epetra_MultiVector> coords = Teuchos::rcp(new Epetra_MultiVector(bmap, 3, false));

    double *v;
    int stride, stride2;

    coords->ExtractView(&v, &stride);
    stride2 = 2 * stride;

410
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
411
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
412
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
413 414 415 416
                 if (bp->isInside(idx, idy, idz)) {
                    v[0] = (double)idx;
                    v[stride] = (double)idy;
                    v[stride2] = (double)idz;
gsell's avatar
gsell committed
417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445
                    v++;
                }
            }
        }
    }

    Teuchos::ParameterList paramlist;
    paramlist.set("Partitioning Method", "RCB");
    Teuchos::ParameterList &sublist = paramlist.sublist("ZOLTAN");
    sublist.set("RCB_RECTILINEAR_BLOCKS", "1");
    sublist.set("DEBUG_LEVEL", "1");

    Teuchos::RCP<Isorropia::Epetra::Partitioner> part = Teuchos::rcp(new Isorropia::Epetra::Partitioner(coords, paramlist));
    Isorropia::Epetra::Redistributor rd(part);
    Teuchos::RCP<Epetra_MultiVector> newcoords = rd.redistribute(*coords);

    newcoords->ExtractView(&v, &stride);
    stride2 = 2 * stride;
    numMyGridPoints = 0;
    std::vector<int> MyGlobalElements;
    for(int i = 0; i < newcoords->MyLength(); i++) {
        MyGlobalElements.push_back(bp->getIdx(v[0], v[stride], v[stride2]));
        v++;
        numMyGridPoints++;
    }

    Map = new Epetra_Map(-1, numMyGridPoints, &MyGlobalElements[0], 0, Comm);
}

446
void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
gsell's avatar
gsell committed
447
    int NumMyElements = 0;
448
    std::vector<int> MyGlobalElements;
gsell's avatar
gsell committed
449

450
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
451
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
452
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
453
                if (bp->isInside(idx, idy, idz)) {
454
                    MyGlobalElements.push_back(bp->getIdx(idx, idy, idz));
gsell's avatar
gsell committed
455 456 457
                    NumMyElements++;
                }
            }
458 459
	}
    }
gsell's avatar
gsell committed
460 461 462 463 464
    Map = new Epetra_Map(-1, NumMyElements, &MyGlobalElements[0], 0, Comm);
}

void MGPoissonSolver::ComputeStencil(Vector_t hr, Teuchos::RCP<Epetra_Vector> RHS) {

465 466
    A->PutScalar(0.0);

gsell's avatar
gsell committed
467 468 469
    int NumMyElements = Map->NumMyElements();
    int *MyGlobalElements = Map->MyGlobalElements();

470 471
    std::vector<double> Values(6);
    std::vector<int> Indices(6);
gsell's avatar
gsell committed
472

473
    for(int i = 0 ; i < NumMyElements ; i++) {
gsell's avatar
gsell committed
474 475 476

        int NumEntries = 0;

477
        double WV, EV, SV, NV, FV, BV, CV, scaleFactor=1.0;
478
        int W, E, S, N, F, B;
gsell's avatar
gsell committed
479

480
        bp->getBoundaryStencil(MyGlobalElements[i], WV, EV, SV, NV, FV, BV, CV, scaleFactor);
481
        RHS->Values()[i] *= scaleFactor;
gsell's avatar
gsell committed
482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508

        bp->getNeighbours(MyGlobalElements[i], W, E, S, N, F, B);
        if(E != -1) {
            Indices[NumEntries] = E;
            Values[NumEntries++] = EV;
        }
        if(W != -1) {
            Indices[NumEntries] = W;
            Values[NumEntries++] = WV;
        }
        if(S != -1) {
            Indices[NumEntries] = S;
            Values[NumEntries++] = SV;
        }
        if(N != -1) {
            Indices[NumEntries] = N;
            Values[NumEntries++] = NV;
        }
        if(F != -1) {
            Indices[NumEntries] = F;
            Values[NumEntries++] = FV;
        }
        if(B != -1) {
            Indices[NumEntries] = B;
            Values[NumEntries++] = BV;
        }

509
    	// if matrix has already been filled (FillComplete()) we can only
gsell's avatar
gsell committed
510
        // replace entries
511

512 513 514 515 516
        if(A->Filled()) {
            // off-diagonal entries
            A->ReplaceGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
            // diagonal entry
            A->ReplaceGlobalValues(MyGlobalElements[i], 1, &CV, MyGlobalElements + i);
517
        } else {
gsell's avatar
gsell committed
518 519 520
            // off-diagonal entries
            A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
            // diagonal entry
521
            A->InsertGlobalValues(MyGlobalElements[i], 1, &CV, MyGlobalElements + i);
522
        }
gsell's avatar
gsell committed
523 524 525
    }

    A->FillComplete();
526

gsell's avatar
gsell committed
527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549
    A->OptimizeStorage();
}

void MGPoissonSolver::printLoadBalanceStats() {

    //compute some load balance statistics
    size_t myNumPart = Map->NumMyElements();
    size_t NumPart = Map->NumGlobalElements() * 1.0 / Comm.NumProc();
    double imbalance = 1.0;
    if(myNumPart >= NumPart)
        imbalance += (myNumPart - NumPart) / NumPart;
    else
        imbalance += (NumPart - myNumPart) / NumPart;

    double max = 0.0, min = 0.0, avg = 0.0;
    int minn = 0, maxn = 0;
    MPI_Reduce(&imbalance, &min, 1, MPI_DOUBLE, MPI_MIN, 0, Ippl::getComm());
    MPI_Reduce(&imbalance, &max, 1, MPI_DOUBLE, MPI_MAX, 0, Ippl::getComm());
    MPI_Reduce(&imbalance, &avg, 1, MPI_DOUBLE, MPI_SUM, 0, Ippl::getComm());
    MPI_Reduce(&myNumPart, &minn, 1, MPI_INT, MPI_MIN, 0, Ippl::getComm());
    MPI_Reduce(&myNumPart, &maxn, 1, MPI_INT, MPI_MAX, 0, Ippl::getComm());

    avg /= Comm.NumProc();
550 551
    if(Comm.MyPID() == 0) *gmsg << "LBAL min = " << min << ", max = " << max << ", avg = " << avg << endl;
    if(Comm.MyPID() == 0) *gmsg << "min nr gridpoints = " << minn << ", max nr gridpoints = " << maxn << endl;
gsell's avatar
gsell committed
552 553 554 555 556 557
}

Inform &MGPoissonSolver::print(Inform &os) const {
    os << "* *************** M G P o i s s o n S o l v e r ************************************ " << endl;
    os << "* h " << hr_m << '\n';
    os << "* ********************************************************************************** " << endl;
558
    return os;
gsell's avatar
gsell committed
559 560
}

561
#endif /* HAVE_SAAMG_SOLVER */