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 266
    NDIndex<3> localId = layout_m->getLocalNDIndex();

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

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

277
    // Allocate the RHS with the new Epetra Map
278 279
    if (Teuchos::is_null(RHS)) 
        RHS = rcp(new Epetra_Vector(*Map));
280
    RHS->PutScalar(0.0);
281

282 283 284
    // get charge densities from IPPL field and store in Epetra vector (RHS)
    IpplTimings::startTimer(FunctionTimer3_m);
    int id = 0;
285
    float scaleFactor = itsBunch_m->getdT();
286
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
287
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
288
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
289
                  if (bp->isInside(idx, idy, idz))
290
                    RHS->Values()[id++] = 4*M_PI*rho[idx][idy][idz].get()/scaleFactor;
gsell's avatar
gsell committed
291 292 293
            }
        }
    }
294
    IpplTimings::stopTimer(FunctionTimer3_m);
295

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

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

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

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

    problem_ptr->setOperator(A);
    problem_ptr->setLHS(LHS);
    problem_ptr->setRHS(RHS);
323 324
    if(Teuchos::is_null(prec_m))
        prec_m = Teuchos::rcp ( new Belos::EpetraPrecOp ( rcp(MLPrec,false)));
325 326 327 328
    problem_ptr->setLeftPrec(prec_m);
    solver_ptr->setProblem( problem_ptr);
    if (!problem_ptr->isProblemSet()) {
        if (problem_ptr->setProblem() == false) {
329
            ERRORMSG("Belos::LinearProblem failed to set up correctly!" << endl);
330
        }
gsell's avatar
gsell committed
331
    }
332
    IpplTimings::stopTimer(FunctionTimer6_m);
gsell's avatar
gsell committed
333 334
    std::ofstream timings;
    char filename[50];
335
    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
336 337 338
    if(Comm.MyPID() == 0) timings.open(filename, std::ios::app);
    double time = MPI_Wtime();

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

    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 <<
350
                                      solver_ptr->getNumIters() << "\t" <<
gsell's avatar
gsell committed
351 352 353 354 355 356
                                      //time <<" "<<
                                      minTime << "\t" <<
                                      maxTime << "\t" <<
                                      avgTime << "\t" <<
                                      numBlocks_m << "\t" <<
                                      recycleBlocks_m << "\t" <<
357
                                      nLHS_m << "\t" <<
gsell's avatar
gsell committed
358 359 360 361 362
                                      //OldLHS.size() <<"\t"<<
                                      std::endl;
    if(Comm.MyPID() == 0) timings.close();

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

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

    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
387 388 389
}


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

    int numMyGridPoints = 0;

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

    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;

412
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
413
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
414
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
415 416 417 418
                 if (bp->isInside(idx, idy, idz)) {
                    v[0] = (double)idx;
                    v[stride] = (double)idy;
                    v[stride2] = (double)idz;
gsell's avatar
gsell committed
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 446 447
                    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);
}

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

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

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

467 468
    A->PutScalar(0.0);

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

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

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

        int NumEntries = 0;

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

482
        bp->getBoundaryStencil(MyGlobalElements[i], WV, EV, SV, NV, FV, BV, CV, scaleFactor);
483
        RHS->Values()[i] *= scaleFactor;
gsell's avatar
gsell committed
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 509 510

        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;
        }

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

514 515 516 517 518
        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);
519
        } else {
gsell's avatar
gsell committed
520 521 522
            // off-diagonal entries
            A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
            // diagonal entry
523
            A->InsertGlobalValues(MyGlobalElements[i], 1, &CV, MyGlobalElements + i);
524
        }
gsell's avatar
gsell committed
525 526 527
    }

    A->FillComplete();
528

gsell's avatar
gsell committed
529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551
    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();
552 553
    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
554 555 556 557 558 559
}

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;
560
    return os;
gsell's avatar
gsell committed
561 562
}

563
#endif /* HAVE_SAAMG_SOLVER */