PetscSolverFeti.cc 53.3 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 "io/VtkWriter.h"
18 19 20 21 22

namespace AMDiS {

  using namespace std;

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

    MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
33
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
34 35 36 37 38 39 40 41 42 43 44
    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)
  {
45 46
    //    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)
47 48 49

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

Thomas Witkowski's avatar
Thomas Witkowski committed
52 53 54
    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
55

Thomas Witkowski's avatar
Thomas Witkowski committed
56
    MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
57
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
58 59 60
    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);
61

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

    return 0;
  }


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

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


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

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

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

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_duals[j] = local_b[i];
93 94

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


98
    // === K_DD - K_DI inv(K_II) K_ID ===
99

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

102
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
103
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
    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);

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_b[i] = local_duals[j];

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


140
    // === Restriction of the B nodes to the boundary nodes. ===
141

142 143 144 145
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
146

147
    PetscScalar *local_b, *local_duals;
148
    VecGetArray(data->tmp_vec_b, &local_b);
149
    VecGetArray(data->tmp_vec_duals0, &local_duals);
150

151 152
    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_duals[j] = local_b[i];
153 154

    VecRestoreArray(data->tmp_vec_b, &local_b);
155 156 157 158 159 160 161
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

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

162

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

165 166 167 168 169
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

    for (int i = nLocalB - nLocalDuals, j = 0; i < nLocalB; i++, j++)
      local_b[i] = local_duals[j];
170

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

174 175

    // Multiply with scaled Lagrange constraint matrix.
176 177 178 179 180 181
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


182 183
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
Thomas Witkowski's avatar
Thomas Witkowski committed
184 185
      nComponents(-1),
      schurPrimalSolver(0)
186 187 188 189 190 191
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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


212
  void PetscSolverFeti::updateDofData()
213 214
  {
    FUNCNAME("PetscSolverFeti::updateDofData()");
215 216 217

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
218
   
219 220 221 222 223
    createPrimals();

    createDuals();

    createLagrange();
224
    
225
    createIndexB();
226 227 228
  }


229
  void PetscSolverFeti::createPrimals()
230
  {
231
    FUNCNAME("PetscSolverFeti::createPrimals()");  
232

233 234
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

235 236 237
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

238 239 240
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
241
    DofContainerSet& vertices = 
242
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
243 244 245 246
    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);
247

248 249 250 251

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

252
    primalDofMap.addFeSpace(feSpace);
253
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
254 255
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
256

257
    primalDofMap[feSpace].update();
258

259
    MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
260 261
	primalDofMap[feSpace].nRankDofs, 
	primalDofMap[feSpace].nOverallDofs);
262

263 264 265 266

    // === Communicate primal's global index from ranks that own the     ===
    // === primals to ranks that contain this primals but are not owning ===
    // === them.                                                         ===
267 268
    
    typedef map<int, map<const FiniteElemSpace*, DofContainer> >::iterator it_type;
269

270
    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
271 272 273
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
274 275 276 277
	if (primalDofMap[feSpace].isSet(it.getDofIndex())) {
	  DegreeOfFreedom d = primalDofMap[feSpace][it.getDofIndex()];
	  stdMpi.getSendData(it.getRank()).push_back(d);
	}
278

279 280
    stdMpi.updateSendDataSize();

281 282
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	 !it.end(); it.nextRank()) {
283
      bool recvFromRank = false;
284 285 286
      for (; !it.endDofIter(); it.nextDof()) {
	if (primals.count(it.getDofIndex()) && 
	    meshDistributor->getIsRankDof(feSpace, it.getDofIndex()) == false) {
287 288 289
	  recvFromRank = true;
	  break;
	}
290
      }
291 292

      if (recvFromRank) 
293
	stdMpi.recv(it.getRank());
294
    }
295

296 297
    stdMpi.startCommunication();

298 299
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
	 !it.end(); it.nextRank()) {
300
      int i = 0;
301 302
      for (; !it.endDofIter(); it.nextDof()) {
	if (primals.count(it.getDofIndex()) && 
303 304 305 306
	    meshDistributor->getIsRankDof(feSpace, it.getDofIndex()) == false) {
	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[i++];
	  primalDofMap[feSpace].insert(it.getDofIndex(), d);
	}
307 308 309
      }
    }

310
    TEST_EXIT_DBG(primals.size() == primalDofMap[feSpace].size())
311
      ("Number of primals %d, but number of global primals on this rank is %d!\n",
312
       primals.size(), primalDofMap[feSpace].size());
313 314 315 316 317 318
  }


  void PetscSolverFeti::createDuals()
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
319 320

    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
321
    
322 323
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
324 325 326

    boundaryDofRanks.clear();

327 328 329
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
330
	// If DOF is not primal, i.e., its a dual node
331
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
332 333
	  boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
334 335
	}
      }
336
	
337 338 339 340

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

341
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
342 343 344 345

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
346
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false)
347
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
348 349 350

    stdMpi.updateSendDataSize();

351 352
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
353
      bool recvFromRank = false;
354
      for (; !it.endDofIter(); it.nextDof()) {
355
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
356 357 358
	  recvFromRank = true;
	  break;
	}
359
      }
360 361

      if (recvFromRank)
362
	stdMpi.recv(it.getRank());
363
    }
364

365 366
    stdMpi.startCommunication();

367 368
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
369
      int i = 0;
370
      for (; !it.endDofIter(); it.nextDof())
371
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false)
372 373
	  boundaryDofRanks[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[i++];
374 375 376 377
    }

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

378
    dualDofMap.addFeSpace(feSpace);
379

380
    int nRankAllDofs = feSpace->getAdmin()->getUsedDofs();
381
    nRankB = nRankAllDofs - primalDofMap[feSpace].size();
382 383 384 385 386
    nOverallB = 0;
    rStartB = 0;
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankB, rStartB, nOverallB);
    DofContainer allBoundaryDofs;
387
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
388 389 390
    int nRankInteriorDofs = nRankAllDofs - allBoundaryDofs.size();

    for (DofContainer::iterator it = allBoundaryDofs.begin();
391 392 393
	 it != allBoundaryDofs.end(); ++it)
      if (primalDofMap[feSpace].isSet(**it) == false)
	dualDofMap[feSpace].insertRankDof(**it);
394

395 396
    dualDofMap[feSpace].update(false);
    dualDofMap[feSpace].addOffset(rStartB + nRankInteriorDofs);
397

398
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
399 400
	dualDofMap[feSpace].nRankDofs, 
	dualDofMap[feSpace].nOverallDofs);
401 402 403 404 405 406 407
  }

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

408 409 410
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

411 412 413 414 415 416 417 418 419
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
    dofFirstLagrange.addFeSpace(feSpace);

    int nRankLagrange = 0;
    DofMapping& dualMap = dualDofMap[feSpace].getMap();
    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
	dofFirstLagrange[feSpace].insert(it->first, nRankLagrange);
	int degree = boundaryDofRanks[it->first].size();
420 421 422
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }
423 424 425
    dofFirstLagrange[feSpace].nRankDofs = nRankLagrange;
    dofFirstLagrange[feSpace].update();
    
426
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
427 428
	dofFirstLagrange[feSpace].nRankDofs, 
	dofFirstLagrange[feSpace].nOverallDofs);
429 430


431
    // === Communicate dofFirstLagrange to all other ranks. ===
432 433

    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
434

435
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
436
	 !it.end(); it.nextRank()) {
437
      for (; !it.endDofIter(); it.nextDof())
438 439
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
	  TEST_EXIT_DBG(dofFirstLagrange[feSpace].isSet(it.getDofIndex()))
440
	    ("Should not happen!\n");
441 442
	  DegreeOfFreedom d = dofFirstLagrange[feSpace][it.getDofIndex()];
	  stdMpi.getSendData(it.getRank()).push_back(d);
443
	}
444
    }
445

446 447
    stdMpi.updateSendDataSize();

448
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
449
	 !it.end(); it.nextRank()) {
450
      bool recvData = false;
451
      for (; !it.endDofIter(); it.nextDof())
452
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
453 454 455 456 457
	  recvData = true;
	  break;
	}
	  
      if (recvData)
458
	stdMpi.recv(it.getRank());
459 460 461 462
    }

    stdMpi.startCommunication();

463
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
464
	 !it.end(); it.nextRank()) {
465
      int counter = 0;
466 467 468 469 470 471
      for (; !it.endDofIter(); it.nextDof()) {
	if (primalDofMap[feSpace].isSet(it.getDofIndex()) == false) {
	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
	  dofFirstLagrange[feSpace].insert(it.getDofIndex(), d); 
	}
      }
472
    }
473 474 475 476 477 478 479
  }


  void PetscSolverFeti::createIndexB()
  {
    FUNCNAME("PetscSolverFeti::createIndeB()");

480
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
481
    localIndexB.clear();
482
    DOFAdmin* admin = feSpace->getAdmin();
483 484 485 486

    // === 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.                                 ===
487 488
    
    for (int i = 0; i < admin->getUsedSize(); i++)
489 490 491 492
      if (admin->isDofFree(i) == false && 
	  primalDofMap[feSpace].isSet(i) == false && 
	  dualDofMap[feSpace].isSet(i) == false)
	localIndexB[i] = -1;
493

494 495 496

    // === Get correct index for all interior nodes. ===

497
    nLocalInterior = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
498 499 500
    for (DofMapping::iterator it = localIndexB.begin(); 
	 it != localIndexB.end(); ++it) {
      it->second = nLocalInterior;
501
      nLocalInterior++;
502
    }
503
    nLocalDuals = dualDofMap[feSpace].size();
504

505 506 507
    TEST_EXIT_DBG(nLocalInterior + 
		  primalDofMap[feSpace].size() + 
		  dualDofMap[feSpace].size() == 
508 509 510
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

511 512 513

    // === And finally, add the global indicies of all dual nodes. ===

514 515 516
    for (DofMapping::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
      localIndexB[it->first] = it->second - rStartB;
517 518 519
  }


520
  void PetscSolverFeti::createMatLagrange()
521 522 523
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

524 525
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

526 527
    // === Create distributed matrix for Lagrange constraints. ===

528
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
529
		    dofFirstLagrange[feSpace].nRankDofs * nComponents, 
530
		    nRankB * nComponents,
531
		    dofFirstLagrange[feSpace].nOverallDofs * nComponents, 
532
		    nOverallB * nComponents,
533 534 535
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

536 537 538 539 540 541 542
    // === 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.                                                     ===

543 544 545 546 547 548
    DofMapping &dualMap = dualDofMap[feSpace].getMap();
    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      TEST_EXIT_DBG(dofFirstLagrange[feSpace].isSet(it->first))
	("Should not happen!\n");
      TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
	("Should not happen!\n");
549

550
      // Global index of the first Lagrange constriant for this node.
551
      int index = dofFirstLagrange[feSpace][it->first];
552
      // Copy set of all ranks that contain this dual node.
553 554
      vector<int> W(boundaryDofRanks[it->first].begin(), 
		    boundaryDofRanks[it->first].end());
555
      // Number of ranks that contain this dual node.
556 557 558 559
      int degree = W.size();

      for (int i = 0; i < degree; i++) {
	for (int j = i + 1; j < degree; j++) {
560 561
	  if (W[i] == mpiRank || W[j] == mpiRank) {
	    // Set the constraint for all components of the system.
562 563
	    for (int k = 0; k < nComponents; k++) {
	      int rowIndex = index * nComponents + k;
564
	      int colIndex = it->second * nComponents + k;
565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
			  INSERT_VALUES);
	    }
	  }

	  index++;
	}
      }
    }

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


581
  void PetscSolverFeti::createSchurPrimalKsp()
582 583 584
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

585 586
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
590 591 592 593 594 595 596 597 598
      schurPrimalData.mat_primal_primal = &mat_primal_primal;
      schurPrimalData.mat_primal_b = &mat_primal_b;
      schurPrimalData.mat_b_primal = &mat_b_primal;
      schurPrimalData.fetiSolver = this;
      
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
599
		   primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
600 601 602
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
603 604
		     primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nRankDofs * nComponents,
		     primalDofMap[feSpace].nOverallDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
605 606 607 608 609 610 611 612 613 614
		     &schurPrimalData, 
		     &mat_schur_primal);
      MatShellSetOperation(mat_schur_primal, MATOP_MULT, 
			   (void(*)(void))petscMultMatSchurPrimal);
      
      KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal);
      KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal, SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_");
      KSPSetFromOptions(ksp_schur_primal);
    } else {
615 616
      MSG("Create direct schur primal solver!\n");

617 618 619 620
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

621 622
      int nRowsRankPrimal = primalDofMap[feSpace].nRankDofs * nComponents;
      int nRowsOverallPrimal = primalDofMap[feSpace].nOverallDofs * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
623 624
      int nRowsRankB = nRankB * nComponents;
      int nRowsOverallB = nOverallB * nComponents;
625

Thomas Witkowski's avatar
Thomas Witkowski committed
626 627 628 629 630 631
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
632

633 634 635
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
636

637 638 639 640 641 642 643 644 645 646
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

      for (int i = 0; i < nRankB * nComponents; i++) {
	PetscInt row = rStartB * nComponents + i;
	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
647
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
648 649 650 651 652 653
      }

      int maxLocalPrimal = mat_b_primal_cols.size();
      mpi::globalMax(maxLocalPrimal);

      TEST_EXIT(mat_b_primal_cols.size() ==
654
		primalDofMap[feSpace].size() * nComponents)("Should not happen!\n");
655 656 657 658 659 660 661 662 663 664 665 666 667 668
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
	VecCreateSeq(PETSC_COMM_SELF, nRankB * nComponents, &tmpVec);

 	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
669
	PetscScalar *tmpValues;
670 671
	VecGetArray(tmpVec, &tmpValues);
	for (int i  = 0; i < nRankB * nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
672 673 674 675 676
	  MatSetValue(matBPi, 
		      rStartB * nComponents + i,
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
677 678 679 680 681
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
682 683 684
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
685 686 687
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
688
      MatDestroy(&matBPi);
689 690 691 692 693 694 695 696 697 698

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

      KSPCreate(PETSC_COMM_WORLD, &ksp_schur_primal);
      KSPSetOperators(ksp_schur_primal, mat_primal_primal, 
		      mat_primal_primal, SAME_NONZERO_PATTERN);
      KSPSetOptionsPrefix(ksp_schur_primal, "solver_sp_");
      KSPSetFromOptions(ksp_schur_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
699

700 701
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
702
    }
703 704 705 706 707 708 709
  }


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

710 711 712 713 714
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
715

716 717
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
718

719 720 721 722 723
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
724 725 726 727 728 729 730
  }


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

731 732
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

733 734
    // === Create FETI-DP solver object. ===

735 736 737
    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
738
    fetiData.fetiSolver = this;
739
    fetiData.ksp_schur_primal = &ksp_schur_primal;
740

741 742
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
743 744
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
745 746
		 dofFirstLagrange[feSpace].nRankDofs * nComponents, 
		 dofFirstLagrange[feSpace].nOverallDofs * nComponents,
Thomas Witkowski's avatar
Thomas Witkowski committed
747
		 &(fetiData.tmp_vec_lagrange));
748
    VecCreateMPI(PETSC_COMM_WORLD, 
749
		 primalDofMap[feSpace].nRankDofs * nComponents, primalDofMap[feSpace].nOverallDofs * nComponents,
750
		 &(fetiData.tmp_vec_primal));
751 752

    MatCreateShell(PETSC_COMM_WORLD,
753 754 755 756
		   dofFirstLagrange[feSpace].nRankDofs * nComponents, 
		   dofFirstLagrange[feSpace].nRankDofs * nComponents,
		   dofFirstLagrange[feSpace].nOverallDofs * nComponents, 
		   dofFirstLagrange[feSpace].nOverallDofs * nComponents,
757
		   &fetiData, &mat_feti);
758 759 760 761 762 763 764
    MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);


    KSPCreate(PETSC_COMM_WORLD, &ksp_feti);
    KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
    KSPSetOptionsPrefix(ksp_feti, "solver_feti_");
    KSPSetFromOptions(ksp_feti);
765 766


767
    // === Create FETI-DP preconditioner object. ===
768

769 770 771 772
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
773

774 775 776
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
777 778
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
779 780 781 782 783 784 785 786 787 788
      KSPSetOptionsPrefix(ksp_interior, "solver_interior_");
      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;
      
789 790 791
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiDirichletPreconData.tmp_vec_b));      
792 793 794 795 796 797 798 799 800 801 802 803 804 805 806
      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));
      
      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;

807 808 809
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiLumpedPreconData.tmp_vec_b));
810 811 812 813 814 815 816 817
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals1));

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
818 819
      break;
    default:
820 821
      break;
    }
822 823 824 825 826 827 828
  }
  

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

829 830
    // === Destroy FETI-DP solver object. ===

831 832 833
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
834
    fetiData.fetiSolver = NULL;
835
    fetiData.ksp_schur_primal = PETSC_NULL;
836

Thomas Witkowski's avatar
Thomas Witkowski committed
837 838
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
839
    VecDestroy(&fetiData.tmp_vec_primal);
840 841
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
842 843


844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871
    // === 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;
872 873
    default:
      break;
874
    }
875 876 877 878 879 880 881 882 883
  }


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

884 885
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

886
    // === Get local part of the solution for B variables. ===
887 888 889 890 891

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


892 893
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
894 895
    
    vector<PetscInt> globalIsIndex, localIsIndex;
896 897
    globalIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
    localIsIndex.reserve(primalDofMap[feSpace].size() * nComponents);
898 899 900

    {
      int counter = 0;
901 902
      for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin();
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932
	for (int i = 0; i < nComponents; i++) {
	  globalIsIndex.push_back(it->second * nComponents + i);
	  localIsIndex.push_back(counter++);
	}
      }
    }
    
    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);

933 934 935
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
936 937 938 939

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

940
    // === And copy from PETSc local vectors to the DOF vectors. ===
941 942 943 944

    for (int i = 0; i < nComponents; i++) {
      DOFVector<double>& dofVec = *(vec.getDOFVector(i));

Thomas Witkowski's avatar
Thomas Witkowski committed
945 946 947
      for (DofMapping::iterator it = localIndexB.begin();
	   it != localIndexB.end(); ++it) {
	int petscIndex = it->second * nComponents + i;
948 949 950 951
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
952 953
      for (DofMapping::iterator it = primalDofMap[feSpace].getMap().begin();
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
954 955 956 957 958