MGPoissonSolver.cpp 19.1 KB
Newer Older
1
#ifdef HAVE_SAAMG_SOLVER
gsell's avatar
gsell committed
2
#define DBG_STENCIL
3
#include "Algorithms/PartBunch.h"
gsell's avatar
gsell committed
4 5 6 7 8 9 10 11 12 13 14
#include "MGPoissonSolver.h"
#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"
#include "Utilities/Options.h"
#include <algorithm>
15
#include "omp.h"
gsell's avatar
gsell committed
16 17 18

using Physics::c;

19 20 21 22 23 24 25 26 27
MGPoissonSolver::MGPoissonSolver ( PartBunch &beam,
                                   Mesh_t *mesh, 
                                   FieldLayout_t *fl, 
                                   std::vector<BoundaryGeometry *> geometries, 
                                   std::string itsolver, 
                                   std::string interpl, 
                                   double tol, 
                                   int maxiters, 
                                   std::string precmode):
28 29 30 31 32 33 34
    itsBunch_m(&beam),
    mesh_m(mesh), 
    layout_m(fl), 
    geometries_m(geometries), 
    tol_m(tol), 
    maxiters_m(maxiters), 
    Comm(Ippl::getComm()) {
35

gsell's avatar
gsell committed
36 37
    domain_m = layout_m->getDomain();

38
    for (int i = 0; i < 3; i++) {
gsell's avatar
gsell committed
39 40 41 42 43 44 45 46 47 48 49 50 51
        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;
52
    else if(precmode == "NO") precmode_m = NO;
gsell's avatar
gsell committed
53 54 55 56 57 58 59

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

    hasParallelDecompositionChanged_m = true;
    useRCB_m = false;
60
    if(Ippl::Info->getOutputLevel() > 1)
gsell's avatar
gsell committed
61 62 63 64
        verbose_m = true;
    else
        verbose_m = false;

65
    // Find CURRENT geometry
gsell's avatar
gsell committed
66 67
    currentGeometry = geometries_m[0];
    if(currentGeometry->getFilename() == "") {
68
        if(currentGeometry->getTopology() == "ELLIPTIC"){
gsell's avatar
gsell committed
69
            bp = new EllipticDomain(currentGeometry->getA(), currentGeometry->getB(), orig_nr_m, hr_m, interpl);
70
	} else if (currentGeometry->getTopology() == "BOXCORNER") {
71
            bp = new BoxCornerDomain(currentGeometry->getA(), currentGeometry->getB(), currentGeometry->getC(), currentGeometry->getLength(),currentGeometry->getL1(), currentGeometry->getL2(), orig_nr_m, hr_m, interpl);
72
            bp->Compute(itsBunch_m->get_hr());
gsell's avatar
gsell committed
73 74 75 76
        } else {
            ERRORMSG("Geometry not known" << endl);
            exit(1);
        }
77
    } else 
78 79
	bp = new ArbitraryDomain(currentGeometry, orig_nr_m, hr_m, interpl);

gsell's avatar
gsell committed
80 81 82 83
    Map = 0;
    A = Teuchos::null;
    LHS = Teuchos::null;
    RHS = Teuchos::null;
84 85
    MLPrec = 0;
    
gsell's avatar
gsell committed
86 87
    numBlocks_m = Options::numBlocks;
    recycleBlocks_m = Options::recycleBlocks;
88
    nLHS_m = Options::nLHS;
gsell's avatar
gsell committed
89 90
    SetupMLList();
    SetupBelosList();
91 92
    problem_ptr = rcp(new Belos::LinearProblem<ST, MV, OP>);
    // setup Belos solver 
gsell's avatar
gsell committed
93
    if(numBlocks_m == 0 || recycleBlocks_m == 0)
94
        solver_ptr = rcp(new Belos::BlockCGSolMgr<double, MV, OP>());
gsell's avatar
gsell committed
95
    else
96 97
        solver_ptr = rcp(new Belos::RCGSolMgr<double, MV, OP>());
    solver_ptr->setParameters(rcp(&belosList, false));
gsell's avatar
gsell committed
98 99 100 101
    convStatusTest = rcp(new Belos::StatusTestGenResNorm<ST, MV, OP> (tol));
    convStatusTest->defineScaleForm(Belos::NormOfRHS, Belos::TwoNorm);

    //all timers used here
102 103 104 105 106 107 108 109
    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
110 111 112 113 114 115 116 117 118 119
}


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() {
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
    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;
    prec_m = Teuchos::null;
gsell's avatar
gsell committed
135 136 137 138 139 140
}

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

141
void MGPoissonSolver::computeMap(NDIndex<3> localId) {
142 143
    if (hasParallelDecompositionChanged_m){
        if (Map !=0 ) delete Map;
gsell's avatar
gsell committed
144
        if(useRCB_m)
145
            redistributeWithRCB(localId);
146 147
        else if ( bp->getType() == "Geometric" )
            IPPLToMap3DGeo(localId);    	    
148
        else if ( bp->getType() == "Elliptic" )
149
            IPPLToMap3D(localId);
150 151

        extrapolateLHS();
152 153
    }
}
gsell's avatar
gsell committed
154

155 156 157 158 159
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.
160 161 162 163 164 165 166 167 168 169 170
    //we also have to redistribute LHS
    
    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
171

172
    //...and all previously saved LHS
173
    std::deque< Epetra_Vector >::iterator it = OldLHS.begin();
174 175 176 177 178 179 180 181 182 183
    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
184

185 186 187
    // extrapolate last OldLHS.size LHS to get a new start vector
    it = OldLHS.begin();
    if(nLHS_m == 0 || OldLHS.size()==0)
188 189 190
        LHS->PutScalar(1.0);
    else if(OldLHS.size() == 1)
        *LHS = *it;
191 192 193
    else if(OldLHS.size() == 2)
        LHS->Update(2.0, *it++, -1.0, *it, 0.0);
    else if(OldLHS.size() > 2) {
194
        int n = OldLHS.size();
195 196 197
        P = rcp(new Epetra_MultiVector(*Map, nLHS_m, false));
        for(int i = 0; i < n; ++i) {
           *(*P)(i) = *it++;
198
        }
199 200 201 202
        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);
           }
203 204
        }
        *LHS = *(*P)(0);
205 206
     } else
        throw OpalException("MGPoissonSolver", "Invalid number of old LHS: " + OldLHS.size());
gsell's avatar
gsell committed
207 208
}

209

gsell's avatar
gsell committed
210 211 212 213 214
// 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) {
215

gsell's avatar
gsell committed
216 217 218
    nr_m[0] = orig_nr_m[0];
    nr_m[1] = orig_nr_m[1];
    nr_m[2] = orig_nr_m[2];
219

220
    //bp->setMinMaxZ(itsBunch_m->get_origin()[2], itsBunch_m->get_maxExtend()[2]);
gsell's avatar
gsell committed
221
    bp->setNr(nr_m);
222 223
    NDIndex<3> localId = layout_m->getLocalNDIndex();

224
    IpplTimings::startTimer(FunctionTimer1_m);
225
    if ( bp->getType() == "Geometric" ) {
226
        bp->Compute(hr, localId, itsBunch_m->getGlobalMeanR(), itsBunch_m->getGlobalToLocalQuaternion());
227
    } else if ( bp->getType() == "Elliptic"  )
228
        bp->Compute(hr);
229
    IpplTimings::stopTimer(FunctionTimer1_m);
230 231

    // Define the Map
232
    IpplTimings::startTimer(FunctionTimer2_m);
233
    computeMap(localId);	
234
    IpplTimings::stopTimer(FunctionTimer2_m);
gsell's avatar
gsell committed
235

236
    // Allocate the RHS with the new Epetra Map
237 238
    if (Teuchos::is_null(RHS)) 
      RHS = rcp(new Epetra_Vector(*Map));
239
    RHS->PutScalar(0.0);
240 241 242
    // get charge densities from IPPL field and store in Epetra vector (RHS)
    IpplTimings::startTimer(FunctionTimer3_m);
    int id = 0;
243
    float scaleFactor = itsBunch_m->getdT(); 
244
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
245
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
246
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
247
                  if (bp->isInside(idx, idy, idz))
248
                    RHS->Values()[id++] = 4*M_PI*rho[idx][idy][idz].get()/scaleFactor;
gsell's avatar
gsell committed
249 250 251
            }
        }
    }
252
    IpplTimings::stopTimer(FunctionTimer3_m);
253

gsell's avatar
gsell committed
254
    // build discretization matrix
255
    IpplTimings::startTimer(FunctionTimer4_m);
256 257
    if (Teuchos::is_null(A))
        A = rcp(new Epetra_CrsMatrix(Copy, *Map,  7, true));
258
    ComputeStencil(hr, RHS);
259
    IpplTimings::stopTimer(FunctionTimer4_m);
gsell's avatar
gsell committed
260 261 262
#ifdef DBG_STENCIL
    EpetraExt::RowMatrixToMatlabFile("DiscrStencil.dat", *A);
#endif
263

264
    IpplTimings::startTimer(FunctionTimer5_m);
265 266 267 268 269 270 271 272
    if ( MLPrec == 0 ){
        MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList_m);
    } else if (precmode_m == REUSE_HIERARCHY) {
        MLPrec->ReComputePreconditioner();
    } else if (precmode_m == REUSE_PREC){
    } else if (precmode_m == STD_PREC) {
        delete MLPrec; MLPrec = 0;
        MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList_m);
273
    }
274
    IpplTimings::stopTimer(FunctionTimer5_m);
275

gsell's avatar
gsell committed
276 277
    // setup preconditioned iterative solver
    // use old LHS solution as initial guess
278
    IpplTimings::startTimer(FunctionTimer6_m);
279 280 281 282 283 284 285 286 287

    problem_ptr->setOperator(A);
    problem_ptr->setLHS(LHS);
    problem_ptr->setRHS(RHS);
    prec_m = Teuchos::rcp ( new Belos::EpetraPrecOp ( rcp(MLPrec,false)));
    problem_ptr->setLeftPrec(prec_m);
    solver_ptr->setProblem( problem_ptr);
    if (!problem_ptr->isProblemSet()) {
        if (problem_ptr->setProblem() == false) {
288
            ERRORMSG("Belos::LinearProblem failed to set up correctly!" << endl);
289
        }
gsell's avatar
gsell committed
290
    }
291
    IpplTimings::stopTimer(FunctionTimer6_m);
gsell's avatar
gsell committed
292 293
    std::ofstream timings;
    char filename[50];
294
    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
295 296 297
    if(Comm.MyPID() == 0) timings.open(filename, std::ios::app);
    double time = MPI_Wtime();

298
    IpplTimings::startTimer(FunctionTimer7_m);
299
    	solver_ptr->solve();
300
    IpplTimings::stopTimer(FunctionTimer7_m);
gsell's avatar
gsell committed
301 302 303 304 305 306 307 308

    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 <<
309
                                      solver_ptr->getNumIters() << "\t" <<
gsell's avatar
gsell committed
310 311 312 313 314 315
                                      //time <<" "<<
                                      minTime << "\t" <<
                                      maxTime << "\t" <<
                                      avgTime << "\t" <<
                                      numBlocks_m << "\t" <<
                                      recycleBlocks_m << "\t" <<
316
                                      nLHS_m << "\t" <<
gsell's avatar
gsell committed
317 318 319 320 321
                                      //OldLHS.size() <<"\t"<<
                                      std::endl;
    if(Comm.MyPID() == 0) timings.close();

    // Store new LHS in OldLHS
322 323
    if(nLHS_m > 1) OldLHS.push_front(*(LHS.get()));
    if(OldLHS.size() > nLHS_m) OldLHS.pop_back();
gsell's avatar
gsell committed
324 325

    //now transfer solution back to IPPL grid
326 327
    IpplTimings::startTimer(FunctionTimer8_m);
    id = 0;
gsell's avatar
gsell committed
328
    rho = 0.0;
329
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
330
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
331
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
332 333
                  NDIndex<3> l(Index(idx, idx), Index(idy, idy), Index(idz, idz));
                  if (bp->isInside(idx, idy, idz))
334
                     rho.localElement(l) = LHS->Values()[id++]*scaleFactor;
gsell's avatar
gsell committed
335 336 337
            }
        }
    }
338
    IpplTimings::stopTimer(FunctionTimer8_m);
339 340
    if ( hasParallelDecompositionChanged_m )
        deletePtr();
gsell's avatar
gsell committed
341 342 343
}


344
void MGPoissonSolver::redistributeWithRCB(NDIndex<3> localId) {
gsell's avatar
gsell committed
345 346 347

    int numMyGridPoints = 0;

348
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
349
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
350
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
351
                 if (bp->isInside(idx, idy, idz))
gsell's avatar
gsell committed
352
                    numMyGridPoints++;
353
             }
gsell's avatar
gsell committed
354
        }
355
     }
gsell's avatar
gsell committed
356 357 358 359 360 361 362 363 364 365

    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;

366
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
367
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
368
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
369 370 371 372
                 if (bp->isInside(idx, idy, idz)) {
                    v[0] = (double)idx;
                    v[stride] = (double)idy;
                    v[stride2] = (double)idz;
gsell's avatar
gsell committed
373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401
                    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);
}

402
void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
gsell's avatar
gsell committed
403 404 405
    int NumMyElements = 0;
    vector<int> MyGlobalElements;

406
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
407
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
408
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
409
                if (bp->isInside(idx, idy, idz)) {
410
                    MyGlobalElements.push_back(bp->getIdx(idx, idy, idz));
gsell's avatar
gsell committed
411 412 413
                    NumMyElements++;
                }
            }
414 415
	}
    }
gsell's avatar
gsell committed
416 417 418
    Map = new Epetra_Map(-1, NumMyElements, &MyGlobalElements[0], 0, Comm);
}

419 420 421
void MGPoissonSolver::IPPLToMap3DGeo(NDIndex<3> localId) {
    int NumMyElements = 0;
    vector<int> MyGlobalElements;
422 423 424
/*   std::cout << Ippl::myNode() << " In IPPLToMap3DGeo localId: " << localId[0].first() << " " << localId[0].last()  << "\n" 
					<< localId[1].first() << " " << localId[1].last()  << "\n" 
					<< localId[2].first() << " " << localId[2].last()  << std::endl; */
425 426 427 428 429 430 431 432 433 434 435
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
                    MyGlobalElements.push_back(bp->getIdx(idx, idy, idz));
                    NumMyElements++;
            }
	}
    }
    Map = new Epetra_Map(-1, NumMyElements, &MyGlobalElements[0], 0, Comm);
}

gsell's avatar
gsell committed
436 437
void MGPoissonSolver::ComputeStencil(Vector_t hr, Teuchos::RCP<Epetra_Vector> RHS) {

438 439
    A->PutScalar(0.0);

gsell's avatar
gsell committed
440 441 442 443 444 445
    int NumMyElements = Map->NumMyElements();
    int *MyGlobalElements = Map->MyGlobalElements();

    vector<double> Values(6);
    vector<int> Indices(6);

446
    for(int i = 0 ; i < NumMyElements ; i++) {
gsell's avatar
gsell committed
447 448 449

        int NumEntries = 0;

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

453
        bp->getBoundaryStencil(MyGlobalElements[i], WV, EV, SV, NV, FV, BV, CV, scaleFactor);
454
//        RHS->Values()[i] *= scaleFactor; 
gsell's avatar
gsell committed
455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481

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

482
    	// if matrix has already been filled (FillComplete()) we can only
gsell's avatar
gsell committed
483
        // replace entries
484 485 486 487 488 489 490
        
        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);
        } else { 
gsell's avatar
gsell committed
491 492 493
            // off-diagonal entries
            A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
            // diagonal entry
494 495
            A->InsertGlobalValues(MyGlobalElements[i], 1, &CV, MyGlobalElements + i);
        } 
gsell's avatar
gsell committed
496 497 498
    }

    A->FillComplete();
499

gsell's avatar
gsell committed
500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522
    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();
523 524
    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
525 526 527 528 529 530
}

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;
531
    return os;	
gsell's avatar
gsell committed
532 533
}

534
#endif /* HAVE_SAAMG_SOLVER */