PetscSolverFeti.cc 43.7 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
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createPrimals(feSpace);      
      createDuals(feSpace);
223
      createLagrange(feSpace);      
224 225
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245

    primalDofMap.update();
    dualDofMap.update();
    lagrangeMap.update();
    localDofMap.update();

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

      MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
	  primalDofMap[feSpace].nRankDofs, primalDofMap[feSpace].nOverallDofs);
      
      MSG("nRankDuals = %d   nOverallDuals = %d\n",
	  dualDofMap[feSpace].nRankDofs, 
	  dualDofMap[feSpace].nOverallDofs);

      MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
	  lagrangeMap[feSpace].nRankDofs, 
	  lagrangeMap[feSpace].nOverallDofs);
    }
246 247 248
  }


249
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
250
  {
251
    FUNCNAME("PetscSolverFeti::createPrimals()");  
252

253 254 255
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

256 257 258
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
259
    DofContainerSet& vertices = 
260
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
261 262 263 264
    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);
265

266 267 268 269

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

270
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
271 272
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
273 274
      else
  	primalDofMap[feSpace].insert(*it);
275

Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
276 277 278
    primalDofMap[feSpace].setOverlap(true);
    primalDofMap[feSpace].setDofComm(meshDistributor->getSendDofs(),
  				     meshDistributor->getRecvDofs());
279 280 281
  }


282
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
283 284
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
285

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

    boundaryDofRanks.clear();

291 292 293
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
294
	if (!isPrimal(feSpace, it.getDofIndex())) {
295 296
	  boundaryDofRanks[it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[it.getDofIndex()].insert(it.getRank());
297 298
	}
      }
299
	
300 301 302 303

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

304
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
305 306 307 308

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
309
	if (!isPrimal(feSpace, it.getDofIndex()))
310
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[it.getDofIndex()]);
311 312 313

    stdMpi.updateSendDataSize();

314 315
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
316
      bool recvFromRank = false;
317
      for (; !it.endDofIter(); it.nextDof()) {
318
	if (!isPrimal(feSpace, it.getDofIndex())) {
319 320 321
	  recvFromRank = true;
	  break;
	}
322
      }
323 324

      if (recvFromRank)
325
	stdMpi.recv(it.getRank());
326
    }
327

328 329
    stdMpi.startCommunication();

330 331
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
332
      int i = 0;
333
      for (; !it.endDofIter(); it.nextDof())
334
	if (!isPrimal(feSpace, it.getDofIndex()))
335 336
	  boundaryDofRanks[it.getDofIndex()] = 
	    stdMpi.getRecvData(it.getRank())[i++];
337 338 339 340 341
    }

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

    DofContainer allBoundaryDofs;
342
    meshDistributor->getAllBoundaryDofs(feSpace, allBoundaryDofs);
343 344

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

  
351
  void PetscSolverFeti::createLagrange(const FiniteElemSpace *feSpace)
352 353 354
  {
    FUNCNAME("PetscSolverFeti::createLagrange()");

355 356 357
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

358
    int nRankLagrange = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
359 360
    map<DegreeOfFreedom, MultiIndex>& dualMap = dualDofMap[feSpace].getMap();
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
361
      if (meshDistributor->getIsRankDof(feSpace, it->first)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
362
	lagrangeMap[feSpace].insertRankDof(it->first, nRankLagrange);
363
	int degree = boundaryDofRanks[it->first].size();
364
	nRankLagrange += (degree * (degree - 1)) / 2;
Thomas Witkowski's avatar
Thomas Witkowski committed
365 366
      } else {
	lagrangeMap[feSpace].insert(it->first);
367 368
      }
    }
369
    lagrangeMap[feSpace].nRankDofs = nRankLagrange;
Thomas Witkowski's avatar
Thomas Witkowski committed
370 371 372
    lagrangeMap[feSpace].setOverlap(true);
    lagrangeMap[feSpace].setDofComm(meshDistributor->getSendDofs(),
				    meshDistributor->getRecvDofs());
373 374 375
  }


376
  void PetscSolverFeti::createIndexB(const FiniteElemSpace *feSpace)
377
  {
378
    FUNCNAME("PetscSolverFeti::createIndexB()");
379

380
    DOFAdmin* admin = feSpace->getAdmin();
381 382 383 384

    // === 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.                                 ===
385 386 387

    nLocalInterior = 0;    
    for (int i = 0; i < admin->getUsedSize(); i++) {
388
      if (admin->isDofFree(i) == false && 
389
	  isPrimal(feSpace, i) == false &&
390
	  dualDofMap[feSpace].isSet(i) == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
391
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);
392 393
	nLocalInterior++;
      }
394
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
395
    
Thomas Witkowski's avatar
Thomas Witkowski committed
396
    TEST_EXIT_DBG(nLocalInterior + primalDofMap[feSpace].size() + 
397
		  dualDofMap[feSpace].size() == 
398 399 400
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
403
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualDofMap[feSpace].getMap().begin();
404
	 it != dualDofMap[feSpace].getMap().end(); ++it)
405
      localDofMap[feSpace].insertRankDof(it->first);
406 407 408
  }


409
  void PetscSolverFeti::createMatLagrange(vector<const FiniteElemSpace*> &feSpaces)
410 411 412
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

413 414
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

415 416
    // === Create distributed matrix for Lagrange constraints. ===

417
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
418 419
		    lagrangeMap.getRankDofs(), localDofMap.getRankDofs(),
		    lagrangeMap.getOverallDofs(), localDofMap.getOverallDofs(),	
420 421 422
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

423 424 425 426 427 428 429
    // === 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
430 431
    map<DegreeOfFreedom, MultiIndex> &dualMap = dualDofMap[feSpace].getMap();
    for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); it != dualMap.end(); ++it) {
432 433 434 435
      TEST_EXIT_DBG(boundaryDofRanks.count(it->first))
	("Should not happen!\n");

      // Global index of the first Lagrange constriant for this node.
Thomas Witkowski's avatar
Thomas Witkowski committed
436
      int index = lagrangeMap[feSpace][it->first].local;
437 438 439 440 441 442 443 444 445 446 447 448
      // 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;
Thomas Witkowski's avatar
Thomas Witkowski committed
449 450 451 452
	      //	      int colIndex = it->second * nComponents + k;
	      int colIndex = 
		(localDofMap[feSpace].rStartDofs +
		 localDofMap[feSpace][it->first].local) * nComponents + k;
453
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
454 455
	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
			  INSERT_VALUES);
456 457
	    }
	  }
458 459

	  index++;
460 461 462 463 464 465 466 467 468
	}
      }
    }

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


469
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
470 471 472
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

473 474
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
478 479 480 481 482 483
      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, 
484
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
485 486
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
487
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
488 489 490
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
491 492
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
493 494 495 496 497 498 499 500 501 502
		     &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 {
503 504
      MSG("Create direct schur primal solver!\n");

505 506 507 508
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

509 510
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
511 512
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
513

Thomas Witkowski's avatar
Thomas Witkowski committed
514 515 516 517 518 519
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
520

521 522 523
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
524

525 526
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

527 528
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
529 530 531 532 533 534
	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
535
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
536 537 538 539 540
      }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
541
      TEST_EXIT(mat_b_primal_cols.size() == primalDofMap.getLocalDofs())
542
	("Should not happen!\n");
543 544 545
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
546
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
547 548 549 550 551 552 553 554 555 556

 	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
557
	PetscScalar *tmpValues;
558
	VecGetArray(tmpVec, &tmpValues);
559
	for (int i  = 0; i < nRowsRankB; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
560
	  MatSetValue(matBPi, 
561
		      localDofMap.getStartDofs() + i,
Thomas Witkowski's avatar
Thomas Witkowski committed
562 563 564
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
565 566 567 568 569
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
570 571 572
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
573 574 575
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
576
      MatDestroy(&matBPi);
577 578 579 580 581 582 583 584 585 586

      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
587

588 589
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
590
    }
591 592 593 594 595 596 597
  }


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

598 599 600 601 602
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
603

604 605
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
606

607 608 609 610 611
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
612 613 614
  }


615
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
616 617 618
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

619 620
    // === Create FETI-DP solver object. ===

621 622 623
    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
624
    fetiData.fetiSolver = this;
625
    fetiData.ksp_schur_primal = &ksp_schur_primal;
626

627
    VecCreateMPI(PETSC_COMM_WORLD, 
628
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
629 630
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
631
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
632
		 &(fetiData.tmp_vec_lagrange));
633
    VecCreateMPI(PETSC_COMM_WORLD, 
634
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
635
		 &(fetiData.tmp_vec_primal));
636 637

    MatCreateShell(PETSC_COMM_WORLD,
638 639
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
640
		   &fetiData, &mat_feti);
641 642 643 644 645 646 647
    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);
648 649


650
    // === Create FETI-DP preconditioner object. ===
651

652 653 654 655
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
656

657 658 659
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
660 661
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
662 663 664 665 666 667 668 669 670 671
      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;
      
672
      VecCreateMPI(PETSC_COMM_WORLD, 
673
		   localDofMap.getRankDofs(),localDofMap.getOverallDofs(),
674
		   &(fetiDirichletPreconData.tmp_vec_b));      
675 676 677 678 679 680
      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));
681 682 683 684 685 686 687 688 689 690 691 692
      
      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;

693
      VecCreateMPI(PETSC_COMM_WORLD, 
694 695
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
696
		   &(fetiLumpedPreconData.tmp_vec_b));
697 698 699 700
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
701 702 703 704 705 706

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
707 708
      break;
    default:
709 710
      break;
    }
711 712 713 714 715 716 717
  }
  

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

718 719
    // === Destroy FETI-DP solver object. ===

720 721 722
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
723
    fetiData.fetiSolver = NULL;
724
    fetiData.ksp_schur_primal = PETSC_NULL;
725

Thomas Witkowski's avatar
Thomas Witkowski committed
726 727
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
728
    VecDestroy(&fetiData.tmp_vec_primal);
729 730
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
731 732


733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760
    // === 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;
761 762
    default:
      break;
763
    }
764 765 766 767 768 769 770 771 772
  }


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

773 774
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);

775
    // === Get local part of the solution for B variables. ===
776 777 778 779 780

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


781 782
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
783 784
    
    vector<PetscInt> globalIsIndex, localIsIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
785 786
    globalIsIndex.reserve(primalDofMap.getLocalDofs());
    localIsIndex.reserve(primalDofMap.getLocalDofs());
787 788 789

    {
      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
790
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpace].getMap().begin();
791
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
792
	for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
793
	  globalIsIndex.push_back(it->second.local * nComponents + i);
794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821
	  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);

822 823 824
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
825 826 827 828

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

829
    // === And copy from PETSc local vectors to the DOF vectors. ===
830 831 832 833

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

Thomas Witkowski's avatar
Thomas Witkowski committed
834
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = localDofMap[feSpace].getMap().begin();
835
	   it != localDofMap[feSpace].getMap().end(); ++it) {
836
	int petscIndex = localDofMap.mapLocal(it->first, i);
837 838 839 840
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
841
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpace].getMap().begin();
842
	   it != primalDofMap[feSpace].getMap().end(); ++it) {
843 844 845 846 847 848 849
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
850
    VecDestroy(&local_sol_primal);
851 852 853
  }


854
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
855 856
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
857

858 859
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);

860
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
861
    nComponents = mat->getSize();
862 863 864

    // === Create all sets and indices. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
865 866 867 868
    primalDofMap.init(mpiComm, feSpaces, true);
    dualDofMap.init(mpiComm, feSpaces, false);
    lagrangeMap.init(mpiComm, feSpaces, true);
    localDofMap.init(mpiComm, feSpaces, false);
869

870 871
    updateDofData();

872 873
    // === Create matrices for the FETI-DP method. ===

874 875
    int nRowsRankB = localDofMap.getRankDofs();
    int nRowsOverallB = localDofMap.getOverallDofs();
876 877 878
    int nRowsRankPrimal = primalDofMap.getRankDofs();
    int nRowsOverallPrimal = primalDofMap.getOverallDofs();
    int nRowsDual = dualDofMap.getRankDofs();
879 880
    int nRowsInterior = nRowsRankB - nRowsDual;

Thomas Witkowski's avatar
Thomas Witkowski committed
881 882
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
883 884

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
885 886
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
887
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
888 889

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
890 891
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
892
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
893 894

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
895 896
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
Thomas Witkowski's avatar
Thomas Witkowski committed
897
		    15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
898

899

900 901 902 903
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
904
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
905 906 907 908
		      &mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatCreateSeqAIJ(PETSC_COMM_SELF,
909
		      nRowsInterior, nRowsInterior, 30, PETSC_NULL,
910 911 912
		      &mat_interior_interior);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
913
		      nRowsInterior, nRowsDual, 30, PETSC_NULL,
914 915 916
		      &mat_interior_duals);
      
      MatCreateSeqAIJ(PETSC_COMM_SELF,
917
		      nRowsDual, nRowsInterior, 30, PETSC_NULL,
918 919
		      &mat_duals_interior);
    }
920

921 922
    
    // === Prepare traverse of sequentially created matrices. ===
923 924 925 926 927 928 929 930 931 932 933 934 935 936 937

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits = mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    typedef traits::range_generator<row, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

    vector<int> cols, colsOther;
    vector<double> values, valuesOther;
    cols.reserve(300);
    colsOther.reserve(300);
    values.reserve(300);
    valuesOther.reserve(300);

938 939 940 941 942 943 944
    vector<int> colsLocal, colsLocalOther;
    vector<double> valuesLocal, valuesLocalOther;
    colsLocal.reserve(300);
    colsLocalOther.reserve(300);
    valuesLocal.reserve(300);
    valuesLocalOther.reserve(300);

945 946 947