PetscSolverFeti.cc 55.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
#include "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
#include "MatrixVector.h"
15
#include "parallel/PetscSolverFeti.h"
16
#include "parallel/PetscSolverFetiStructs.h"
17
18
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
19
#include "parallel/PetscSolverGlobalMatrix.h"
20
#include "io/VtkWriter.h"
21
22
23
24
25

namespace AMDiS {

  using namespace std;

26
27
28
29
30
31
  double FetiTimings::fetiSolve = 0.0;
  double FetiTimings::fetiSolve01 = 0.0;
  double FetiTimings::fetiSolve02 = 0.0;
  double FetiTimings::fetiPreconditioner = 0.0;


32
33
34
35
36
37
38
  // y = mat * x
  int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y)
  {
    // S_PiPi = K_PiPi - K_PiB inv(K_BB) K_BPi

    void *ctx;
    MatShellGetContext(mat, &ctx);
39
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
40

41
    MatMult(data->subSolver->getMatInteriorCoarse(), x, data->tmp_vec_b);
42
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
43
44
45
    MatMult(data->subSolver->getMatCoarseInterior(), data->tmp_vec_b, 
	    data->tmp_vec_primal);
    MatMult(data->subSolver->getMatCoarse(), x, y);
46
47
48
49
50
51
52
53
54
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
55
56
    FUNCNAME("petscMultMatFeti()");

57
58
    //    F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)
    // => F = L [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(L)
59

60
61
    double wtime = MPI::Wtime();

62
63
    void *ctx;
    MatShellGetContext(mat, &ctx);
64
    FetiData* data = static_cast<FetiData*>(ctx);
65

Thomas Witkowski's avatar
Thomas Witkowski committed
66
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
67
68

    double wtime01 = MPI::Wtime();
69
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
70
71
72

    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

Thomas Witkowski's avatar
Thomas Witkowski committed
73
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
74

75
76
    MatMult(data->subSolver->getMatCoarseInterior(), 
	    data->tmp_vec_b, data->tmp_vec_primal);
77
78

    wtime01 = MPI::Wtime();
79
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
80
81
    FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);

82
83
    MatMult(data->subSolver->getMatInteriorCoarse(), 
	    data->tmp_vec_primal, data->tmp_vec_b);
84
85

    wtime01 = MPI::Wtime();
86
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
87
88
    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

Thomas Witkowski's avatar
Thomas Witkowski committed
89
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
90

Thomas Witkowski's avatar
Thomas Witkowski committed
91
    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
92

93
94
    FetiTimings::fetiSolve += (MPI::Wtime() - wtime);

95
96
97
98
    return 0;
  }


99
  // y = PC * x
100
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
101
  {
102
103
    double wtime = MPI::Wtime();

104
    // Get data for the preconditioner
105
106
    void *ctx;
    PCShellGetContext(pc, &ctx);
107
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
108

109
    // Multiply with scaled Lagrange constraint matrix.
110
111
112
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


113
    // === Restriction of the B nodes to the boundary nodes. ===
114

115
116
117
118
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
119

120
121
122
123
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

124
125
126
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
127
128

    VecRestoreArray(data->tmp_vec_b, &local_b);
129
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
130
131


132
    // === K_DD - K_DI inv(K_II) K_ID ===
133

134
    MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);
135

136
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
137
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
138
139
140
141
142
143
144
145
146
147
    MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0);

    VecAXPBY(data->tmp_vec_duals0, 1.0, -1.0, data->tmp_vec_duals1);


    // === Prolongation from local dual nodes to B nodes.

    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

148
149
150
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
151
152
153
154
155
156
157
158

    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // Multiply with scaled Lagrange constraint matrix.
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

159
160
    FetiTimings::fetiPreconditioner += (MPI::Wtime() - wtime);

161
162
163
164
165
166
167
168
169
170
171
172
173
174
    return 0;
  }


  // y = PC * x
  PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y)
  {
    // Get data for the preconditioner
    void *ctx;
    PCShellGetContext(pc, &ctx);
    FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);

    // Multiply with scaled Lagrange constraint matrix.
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
175
176


177
    // === Restriction of the B nodes to the boundary nodes. ===
178

179
180
181
182
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
183

184
    PetscScalar *local_b, *local_duals;
185
    VecGetArray(data->tmp_vec_b, &local_b);
186
    VecGetArray(data->tmp_vec_duals0, &local_duals);
187

Thomas Witkowski's avatar
Thomas Witkowski committed
188
189
190
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
191
192

    VecRestoreArray(data->tmp_vec_b, &local_b);
193
194
195
196
197
198
199
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

    MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);

200

201
    // === Prolongation from local dual nodes to B nodes.
202

203
204
205
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

Thomas Witkowski's avatar
Thomas Witkowski committed
206
207
208
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
209

210
211
    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
212

213
214

    // Multiply with scaled Lagrange constraint matrix.
215
216
217
218
219
220
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


221
222
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
223
      schurPrimalSolver(0),
224
      multiLevelTest(false),
Thomas Witkowski's avatar
Thomas Witkowski committed
225
      subdomain(NULL),
226
227
      meshLevel(0),
      rStartInterior(0),
Thomas Witkowski's avatar
Thomas Witkowski committed
228
229
      nGlobalOverallInterior(0),
      printTimings(false)
230
231
232
233
234
235
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
236
      MSG("Create FETI-DP solver with no preconditioner!\n");
237
238
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
239
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
240
241
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
242
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
243
244
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
245
246
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
247
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
248
249
250
251
252

    Parameters::get("parallel->feti->schur primal solver", schurPrimalSolver);
    TEST_EXIT(schurPrimalSolver == 0 || schurPrimalSolver == 1)
      ("Wrong solver \"%d\"for the Schur primal complement!\n", 
       schurPrimalSolver);
253
254

    Parameters::get("parallel->multi level test", multiLevelTest);
255
256
    if (multiLevelTest)
      meshLevel = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
257
258

    Parameters::get("parallel->print timings", printTimings);
259
260
261
  }


262
  void PetscSolverFeti::initialize(vector<const FiniteElemSpace*> feSpaces)
263
  {
264
265
    FUNCNAME("PetscSolverFeti::initialize()");

266
267
268
    TEST_EXIT_DBG(meshLevel + 1 == meshDistributor->getMeshLevelData().getLevelNumber())
      ("Mesh hierarchy does not contain %d levels!\n", meshLevel + 1);

269
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();
270
    vector<const FiniteElemSpace*>& uniqueFe = meshDistributor->getFeSpaces();
271

Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
    if (subdomain == NULL) {
      subdomain = new PetscSolverGlobalMatrix();
274

275
      if (meshLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
276
	subdomain->setMeshDistributor(meshDistributor, 
277
				      mpiCommGlobal, mpiCommLocal);
278
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
279
	subdomain->setMeshDistributor(meshDistributor, 
280
281
				      levelData.getMpiComm(meshLevel - 1),
				      levelData.getMpiComm(meshLevel));
Thomas Witkowski's avatar
Thomas Witkowski committed
282
	subdomain->setLevel(meshLevel);
283
284
      }
    }
285

286
287
288
289
    primalDofMap.init(levelData, feSpaces, uniqueFe);   
    dualDofMap.init(levelData, feSpaces, uniqueFe, false);
    localDofMap.init(levelData, feSpaces, uniqueFe, meshLevel != 0);
    lagrangeMap.init(levelData, feSpaces, uniqueFe);
Thomas Witkowski's avatar
Thomas Witkowski committed
290

291
292
    if (fullInterface != NULL)
      interfaceDofMap.init(levelData, fullInterface);
Thomas Witkowski's avatar
Thomas Witkowski committed
293

294
    if (fetiPreconditioner != FETI_NONE) {
295
296
297
      TEST_EXIT(meshLevel == 0)
	("Dirichlet preconditioner not yet implemented for multilevel FETI-DP\n");

298
      interiorDofMap.init(levelData, feSpaces, uniqueFe, false);
299
    }
300
301
302
303
304
305
306
  }


  void PetscSolverFeti::createFetiData()
  {
    FUNCNAME("PetscSolverFeti::createFetiData()");

Thomas Witkowski's avatar
Thomas Witkowski committed
307
308
    double timeCounter = MPI::Wtime();

309
310
311
    TEST_EXIT(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT(meshDistributor->getFeSpaces().size() > 0)
      ("No FE space defined in mesh distributor!\n");
312

313
314
    MeshLevelData& levelData = meshDistributor->getMeshLevelData();

315
316
317
318
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
319
    if (fetiPreconditioner != FETI_NONE)
320
321
      interiorDofMap.clear();

322
323
    primalDofMap.setDofComm(meshDistributor->getDofComm());
    lagrangeMap.setDofComm(meshDistributor->getDofComm());
324

325
326
    primalDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    dualDofMap.setMpiComm(levelData.getMpiComm(0), 0);
327
    lagrangeMap.setMpiComm(levelData.getMpiComm(0), 0);
328
    localDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
329
330
    if (fetiPreconditioner != FETI_NONE)
      interiorDofMap.setMpiComm(levelData.getMpiComm(meshLevel), meshLevel);
331

332
333
334
    if (meshLevel == 0)
      localDofMap.setDofComm(meshDistributor->getDofComm());
    else
335
336
      localDofMap.setDofComm(meshDistributor->getDofCommSd());

337
338
339
340
341
342
    if (fullInterface != NULL) {
      interfaceDofMap.clear();
      interfaceDofMap.setDofComm(meshDistributor->getDofComm());
      interfaceDofMap.setMpiComm(levelData.getMpiComm(0), 0);
    }

343
344
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
345
    
346
347
      createPrimals(feSpace);  

348
349
350
      createDuals(feSpace);

      createInterfaceNodes(feSpace);
351

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
352
      createIndexB(feSpace);     
353
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
354
355
356
357

    primalDofMap.update();
    dualDofMap.update();
    localDofMap.update();
358
    if (fetiPreconditioner != FETI_NONE)
359
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
360

361
362
363
    if (fullInterface != NULL)
      interfaceDofMap.update();

364
365
366
367
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createLagrange(feSpace);
    }
368

369
370
371
    lagrangeMap.update();


372
373
374
375
376
377
378
379
380
381
382
383
    // === ===

    if (meshLevel == 0) {
      rStartInterior = 0;
      nGlobalOverallInterior = localDofMap.getOverallDofs();
    } else {
      MeshLevelData& levelData = meshDistributor->getMeshLevelData();

      int groupRowsInterior = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	groupRowsInterior = localDofMap.getOverallDofs();

384
      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
385
386
387
388
389
390
391
392
393
			   rStartInterior, nGlobalOverallInterior);

      int tmp = 0;
      if (levelData.getMpiComm(1).Get_rank() == 0)
	tmp = rStartInterior;

      levelData.getMpiComm(1).Allreduce(&tmp, &rStartInterior, 1, MPI_INT, MPI_SUM);
    }

394
    MSG("FETI-DP data created on mesh level %d\n", meshLevel);
Thomas Witkowski's avatar
Thomas Witkowski committed
395
396
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
397
      MSG("FETI-DP data for %d-ith FE space:\n", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
398

399
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
400
401
	  primalDofMap[feSpace].nRankDofs, 
	  primalDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
402
      
403
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
404
405
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
406

407
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
408
409
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
410

411
      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
412
413
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
414
    }
415

Thomas Witkowski's avatar
Thomas Witkowski committed
416
417
    subdomain->setDofMapping(&localDofMap);
    subdomain->setCoarseSpaceDofMapping(&primalDofMap);
Thomas Witkowski's avatar
Thomas Witkowski committed
418
419

    if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
420
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
421
      timeCounter = MPI::Wtime() - timeCounter;
422
      MSG("FETI-DP timing 01: %.5f seconds (creation of basic data structures)\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
423
424
	  timeCounter);
    }
425
426
427
  }


428
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
429
  {
430
    FUNCNAME("PetscSolverFeti::createPrimals()");  
431

432
433
434
    if (feSpace == fullInterface)
      return;

435
436
437
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

Thomas Witkowski's avatar
Thomas Witkowski committed
438
    // Set of DOF indices that are considered to be primal variables.
439
    DofContainerSet& vertices = 
440
      meshDistributor->getBoundaryDofInfo(feSpace, meshLevel).geoDofs[VERTEX];
441
442

    DofIndexSet primals;
443
    for (DofContainerSet::iterator it = vertices.begin(); 
444
445
	 it != vertices.end(); ++it) {
      double e = 1e-8;
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
446
      if (meshLevel == 0) {
447
	primals.insert(**it);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
448
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
449
450
451
	WorldVector<double> c;
	feSpace->getMesh()->getDofIndexCoords(*it, feSpace, c);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
452
453
454
455
456
457
458
	if (fabs(c[0]) < e || fabs(c[1]) < e ||
	    fabs(c[0] - 25.0) < e || fabs(c[1] - 25.0) < e ||
	    (fabs(c[0] - 12.5) < e && fabs(c[1] - 12.5) < e)) {
	  MSG("PRIMAL COORD %f %f\n", c[0], c[1]);
	  primals.insert(**it);
	}
      }
459
    }
460

461
462
463
464

    // === Calculate the number of primals that are owned by the rank and ===
    // === create local indices of the primals starting at zero.          ===

465
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
466
467
      if (meshDistributor->getDofMap()[feSpace].isRankDof(*it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
468
      else
469
  	primalDofMap[feSpace].insertNonRankDof(*it);
470
471
472
  }


473
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
474
475
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
476

477
478
479
    if (feSpace == fullInterface)
      return;

480
481
482
    // === Create global index of the dual nodes on each rank. ===

    DofContainer allBoundaryDofs;
483
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);
484
485
486
487

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
488
489
	if (meshLevel == 0) {
	  dualDofMap[feSpace].insertRankDof(**it);
490
491
492
493
 	} else {
 	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(**it))
 	    dualDofMap[feSpace].insertRankDof(**it);
 	}	  
494
495
496
  }

  
497
498
499
500
501
502
503
504
505
506
507
508
  void PetscSolverFeti::createInterfaceNodes(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createInterfaceNodes()");

    if (feSpace != fullInterface)
      return;

    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(feSpace, meshLevel, allBoundaryDofs);

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
509
510
      if (meshDistributor->getDofMap()[feSpace].isRankDof(**it))
	interfaceDofMap[feSpace].insertRankDof(**it);
511
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
512
	interfaceDofMap[feSpace].insertNonRankDof(**it);
513
514
515
  }


516
517
518
519
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

520
521
522
    if (feSpace == fullInterface)
      return;

523
524
    boundaryDofRanks[feSpace].clear();

Thomas Witkowski's avatar
Thomas Witkowski committed
525
526
527
    // Stores for all rank owned communication DOFs, if the counterpart is
    // a rank owned DOF in its subdomain. Thus, the following map stores to
    // each rank number all DOFs that fulfill this requirenment.
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
528
    map<int, std::set<DegreeOfFreedom> > sdRankDofs;
529

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
530
    if (meshLevel > 0) {
531
      StdMpi<vector<int> > stdMpi(mpiCommGlobal);
532
533
534
535
536
537
538
539
540

      for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank()) {

	vector<int> subdomainRankDofs;
	subdomainRankDofs.reserve(it.getDofs().size());

	for (; !it.endDofIter(); it.nextDof()) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
541
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))
542
	    subdomainRankDofs.push_back(1);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
543
	  else
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
	    subdomainRankDofs.push_back(0);
	}

	stdMpi.send(it.getRank(), subdomainRankDofs);
      }	     

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace);
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
				meshLevel, feSpace); 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
559
560
561
562
563
564
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof())
	  if (!isPrimal(feSpace, it.getDofIndex()))
	    if (stdMpi.getRecvData(it.getRank())[it.getDofCounter()] == 1)
	      sdRankDofs[it.getRank()].insert(it.getDofIndex());
    }
565

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
566
567
568
569
570
571
572
    if (dualDofMap[feSpace].nLocalDofs == 0)
      return;


    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===

573
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
574
575
576
577
578
579
580
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), 
			      meshLevel, feSpace); 
	 !it.end(); it.nextRank()) {
      for (; !it.endDofIter(); it.nextDof()) {
	if (!isPrimal(feSpace, it.getDofIndex())) {
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);

581
 	  if (meshLevel == 0 ||
Thomas Witkowski's avatar
Thomas Witkowski committed
582
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
583
	    boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());	  
584
585
586
587
	}
      }
    }

588
589
590
591

    // === Communicate these sets for all rank owned dual nodes to other ===
    // === ranks that also have this node.                               ===

592
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
593

594
    for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), meshLevel, feSpace);
595
596
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
597
	if (!isPrimal(feSpace, it.getDofIndex()))
598
599
600
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && sdRankDofs[it.getRank()].count(it.getDofIndex())))
	    stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
601
602
603

    stdMpi.updateSendDataSize();

604
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
605
	 !it.end(); it.nextRank()) {
606
      bool recvFromRank = false;
607
      for (; !it.endDofIter(); it.nextDof()) {
608
	if (!isPrimal(feSpace, it.getDofIndex())) {
609
610
611
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex()))) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
612
613
614
	    recvFromRank = true;
	    break;
	  }
615
	}
616
      }
617
618

      if (recvFromRank)
619
	stdMpi.recv(it.getRank());
620
    }
621

622
623
    stdMpi.startCommunication();

624
    for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), meshLevel, feSpace); 
625
	 !it.end(); it.nextRank()) {
626
      int i = 0;
627
      for (; !it.endDofIter(); it.nextDof())
628
	if (!isPrimal(feSpace, it.getDofIndex()))
629
630
631
 	  if (meshLevel == 0 ||
 	      (meshLevel > 0 && 
 	       meshDistributor->getDofMapSd()[feSpace].isRankDof(it.getDofIndex())))	    
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
632
633
	    boundaryDofRanks[feSpace][it.getDofIndex()] = 
	      stdMpi.getRecvData(it.getRank())[i++];
634
635
	  else
	    lagrangeMap[feSpace].insertNonRankDof(it.getDofIndex());
636
637
    }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
638

639
640
641
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

642
    int nRankLagrange = 0;
643
    DofMap& dualMap = dualDofMap[feSpace].getMap();
644
    for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
645
646
      if (meshDistributor->getDofMap()[feSpace].isRankDof(it->first)) {
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
647
	int degree = boundaryDofRanks[feSpace][it->first].size();
648
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
649
      } else {
650
	lagrangeMap[feSpace].insertNonRankDof(it->first);
651
652
      }
    }
653
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
654
655
656
  }


657
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
658
  {
659
    FUNCNAME("PetscSolverFeti::createIndexB()");
660

661
    DOFAdmin* admin = feSpace->getAdmin();
662
663
664
665

    // === To ensure that all interior node on each rank are listen first in ===
    // === the global index of all B nodes, insert all interior nodes first, ===
    // === without defining a correct index.                                 ===
666

667
    int nLocalInterior = 0;    
668
    for (int i = 0; i < admin->getUsedSize(); i++) {
669
      if (admin->isDofFree(i) == false && 
670
	  isPrimal(feSpace, i) == false &&
671
672
	  isDual(feSpace, i) == false &&
	  isInterface(feSpace, i) == false) {
673
674
	if (meshLevel == 0) {
	  localDofMap[feSpace].insertRankDof(i, nLocalInterior);
675

676
677
678
679
680
681
682
683
684
	  if (fetiPreconditioner != FETI_NONE)
	    interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);
	  
	  nLocalInterior++;	  
	} else {
	  if (meshDistributor->getDofMapSd()[feSpace].isRankDof(i))
	    localDofMap[feSpace].insertRankDof(i);
	  else
	    localDofMap[feSpace].insertNonRankDof(i);
685
686
687

	  TEST_EXIT_DBG(fetiPreconditioner == FETI_NONE)
	    ("Not yet implemnted!\n");
688
	}
689
      }
690
    }
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
691
    
692
693
    // === And finally, add the global indicies of all dual nodes. ===

694
695
696
697
698
699
700
701
702
703
    for (DofMap::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
      if (meshLevel == 0)
	localDofMap[feSpace].insertRankDof(it->first);
      else {
	if (meshDistributor->getDofMapSd()[feSpace].isRankDof(it->first))
	  localDofMap[feSpace].insertRankDof(it->first);
	else 
	  localDofMap[feSpace].insertNonRankDof(it->first);
      }
704
705
706
  }


707
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
708
709
710
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

Thomas Witkowski's avatar
Thomas Witkowski committed
711
    double wtime = MPI::Wtime();
712
    int mpiRank = meshDistributor->getMpiRank();
Thomas Witkowski's avatar
Thomas Witkowski committed
713

714
715
    // === Create distributed matrix for Lagrange constraints. ===

716
717
718
719
720
    MatCreateAIJ(mpiCommGlobal,
		 lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		 lagrangeMap.getOverallDofs(), nGlobalOverallInterior,
		 2, PETSC_NULL, 2, PETSC_NULL,
		 &mat_lagrange);
721
    MatSetOption(mat_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
722

723
724
725
726
727
728
729
    // === Create for all duals the corresponding Lagrange constraints. On ===
    // === each rank we traverse all pairs (n, m) of ranks, with n < m,    ===
    // === that contain this node. If the current rank number is r, and    ===
    // === n == r, the rank sets 1.0 for the corresponding constraint, if  ===
    // === m == r, than the rank sets -1.0 for the corresponding           ===
    // === constraint.                                                     ===

Thomas Witkowski's avatar
Thomas Witkowski committed
730
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
731
      DofMap &dualMap = dualDofMap[feSpaces[k]].getMap();
732

733
      for (DofMap::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
734
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
735
736
737
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
738
	int index = lagrangeMap.getMatIndex(k, it->first);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
739

Thomas Witkowski's avatar
Thomas Witkowski committed
740
	// Copy set of all ranks that contain this dual node.
741
742
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
743
744
	// Number of ranks that contain this dual node.
	int degree = W.size();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
745
746

	TEST_EXIT_DBG(degree > 1)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
749
750
	
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
751
752
	      int colIndex = 
		localDofMap.getMatIndex(k, it->first) + rStartInterior;
753

754
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
755
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
756
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
757
	    index++;	      
758
759
760
761
762
763
764
	  }
	}
      }
    }

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
765

Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
766
767
    if (printTimings) {
      MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
768
769
      MSG("FETI-DP timing 05: %.5f seconds (creation of lagrange constraint matrix)\n", 
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
770
    }
771
772
773
  }


774
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
775
  {
776
    FUNCNAME("PetscSolverFeti::createSchurPrimalKsp()");
777

Thomas Witkowski's avatar
Thomas Witkowski committed
778
    if (schurPrimalSolver == 0) {
779
780
      MSG("Create iterative schur primal solver!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
781
      schurPrimalData.subSolver = subdomain;
Thomas Witkowski's avatar
Thomas Witkowski committed
782
      
783
      VecCreateMPI(mpiCommGlobal, 
784
		   localDofMap.getRankDofs(), 
785
		   nGlobalOverallInterior,
Thomas Witkowski's avatar
Thomas Witkowski committed
786
		   &(schurPrimalData.tmp_vec_b));
787
      VecCreateMPI(mpiCommGlobal, 
788
789
		   primalDofMap.getRankDofs(), 
		   primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
790
791
		   &(schurPrimalData.tmp_vec_primal));

792
      MatCreateShell(mpiCommGlobal,
793
794
795
796
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), 
		     primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
797
798
799
800
801
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
802
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
803
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
804
805
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
806
807
      KSPSetFromOptions(ksp_schur_primal);
    } else {
808
809
      MSG("Create direct schur primal solver!\n");

810
811
      double wtime = MPI::Wtime();

812
813
814
      TEST_EXIT_DBG(meshLevel == 0)
	("Does not support for multilevel, check usage of localDofMap.\n");

815
816
817
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
      int nRowsRankB = localDofMap.getRankDofs();
818

Thomas Witkowski's avatar
Thomas Witkowski committed
819
      Mat matBPi;
820
821
822
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankB, nRowsRankPrimal, 
		   nGlobalOverallInterior, nRowsOverallPrimal,
823
824
825
		   150, PETSC_NULL, 150, PETSC_NULL, &matBPi);
      MatSetOption(matBPi, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

826

827
828
829
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
830

831
832
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

833
      for (int i = 0; i < nRowsRankB; i++) {
834
	PetscInt row = localDofMap.getStartDofs() + i;
835
	MatGetRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
836
837
838
839
840

	for (int j = 0; j < nCols; j++)
	  if (values[j] != 0.0)
	    mat_b_primal_cols[cols[j]].push_back(make_pair(i, values[j]));

841
	MatRestoreRow(subdomain->getMatInteriorCoarse(), row, &nCols, &cols, &values);
842
843
      }

844
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
845
		primalDofMap.getLocalDofs())
846
	("Should not happen!\n");
847

848
849
850
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
851
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
852
853
854
855
856
857
858
859

 	for (unsigned int i = 0; i < it->second.size(); i++) 
 	  VecSetValue(tmpVec, 
 		      it->second[i].first, it->second[i].second, INSERT_VALUES);

	VecAssemblyBegin(tmpVec);
	VecAssemblyEnd(tmpVec);

Thomas Witkowski's avatar
Thomas Witkowski committed
860
	subdomain->solve(tmpVec, tmpVec);
861

Thomas Witkowski's avatar
Thomas Witkowski committed
862
	PetscScalar *tmpValues;
863
	VecGetArray(tmpVec, &tmpValues);
864
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
865
	  MatSetValue(matBPi, 
866
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
867
868
869
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
870
871
872
873
874
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
875
876
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
877

878
879
      MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES, 
		   &mat_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
880
881

      Mat matPrimal;
882
883
      MatMatMult(subdomain->getMatCoarseInterior(), matBPi, MAT_INITIAL_MATRIX, 
		 PETSC_DEFAULT, &matPrimal);
884
      MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
885
886

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
887
      MatDestroy(&matBPi);
888

889
      KSPCreate(mpiCommGlobal, &ksp_schur_primal);
890
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
891
		      SAME_NONZERO_PATTERN);
892
893
894
895
896
897
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPPREONLY);
      PC pc_schur_primal;      
      KSPGetPC(ksp_schur_primal, &pc_schur_primal);
      PCSetType(pc_schur_primal, PCLU);
      PCFactorSetMatSolverPackage(pc_schur_primal, MATSOLVERMUMPS);
898
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
899

Thomas Witkowski's avatar
Thomas Witkowski committed
900
      if (printTimings) {
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
901
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
902
903
904
905
906
907
908
909
910
911
	MatInfo minfo;
	MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
	MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);
	
	MSG("FETI-DP timing 06: %.5f seconds (creation of schur primal matrix)\n",
	    MPI::Wtime() - wtime);

	wtime = MPI::Wtime();
	KSPSetUp(ksp_schur_primal);
	KSPSetUpOnBlocks(ksp_schur_primal);
Thomas Witkowski's avatar
BLUB    
Thomas Witkowski committed
912
	MPI::COMM_WORLD.Barrier();
Thomas Witkowski's avatar
Thomas Witkowski committed
913
914
915
	MSG("FETI-DP timing 07: %.5f seconds (factorization of primal schur matrix).\n",
	    MPI::Wtime() - wtime);
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
916
    }
917
918
919