Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

PetscSolverFeti.cc 46.9 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
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
221

222 223
      createPrimals(feSpace);      
      createDuals(feSpace);
224
      createLagrange(feSpace);      
225 226
      createIndexB(feSpace);
    }
227 228 229
  }


230
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
231
  {
232
    FUNCNAME("PetscSolverFeti::createPrimals()");  
233

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

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

247 248 249 250

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

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

255
    primalDofMap[feSpace].update();
256

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

261 262 263 264

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

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

277 278
    stdMpi.updateSendDataSize();

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

      if (recvFromRank) 
291
	stdMpi.recv(it.getRank());
292
    }
293

294 295
    stdMpi.startCommunication();

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

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


314
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
315 316
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
317

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

    boundaryDofRanks.clear();

323 324 325
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
326
	if (!isPrimal(feSpace, it.getDofIndex())) {
327 328
	  boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
329 330
	}
      }
331
	
332 333 334 335

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

336
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
337 338 339 340

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
341
	if (!isPrimal(feSpace, it.getDofIndex()))
342
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
343 344 345

    stdMpi.updateSendDataSize();

346 347
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
348
      bool recvFromRank = false;
349
      for (; !it.endDofIter(); it.nextDof()) {
350
	if (!isPrimal(feSpace, it.getDofIndex())) {
351 352 353
	  recvFromRank = true;
	  break;
	}
354
      }
355 356

      if (recvFromRank)
357
	stdMpi.recv(it.getRank());
358
    }
359

360 361
    stdMpi.startCommunication();

362 363
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
364
      int i = 0;
365
      for (; !it.endDofIter(); it.nextDof())
366
	if (!isPrimal(feSpace, it.getDofIndex()))
367 368
	  boundaryDofRanks[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[i++];
369 370 371 372 373
    }

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

    DofContainer allBoundaryDofs;
374
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
375 376

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

381
    dualDofMap[feSpace].update(false);
382

383
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
384 385
	dualDofMap[feSpace].nRankDofs, 
	dualDofMap[feSpace].nOverallDofs);
386 387 388
  }

  
389
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
390 391 392
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

393 394 395
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

396 397 398 399
    int nRankLagrange = 0;
    DofMapping& dualMap = dualDofMap[feSpace].getMap();
    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
400
	lagrangeMap[feSpace].insert(it->first, nRankLagrange);
401
	int degree = boundaryDofRanks[it->first].size();
402 403 404
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }
405
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
406
    lagrangeMap[feSpace].update();
407
    
408
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
409 410
	lagrangeMap[feSpace].nRankDofs, 
	lagrangeMap[feSpace].nOverallDofs);
411 412


413
    // === Communicate lagrangeMap to all other ranks. ===
414 415

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

417
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
418
	 !it.end(); it.nextRank()) {
419
      for (; !it.endDofIter(); it.nextDof())
420
	if (!isPrimal(feSpace, it.getDofIndex())) {
421
	  DegreeOfFreedom d = lagrangeMap[feSpace][it.getDofIndex()];
422
	  stdMpi.getSendData(it.getRank()).push_back(d);
423
	}
424
    }
425

426 427
    stdMpi.updateSendDataSize();

428
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
429
	 !it.end(); it.nextRank()) {
430
      bool recvData = false;
431
      for (; !it.endDofIter(); it.nextDof())
432
	if (!isPrimal(feSpace, it.getDofIndex())) {
433 434 435 436 437
	  recvData = true;
	  break;
	}
	  
      if (recvData)
438
	stdMpi.recv(it.getRank());
439 440 441 442
    }

    stdMpi.startCommunication();

443
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace);
444
	 !it.end(); it.nextRank()) {
445
      int counter = 0;
446
      for (; !it.endDofIter(); it.nextDof()) {
447
	if (!isPrimal(feSpace, it.getDofIndex())) {
448
	  DegreeOfFreedom d = stdMpi.getRecvData(it.getRank())[counter++];
449
	  lagrangeMap[feSpace].insert(it.getDofIndex(), d); 
450 451
	}
      }
452
    }
453 454 455
  }


456
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
457
  {
458
    FUNCNAME("PetscSolverFeti::createIndexB()");
459

460
    DOFAdmin* admin = feSpace->getAdmin();
461 462 463 464

    // === 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.                                 ===
465 466 467

    nLocalInterior = 0;    
    for (int i = 0; i < admin->getUsedSize(); i++) {
468
      if (admin->isDofFree(i) == false && 
469
	  isPrimal(feSpace, i) == false &&
470 471 472 473
	  dualDofMap[feSpace].isSet(i) == false) {
	localDofMap[feSpace].insertRankDof(i);
	nLocalInterior++;
      }
474 475
    }

476 477
    localDofMap[feSpace].update(false);

478 479
    TEST_EXIT_DBG(nLocalInterior + 
		  primalDofMap[feSpace].size() + 
480
		  dualDofMap[feSpace].size() == 
481 482 483
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

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

486 487
    for (DofMapping::iterator it = dualDofMap[feSpace].getMap().begin();
	 it != dualDofMap[feSpace].getMap().end(); ++it)
488 489 490 491 492 493
      localDofMap[feSpace].insertRankDof(it->first);

    localDofMap[feSpace].update(false);

    dualDofMap[feSpace].addOffset(localDofMap[feSpace].rStartDofs + 
				  nLocalInterior);
494 495 496
  }


497
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
498 499 500
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

501 502
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

503 504
    // === Create distributed matrix for Lagrange constraints. ===

505
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
506 507
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),	
508 509 510
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

511 512 513 514 515 516 517
    // === 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.                                                     ===

518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536
    DofMapping &dualMap = dualDofMap[feSpace].getMap();
    for (DofMapping::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
      TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
	("Should not happen!\n");

      // Global index of the first Lagrange constriant for this node.
      int index = lagrangeMap[feSpace][it->first];
      // Copy set of all ranks that contain this dual node.
      vector<int> W(boundaryDofRanks[it->first].begin(), 
		    boundaryDofRanks[it->first].end());
      // 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) {
	    // Set the constraint for all components of the system.
	    for (int k = 0; k < nComponents; k++) {
	      int rowIndex = index * nComponents + k;
537
	      int colIndex = it->second * nComponents + k;
538
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
539 540
	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
			  INSERT_VALUES);
541 542
	    }
	  }
543 544

	  index++;
545 546 547 548 549 550 551 552 553
	}
      }
    }

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


554
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
555 556 557
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

558 559
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
563 564 565 566 567 568
      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, 
569
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
570 571
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
572
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
573 574 575
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
576 577
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
578 579 580 581 582 583 584 585 586 587
		     &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 {
588 589
      MSG("Create direct schur primal solver!\n");

590 591 592 593
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

594 595
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
596 597
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
598

Thomas Witkowski's avatar
Thomas Witkowski committed
599 600 601 602 603 604
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
605

606 607 608
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
609

610 611
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

612 613
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
614 615 616 617 618 619
	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
620
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
621 622 623 624 625 626
      }

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

      TEST_EXIT(mat_b_primal_cols.size() ==
627 628
		primalDofMap[feSpace].size() * nComponents)
	("Should not happen!\n");
629 630 631
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
632
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
633 634 635 636 637 638 639 640 641 642

 	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
643
	PetscScalar *tmpValues;
644
	VecGetArray(tmpVec, &tmpValues);
645
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
646
	  MatSetValue(matBPi, 
647
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
648 649 650
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
651 652 653 654 655
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
656 657 658
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
659 660 661
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
662
      MatDestroy(&matBPi);
663 664 665 666 667 668 669 670 671 672

      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
673

674 675
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
676
    }
677 678 679 680 681 682 683
  }


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

684 685 686 687 688
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
689

690 691
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
692

693 694 695 696 697
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
698 699 700
  }


701
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
702 703 704
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

705 706
    // === Create FETI-DP solver object. ===

707 708 709
    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
710
    fetiData.fetiSolver = this;
711
    fetiData.ksp_schur_primal = &ksp_schur_primal;
712

713
    VecCreateMPI(PETSC_COMM_WORLD, 
714
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
715 716
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
717
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
718
		 &(fetiData.tmp_vec_lagrange));
719
    VecCreateMPI(PETSC_COMM_WORLD, 
720
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
721
		 &(fetiData.tmp_vec_primal));
722 723

    MatCreateShell(PETSC_COMM_WORLD,
724 725
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
726
		   &fetiData, &mat_feti);
727 728 729 730 731 732 733
    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);
734 735


736
    // === Create FETI-DP preconditioner object. ===
737

738 739 740 741
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
742

743 744 745
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
746 747
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
748 749 750 751 752 753 754 755 756 757
      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;
      
758
      VecCreateMPI(PETSC_COMM_WORLD, 
759
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
760
		   &(fetiDirichletPreconData.tmp_vec_b));      
761 762 763 764 765 766
      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));
767 768 769 770 771 772 773 774 775 776 777 778
      
      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;

779
      VecCreateMPI(PETSC_COMM_WORLD, 
780 781
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
782
		   &(fetiLumpedPreconData.tmp_vec_b));
783 784 785 786
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
787 788 789 790 791 792

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
793 794
      break;
    default:
795 796
      break;
    }
797 798 799 800 801 802 803
  }
  

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

804 805
    // === Destroy FETI-DP solver object. ===

806 807 808
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
809
    fetiData.fetiSolver = NULL;
810
    fetiData.ksp_schur_primal = PETSC_NULL;
811

Thomas Witkowski's avatar
Thomas Witkowski committed
812 813
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
814
    VecDestroy(&fetiData.tmp_vec_primal);
815 816
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
817 818


819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846
    // === 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;
847 848
    default:
      break;
849
    }
850 851 852 853 854 855 856 857 858
  }


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

859 860
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

861
    // === Get local part of the solution for B variables. ===
862 863 864 865 866

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


867 868
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===