PetscSolverFeti.cc 45.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13
//
// 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.


#include "parallel/PetscSolverFeti.h"
14
#include "parallel/PetscSolverFetiStructs.h"
15 16
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
17
#include "parallel/SubDomainSolver.h"
18
#include "io/VtkWriter.h"
19 20 21 22 23

namespace AMDiS {

  using namespace std;

24 25 26 27 28 29 30
  // 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);
31
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
32 33

    MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
34
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
35 36 37 38 39 40 41 42 43 44 45
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
    MatMult(*(data->mat_primal_primal), x, y);
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
46 47
    //    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)
48 49 50

    void *ctx;
    MatShellGetContext(mat, &ctx);
51
    FetiData* data = static_cast<FetiData*>(ctx);
52

Thomas Witkowski's avatar
Thomas Witkowski committed
53 54 55
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
56

Thomas Witkowski's avatar
Thomas Witkowski committed
57
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
58
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
59 60 61
    MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b);
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
62

Thomas Witkowski's avatar
Thomas Witkowski committed
63
    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
64 65 66 67 68

    return 0;
  }


69
  // y = PC * x
70
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
71
  {
72
    // Get data for the preconditioner
73 74
    void *ctx;
    PCShellGetContext(pc, &ctx);
75
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
76

77
    // Multiply with scaled Lagrange constraint matrix.
78 79 80
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


81
    // === Restriction of the B nodes to the boundary nodes. ===
82

83 84 85 86
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
87

88 89 90 91
    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

92 93 94
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
95 96

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


100
    // === K_DD - K_DI inv(K_II) K_ID ===
101

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

104
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
105
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
106 107 108 109 110 111 112 113 114 115
    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);

116 117 118
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140

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

    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);
141 142


143
    // === Restriction of the B nodes to the boundary nodes. ===
144

145 146 147 148
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
149

150
    PetscScalar *local_b, *local_duals;
151
    VecGetArray(data->tmp_vec_b, &local_b);
152
    VecGetArray(data->tmp_vec_duals0, &local_duals);
153

Thomas Witkowski's avatar
Thomas Witkowski committed
154 155 156
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
157 158

    VecRestoreArray(data->tmp_vec_b, &local_b);
159 160 161 162 163 164 165
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

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

166

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

169 170 171
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

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

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

179 180

    // Multiply with scaled Lagrange constraint matrix.
181 182 183 184 185 186
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


187 188
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
189 190
      schurPrimalSolver(0),
      multiLevelTest(false)
191 192 193 194 195 196
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
197
      MSG("Create FETI-DP solver with no preconditioner!\n");
198 199
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
200
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
201 202
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
203
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
204 205
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
206 207
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
208
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
209 210 211 212 213

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

    Parameters::get("parallel->multi level test", multiLevelTest);

217 218
    if (multiLevelTest) {      
      subDomainSolver = new SubDomainSolver(meshDistributor, mpiComm, mpiSelfComm);
219
    }
220 221 222
  }


223
  void PetscSolverFeti::updateDofData()
224 225
  {
    FUNCNAME("PetscSolverFeti::updateDofData()");
226

227 228 229 230 231 232
    primalDofMap.clear();
    dualDofMap.clear();
    lagrangeMap.clear();
    localDofMap.clear();
    interiorDofMap.clear();

233 234
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
235

236 237
      createPrimals(feSpace);      
      createDuals(feSpace);
238
      createLagrange(feSpace);      
239 240
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
241

242 243 244 245 246
    primalDofMap.setDofComm(meshDistributor->getSendDofs(),
			    meshDistributor->getRecvDofs());
    lagrangeMap.setDofComm(meshDistributor->getSendDofs(), 
			   meshDistributor->getRecvDofs());

Thomas Witkowski's avatar
Thomas Witkowski committed
247 248 249 250
    primalDofMap.update();
    dualDofMap.update();
    lagrangeMap.update();
    localDofMap.update();
251
    if (fetiPreconditioner != FETI_NONE)
252
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
253 254 255

    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
256
      MSG("FETI-DP data for %d-ith FE space:\n", i);
Thomas Witkowski's avatar
Thomas Witkowski committed
257

258
      MSG("  nRankPrimals = %d   nOverallPrimals = %d\n", 
Thomas Witkowski's avatar
Thomas Witkowski committed
259 260
	  primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs);
      
261
      MSG("  nRankDuals = %d  nOverallDuals = %d\n",
262
	  dualDofMap[feSpace].nRankDofs, dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
263

264
      MSG("  nRankLagrange = %d  nOverallLagrange = %d\n",
265
	  lagrangeMap[feSpace].nRankDofs, lagrangeMap[feSpace].nOverallDofs);
266 267 268 269

      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
270
    }
271 272 273 274 275 276


    // If multi level test, inform sub domain solver about coarse space.
    if (multiLevelTest) {
      subDomainSolver->setCoarseSpace(&primalDofMap);
    }
277 278 279
  }


280
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
281
  {
282
    FUNCNAME("PetscSolverFeti::createPrimals()");  
283

284 285 286
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

287 288 289
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
290
    DofContainerSet& vertices = 
291
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
292 293 294 295
    TEST_EXIT_DBG(vertices.size())("No primal vertices on this rank!\n");
    for (DofContainerSet::iterator it = vertices.begin(); 
	 it != vertices.end(); ++it)
      primals.insert(**it);
296

297 298 299 300

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

301
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
302 303
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
304 305
      else
  	primalDofMap[feSpace].insert(*it);
306 307 308
  }


309
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
310 311
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
312

313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328
    // === Create global index of the dual nodes on each rank. ===

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

    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it)
      if (!isPrimal(feSpace, **it))
	dualDofMap[feSpace].insertRankDof(**it);
  }

  
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

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

332
    boundaryDofRanks[feSpace].clear();
333

334 335 336
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
337
	if (!isPrimal(feSpace, it.getDofIndex())) {
338 339
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
340 341
	}
      }
342
	
343 344 345 346

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

347
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
348 349 350 351

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
352
	if (!isPrimal(feSpace, it.getDofIndex()))
353
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
354 355 356

    stdMpi.updateSendDataSize();

357 358
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
359
      bool recvFromRank = false;
360
      for (; !it.endDofIter(); it.nextDof()) {
361
	if (!isPrimal(feSpace, it.getDofIndex())) {
362 363 364
	  recvFromRank = true;
	  break;
	}
365
      }
366 367

      if (recvFromRank)
368
	stdMpi.recv(it.getRank());
369
    }
370

371 372
    stdMpi.startCommunication();

373 374
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
375
      int i = 0;
376
      for (; !it.endDofIter(); it.nextDof())
377
	if (!isPrimal(feSpace, it.getDofIndex()))
378
	  boundaryDofRanks[feSpace][it.getDofIndex()] = 
379
	    stdMpi.getRecvData(it.getRank())[i++];
380 381 382
    }


383 384 385
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

386
    int nRankLagrange = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
387 388
    map<DegreeOfFreedom, MultiIndex>& dualMap = dualDofMap[feSpace].getMap();
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
389
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
390
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
391
	int degree = boundaryDofRanks[feSpace][it->first].size();
392
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
393 394
      } else {
	lagrangeMap[feSpace].insert(it->first);
395 396
      }
    }
397
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
398 399 400
  }


401
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
402
  {
403
    FUNCNAME("PetscSolverFeti::createIndexB()");
404

405
    DOFAdmin* admin = feSpace->getAdmin();
406 407 408 409

    // === 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.                                 ===
410

411
    int nLocalInterior = 0;    
412
    for (int i = 0; i < admin->getUsedSize(); i++) {
413
      if (admin->isDofFree(i) == false && 
414
	  isPrimal(feSpace, i) == false &&
415
	  dualDofMap[feSpace].isSet(i) == false) {
416 417
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);

418
	if (fetiPreconditioner != FETI_NONE)
419 420 421
	  interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);

	nLocalInterior++;
422
      }
423
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
424
    
425 426
    // === And finally, add the global indicies of all dual nodes. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
427
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualDofMap[feSpace].getMap().begin();
428
	 it != dualDofMap[feSpace].getMap().end(); ++it)
429
      localDofMap[feSpace].insertRankDof(it->first);
430 431 432
  }


433
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
434 435 436
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

437 438
    // === Create distributed matrix for Lagrange constraints. ===

439
    MatCreateMPIAIJ(mpiComm,
440 441
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),	
442 443 444
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

445 446 447 448 449 450 451
    // === 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
452 453 454
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
      map<DegreeOfFreedom, MultiIndex> &dualMap = 
	dualDofMap[feSpaces[k]].getMap();
455

Thomas Witkowski's avatar
Thomas Witkowski committed
456 457
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	   it != dualMap.end(); ++it) {
458
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
459 460 461 462 463 464
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
	int index = lagrangeMap.getMatIndex(k, it->first);
	
	// Copy set of all ranks that contain this dual node.
465 466
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
467 468 469 470 471 472
	// Number of ranks that contain this dual node.
	int degree = W.size();
	
	for (int i = 0; i < degree; i++) {
	  for (int j = i + 1; j < degree; j++) {
	    if (W[i] == mpiRank || W[j] == mpiRank) {
473
	      int colIndex = localDofMap.getMatIndex(k, it->first);
474
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
475
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
476
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
477
	    index++;	      
478 479 480 481 482 483 484 485 486 487
	  }
	}
      }
    }

    MatAssemblyBegin(mat_lagrange, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_lagrange, MAT_FINAL_ASSEMBLY);
  }


488
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
489 490 491
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
495 496 497 498 499
      schurPrimalData.mat_primal_primal = &mat_primal_primal;
      schurPrimalData.mat_primal_b = &mat_primal_b;
      schurPrimalData.mat_b_primal = &mat_b_primal;
      schurPrimalData.fetiSolver = this;
      
500
      VecCreateMPI(mpiComm, 
501
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
502
		   &(schurPrimalData.tmp_vec_b));
503
      VecCreateMPI(mpiComm, 
504
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
505 506
		   &(schurPrimalData.tmp_vec_primal));

507
      MatCreateShell(mpiComm,
508 509
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
510 511 512 513 514
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
515
      KSPCreate(mpiComm, &ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
516
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
517 518
      KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
      KSPSetType(ksp_schur_primal, KSPGMRES);
Thomas Witkowski's avatar
Thomas Witkowski committed
519 520
      KSPSetFromOptions(ksp_schur_primal);
    } else {
521 522
      MSG("Create direct schur primal solver!\n");

523 524 525 526
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

527 528
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
529 530
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
531

Thomas Witkowski's avatar
Thomas Witkowski committed
532
      Mat matBPi;
533
      MatCreateMPIAIJ(mpiComm,
Thomas Witkowski's avatar
Thomas Witkowski committed
534 535 536 537
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
538

539 540 541
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
542

543 544
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

545 546
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
547 548 549 550 551 552
	MatGetRow(mat_b_primal, row, &nCols, &cols, &values);

	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]));

Thomas Witkowski's avatar
Thomas Witkowski committed
553
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
554 555
      }

556 557
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
		primalDofMap.getLocalDofs())
558
	("Should not happen!\n");
559

560 561 562
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
563
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
564 565 566 567 568 569 570 571 572 573

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

       	KSPSolve(ksp_b, tmpVec, tmpVec);

Thomas Witkowski's avatar
Thomas Witkowski committed
574
	PetscScalar *tmpValues;
575
	VecGetArray(tmpVec, &tmpValues);
576
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
577
	  MatSetValue(matBPi, 
578
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
579 580 581
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
582 583 584 585 586
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
587 588 589
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
590 591 592
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
593
      MatDestroy(&matBPi);
594 595 596 597 598

      MatInfo minfo;
      MatGetInfo(mat_primal_primal, MAT_GLOBAL_SUM, &minfo);
      MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);

599
      KSPCreate(mpiComm, &ksp_schur_primal);
600 601
      KSPSetOperators(ksp_schur_primal, mat_primal_primal, 
		      mat_primal_primal, SAME_NONZERO_PATTERN);
602 603 604 605 606 607
      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);
608
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
609

610 611
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
612
    }
613 614 615 616 617 618 619
  }


  void PetscSolverFeti::destroySchurPrimalKsp()
  {
    FUNCNAME("PetscSolverFeti::destroySchurPrimal()");

620 621 622 623 624
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
625

626 627
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
628

629 630 631 632 633
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
634 635 636
  }


637
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
638 639 640
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

641 642
    // === Create FETI-DP solver object. ===

643 644 645
    fetiData.mat_primal_b = &mat_primal_b;
    fetiData.mat_b_primal = &mat_b_primal;
    fetiData.mat_lagrange = &mat_lagrange;
Thomas Witkowski's avatar
Thomas Witkowski committed
646
    fetiData.fetiSolver = this;
647
    fetiData.ksp_schur_primal = &ksp_schur_primal;
648

649
    VecCreateMPI(mpiComm, 
650
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
651
		 &(fetiData.tmp_vec_b));
652
    VecCreateMPI(mpiComm,
653
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
654
		 &(fetiData.tmp_vec_lagrange));
655
    VecCreateMPI(mpiComm, 
656
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
657
		 &(fetiData.tmp_vec_primal));
658

659
    MatCreateShell(mpiComm,
660 661
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
662
		   &fetiData, &mat_feti);
663 664 665
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);


666
    KSPCreate(mpiComm, &ksp_feti);
667
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
668 669 670
    KSPSetOptionsPrefix(ksp_feti, "feti_");
    KSPSetType(ksp_feti, KSPGMRES);
    KSPSetTolerances(ksp_feti, 0, 1e-8, 1e+3, 1000);
671
    KSPSetFromOptions(ksp_feti);
672 673


674
    // === Create FETI-DP preconditioner object. ===
675

676 677 678 679
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
680

681 682 683
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
684 685
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
686 687 688 689 690 691
      KSPSetOptionsPrefix(ksp_interior, "precon_interior_");
      KSPSetType(ksp_interior, KSPPREONLY);
      PC pc_interior;
      KSPGetPC(ksp_interior, &pc_interior);
      PCSetType(pc_interior, PCLU);
      PCFactorSetMatSolverPackage(pc_interior, MATSOLVERUMFPACK);
692 693 694 695 696 697 698 699 700
      KSPSetFromOptions(ksp_interior);
            
      fetiDirichletPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiDirichletPreconData.mat_interior_interior = &mat_interior_interior;
      fetiDirichletPreconData.mat_duals_duals = &mat_duals_duals;
      fetiDirichletPreconData.mat_interior_duals = &mat_interior_duals;
      fetiDirichletPreconData.mat_duals_interior = &mat_duals_interior;
      fetiDirichletPreconData.ksp_interior = &ksp_interior;
      
701
      VecCreateMPI(mpiComm, 
702
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
703
		   &(fetiDirichletPreconData.tmp_vec_b));      
704 705 706 707 708 709
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_duals1));
      MatGetVecs(mat_interior_interior, PETSC_NULL, 
		 &(fetiDirichletPreconData.tmp_vec_interior));
710 711 712 713 714 715 716 717 718 719 720 721 722

      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	map<DegreeOfFreedom, MultiIndex> &dualMap = 
	  dualDofMap[feSpaces[i]].getMap();
	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	     it != dualMap.end(); ++it) {
	  DegreeOfFreedom d = it->first;
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	  fetiDirichletPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

723 724 725 726 727 728 729 730 731 732 733
      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiDirichletPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiDirichletPrecon);
      
      break;

    case FETI_LUMPED:
      fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
      fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals;

Thomas Witkowski's avatar
Thomas Witkowski committed
734 735 736 737 738 739 740 741 742 743 744 745
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	map<DegreeOfFreedom, MultiIndex> &dualMap = 
	  dualDofMap[feSpaces[i]].getMap();
	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	     it != dualMap.end(); ++it) {
	  DegreeOfFreedom d = it->first;
	  int matIndexLocal = localDofMap.getLocalMatIndex(i, d);
	  int matIndexDual = dualDofMap.getLocalMatIndex(i, d);
	  fetiLumpedPreconData.localToDualMap[matIndexLocal] = matIndexDual;
	}
      }

746
      VecCreateMPI(mpiComm, 
747 748
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
749
		   &(fetiLumpedPreconData.tmp_vec_b));
750 751 752 753
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
754 755 756 757 758 759

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
760 761
      break;
    default:
762 763
      break;
    }
764 765 766 767 768 769 770
  }
  

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

771 772
    // === Destroy FETI-DP solver object. ===

773 774 775
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
776
    fetiData.fetiSolver = NULL;
777
    fetiData.ksp_schur_primal = PETSC_NULL;
778

Thomas Witkowski's avatar
Thomas Witkowski committed
779 780
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
781
    VecDestroy(&fetiData.tmp_vec_primal);
782 783
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
784 785


786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813
    // === Destroy FETI-DP preconditioner object. ===

    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPDestroy(&ksp_interior);

      fetiDirichletPreconData.mat_lagrange_scaled = NULL;
      fetiDirichletPreconData.mat_interior_interior = NULL;
      fetiDirichletPreconData.mat_duals_duals = NULL;
      fetiDirichletPreconData.mat_interior_duals = NULL;
      fetiDirichletPreconData.mat_duals_interior = NULL;
      fetiDirichletPreconData.ksp_interior = NULL;
      
      VecDestroy(&fetiDirichletPreconData.tmp_vec_b);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_duals0);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_duals1);
      VecDestroy(&fetiDirichletPreconData.tmp_vec_interior);
      MatDestroy(&mat_lagrange_scaled);
      break;

    case FETI_LUMPED:
      fetiLumpedPreconData.mat_lagrange_scaled = NULL;
      fetiLumpedPreconData.mat_duals_duals = NULL;

      VecDestroy(&fetiLumpedPreconData.tmp_vec_b);
      VecDestroy(&fetiLumpedPreconData.tmp_vec_duals0);
      VecDestroy(&fetiLumpedPreconData.tmp_vec_duals1);
      break;
814 815
    default:
      break;
816
    }
817 818 819 820 821 822 823 824 825
  }


  void PetscSolverFeti::recoverSolution(Vec &vec_sol_b,
					Vec &vec_sol_primal,
					SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::recoverSolution()");

826
    // === Get local part of the solution for B variables. ===
827 828 829 830 831

    PetscScalar *localSolB;
    VecGetArray(vec_sol_b, &localSolB);


832 833
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
834
    
Thomas Witkowski's avatar
Thomas Witkowski committed
835 836
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(&vec);

837
    vector<PetscInt> globalIsIndex, localIsIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
838 839
    globalIsIndex.reserve(primalDofMap.getLocalDofs());
    localIsIndex.reserve(primalDofMap.getLocalDofs());
840 841

    {
Thomas Witkowski's avatar
Thomas Witkowski committed
842 843 844 845 846 847 848 849
      int cnt = 0;
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	map<DegreeOfFreedom, MultiIndex>& dofMap = 
	  primalDofMap[feSpaces[i]].getMap();
	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin();
	     it != dofMap.end(); ++it) {
	  globalIsIndex.push_back(primalDofMap.getMatIndex(i, it->first));
	  localIsIndex.push_back(cnt++);
850 851
	}
      }
852
      
Thomas Witkowski's avatar
Thomas Witkowski committed
853
      TEST_EXIT_DBG(cnt == primalDofMap.getLocalDofs())("Should not happen!\n");
854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878
    }
    
    IS globalIs, localIs;
    ISCreateGeneral(PETSC_COMM_SELF, 
		    globalIsIndex.size(), 
		    &(globalIsIndex[0]),
		    PETSC_USE_POINTER,
		    &globalIs);

    ISCreateGeneral(PETSC_COMM_SELF, 
		    localIsIndex.size(), 
		    &(localIsIndex[0]),
		    PETSC_USE_POINTER,
		    &localIs);

    Vec local_sol_primal;
    VecCreateSeq(PETSC_COMM_SELF, localIsIndex.size(), &local_sol_primal);

    VecScatter primalScatter;
    VecScatterCreate(vec_sol_primal, globalIs, local_sol_primal, localIs, &primalScatter);
    VecScatterBegin(primalScatter, vec_sol_primal, local_sol_primal, 
		    INSERT_VALUES, SCATTER_FORWARD);
    VecScatterEnd(primalScatter, vec_sol_primal, local_sol_primal, 
		  INSERT_VALUES, SCATTER_FORWARD);

879 880 881
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
882 883 884 885

    PetscScalar *localSolPrimal;
    VecGetArray(local_sol_primal, &localSolPrimal);

886
    // === And copy from PETSc local vectors to the DOF vectors. ===
887

Thomas Witkowski's avatar
Thomas Witkowski committed
888 889
    int cnt = 0;
    for (int i = 0; i < vec.getSize(); i++) {
890 891
      DOFVector<double>& dofVec = *(vec.getDOFVector(i));

Thomas Witkowski's avatar
Thomas Witkowski committed
892 893
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = localDofMap[feSpaces[i]].getMap().begin();
	   it != localDofMap[feSpaces[i]].getMap().end(); ++it) {
894
	int petscIndex = localDofMap.getLocalMatIndex(i, it->first);
895 896 897
	dofVec[it->first] = localSolB[petscIndex];
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
898 899 900
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpaces[i]].getMap().begin();
	   it != primalDofMap[feSpaces[i]].getMap().end(); ++it)
	dofVec[it->first] = localSolPrimal[cnt++];
901 902 903 904
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
905
    VecDestroy(&local_sol_primal);
906 907 908
  }


909
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
910 911
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");