PetscSolverFeti.cc 43.1 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
      schurPrimalSolver(0)
185 186 187 188 189 190
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

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

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


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

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
217
   
218 219 220 221
    for (unsigned int i = 0; i < meshDistributor->getFeSpaces().size(); i++) {
      const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(i);
      createPrimals(feSpace);      
      createDuals(feSpace);
222
      createLagrange(feSpace);      
223 224
      createIndexB(feSpace);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
225

226 227 228 229 230
    primalDofMap.setDofComm(meshDistributor->getSendDofs(),
			    meshDistributor->getRecvDofs());
    lagrangeMap.setDofComm(meshDistributor->getSendDofs(), 
			   meshDistributor->getRecvDofs());

Thomas Witkowski's avatar
Thomas Witkowski committed
231 232 233 234
    primalDofMap.update();
    dualDofMap.update();
    lagrangeMap.update();
    localDofMap.update();
235 236
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.update();
Thomas Witkowski's avatar
Thomas Witkowski committed
237 238 239 240 241 242 243 244

    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",
245
	  dualDofMap[feSpace].nRankDofs, dualDofMap[feSpace].nOverallDofs);
Thomas Witkowski's avatar
Thomas Witkowski committed
246 247

      MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
248
	  lagrangeMap[feSpace].nRankDofs, lagrangeMap[feSpace].nOverallDofs);
249 250 251 252

      TEST_EXIT_DBG(localDofMap[feSpace].size() + primalDofMap[feSpace].size() == 
		    static_cast<unsigned int>(feSpace->getAdmin()->getUsedDofs()))
	("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
253
    }
254 255 256
  }


257
  void PetscSolverFeti::createPrimals(const FiniteElemSpace *feSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
258
  {
259
    FUNCNAME("PetscSolverFeti::createPrimals()");  
260

261 262 263
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

264 265 266
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
267
    DofContainerSet& vertices = 
268
      meshDistributor->getBoundaryDofInfo(feSpace).geoDofs[VERTEX];
269 270 271 272
    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);
273

274 275 276 277

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

278
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
279 280
      if (meshDistributor->getIsRankDof(feSpace, *it))
	primalDofMap[feSpace].insertRankDof(*it);
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
281 282
      else
  	primalDofMap[feSpace].insert(*it);
283 284 285
  }


286
  void PetscSolverFeti::createDuals(const FiniteElemSpace *feSpace)
287 288
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
289

290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305
    // === Create global index of the dual nodes on each rank. ===

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

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

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

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

309
    boundaryDofRanks[feSpace].clear();
310

311 312 313
    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace); 
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof()) {
314
	if (!isPrimal(feSpace, it.getDofIndex())) {
315 316
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(mpiRank);
	  boundaryDofRanks[feSpace][it.getDofIndex()].insert(it.getRank());
317 318
	}
      }
319
	
320 321 322 323

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

324
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
325 326 327 328

    for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpace);
	 !it.end(); it.nextRank())
      for (; !it.endDofIter(); it.nextDof())
329
	if (!isPrimal(feSpace, it.getDofIndex()))
330
	  stdMpi.getSendData(it.getRank()).push_back(boundaryDofRanks[feSpace][it.getDofIndex()]);
331 332 333

    stdMpi.updateSendDataSize();

334 335
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
336
      bool recvFromRank = false;
337
      for (; !it.endDofIter(); it.nextDof()) {
338
	if (!isPrimal(feSpace, it.getDofIndex())) {
339 340 341
	  recvFromRank = true;
	  break;
	}
342
      }
343 344

      if (recvFromRank)
345
	stdMpi.recv(it.getRank());
346
    }
347

348 349
    stdMpi.startCommunication();

350 351
    for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpace); 
	 !it.end(); it.nextRank()) {
352
      int i = 0;
353
      for (; !it.endDofIter(); it.nextDof())
354
	if (!isPrimal(feSpace, it.getDofIndex()))
355
	  boundaryDofRanks[feSpace][it.getDofIndex()] = 
356
	    stdMpi.getRecvData(it.getRank())[i++];
357 358 359
    }


360 361 362
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

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


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

382
    DOFAdmin* admin = feSpace->getAdmin();
383 384 385 386

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

388
    int nLocalInterior = 0;    
389
    for (int i = 0; i < admin->getUsedSize(); i++) {
390
      if (admin->isDofFree(i) == false && 
391
	  isPrimal(feSpace, i) == false &&
392
	  dualDofMap[feSpace].isSet(i) == false) {
393 394 395 396 397 398
	localDofMap[feSpace].insertRankDof(i, nLocalInterior);

	if (fetiPreconditioner == FETI_DIRICHLET)
	  interiorDofMap[feSpace].insertRankDof(i, nLocalInterior);

	nLocalInterior++;
399
      }
400
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
401
    
402 403
    // === And finally, add the global indicies of all dual nodes. ===

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


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

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

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

422 423 424 425 426 427 428
    // === 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
429 430 431
    for (unsigned int k = 0; k < feSpaces.size(); k++) {
      map<DegreeOfFreedom, MultiIndex> &dualMap = 
	dualDofMap[feSpaces[k]].getMap();
432

Thomas Witkowski's avatar
Thomas Witkowski committed
433 434
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = dualMap.begin(); 
	   it != dualMap.end(); ++it) {
435
	TEST_EXIT_DBG(boundaryDofRanks[feSpaces[k]].count(it->first))
Thomas Witkowski's avatar
Thomas Witkowski committed
436 437 438 439 440 441
	  ("Should not happen!\n");
	
	// Global index of the first Lagrange constriant for this node.
	int index = lagrangeMap.getMatIndex(k, it->first);
	
	// Copy set of all ranks that contain this dual node.
442 443
	vector<int> W(boundaryDofRanks[feSpaces[k]][it->first].begin(), 
		      boundaryDofRanks[feSpaces[k]][it->first].end());
Thomas Witkowski's avatar
Thomas Witkowski committed
444 445 446 447 448 449
	// 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) {
450
	      int colIndex = localDofMap.getMatIndex(k, it->first);
451
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
452
	      MatSetValue(mat_lagrange, index, colIndex, value, INSERT_VALUES);
453
	    }
Thomas Witkowski's avatar
Thomas Witkowski committed
454
	    index++;	      
455 456 457 458 459 460 461 462 463 464
	  }
	}
      }
    }

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


465
  void PetscSolverFeti::createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces)
466 467 468
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
472 473 474 475 476 477
      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, 
478
		   localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
479 480
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
481
		   primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
482 483 484
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
485 486
		     primalDofMap.getRankDofs(), primalDofMap.getRankDofs(), 
		     primalDofMap.getOverallDofs(), primalDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
487 488 489 490 491 492 493 494 495 496
		     &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 {
497 498
      MSG("Create direct schur primal solver!\n");

499 500 501 502
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

503 504
      int nRowsRankPrimal = primalDofMap.getRankDofs();
      int nRowsOverallPrimal = primalDofMap.getOverallDofs();
505 506
      int nRowsRankB = localDofMap.getRankDofs();
      int nRowsOverallB = localDofMap.getOverallDofs();
507

Thomas Witkowski's avatar
Thomas Witkowski committed
508 509 510 511 512 513
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
514

515 516 517
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
518

519 520
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

521 522
      for (int i = 0; i < nRowsRankB; i++) {
	PetscInt row = localDofMap.getStartDofs() + i;
523 524 525 526 527 528
	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
529
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
530 531 532 533 534
      }

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

535 536
      TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) == 
		primalDofMap.getLocalDofs())
537
	("Should not happen!\n");
538

539 540 541
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
542
	VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
543 544 545 546 547 548 549 550 551 552

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

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
566 567 568
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
569 570 571
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
572
      MatDestroy(&matBPi);
573 574 575 576 577 578 579 580 581 582

      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
583

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


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

594 595 596 597 598
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
599

600 601
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
602

603 604 605 606 607
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
608 609 610
  }


611
  void PetscSolverFeti::createFetiKsp(vector<const FiniteElemSpace*> &feSpaces)
612 613 614
  {
    FUNCNAME("PetscSolverFeti::createFetiKsp()");

615 616
    // === Create FETI-DP solver object. ===

617 618 619
    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
620
    fetiData.fetiSolver = this;
621
    fetiData.ksp_schur_primal = &ksp_schur_primal;
622

623
    VecCreateMPI(PETSC_COMM_WORLD, 
624
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
625 626
		 &(fetiData.tmp_vec_b));
    VecCreateMPI(PETSC_COMM_WORLD,
627
		 lagrangeMap.getRankDofs(), lagrangeMap.getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
628
		 &(fetiData.tmp_vec_lagrange));
629
    VecCreateMPI(PETSC_COMM_WORLD, 
630
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
631
		 &(fetiData.tmp_vec_primal));
632 633

    MatCreateShell(PETSC_COMM_WORLD,
634 635
		   lagrangeMap.getRankDofs(), lagrangeMap.getRankDofs(),
		   lagrangeMap.getOverallDofs(), lagrangeMap.getOverallDofs(),
636
		   &fetiData, &mat_feti);
637 638 639 640 641 642 643
    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);
644 645


646
    // === Create FETI-DP preconditioner object. ===
647

648 649 650 651
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
652

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

689
      VecCreateMPI(PETSC_COMM_WORLD, 
690 691
		   localDofMap.getRankDofs(),
		   localDofMap.getOverallDofs(),
692
		   &(fetiLumpedPreconData.tmp_vec_b));
693 694 695 696
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, 
		 &(fetiLumpedPreconData.tmp_vec_duals1));
697 698 699 700 701 702

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
703 704
      break;
    default:
705 706
      break;
    }
707 708 709 710 711 712 713
  }
  

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

714 715
    // === Destroy FETI-DP solver object. ===

716 717 718
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
719
    fetiData.fetiSolver = NULL;
720
    fetiData.ksp_schur_primal = PETSC_NULL;
721

Thomas Witkowski's avatar
Thomas Witkowski committed
722 723
    VecDestroy(&fetiData.tmp_vec_b);
    VecDestroy(&fetiData.tmp_vec_lagrange);
724
    VecDestroy(&fetiData.tmp_vec_primal);
725 726
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
727 728


729 730 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
    // === 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;
757 758
    default:
      break;
759
    }
760 761 762 763 764 765 766 767 768
  }


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

769
    // === Get local part of the solution for B variables. ===
770 771 772 773 774

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


775 776
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
777
    
Thomas Witkowski's avatar
Thomas Witkowski committed
778 779
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(&vec);

780
    vector<PetscInt> globalIsIndex, localIsIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
781 782
    globalIsIndex.reserve(primalDofMap.getLocalDofs());
    localIsIndex.reserve(primalDofMap.getLocalDofs());
783 784

    {
Thomas Witkowski's avatar
Thomas Witkowski committed
785 786 787 788 789 790 791 792
      int cnt = 0;
      for (unsigned int i = 0; i < feSpaces.size(); i++) {
	map<DegreeOfFreedom, MultiIndex>& dofMap = 
	  primalDofMap[feSpaces[i]].getMap();
	for (map<DegreeOfFreedom, MultiIndex>::iterator it = dofMap.begin();
	     it != dofMap.end(); ++it) {
	  globalIsIndex.push_back(primalDofMap.getMatIndex(i, it->first));
	  localIsIndex.push_back(cnt++);
793 794
	}
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
795 796

      TEST_EXIT_DBG(cnt == primalDofMap.getLocalDofs())("Should not happen!\n");
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
    }
    
    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

Thomas Witkowski's avatar
Thomas Witkowski committed
831 832
    int cnt = 0;
    for (int i = 0; i < vec.getSize(); i++) {
833 834
      DOFVector<double>& dofVec = *(vec.getDOFVector(i));

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

Thomas Witkowski's avatar
Thomas Witkowski committed
841 842 843
      for (map<DegreeOfFreedom, MultiIndex>::iterator it = primalDofMap[feSpaces[i]].getMap().begin();
	   it != primalDofMap[feSpaces[i]].getMap().end(); ++it)
	dofVec[it->first] = localSolPrimal[cnt++];
844 845 846 847
    }

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
848
    VecDestroy(&local_sol_primal);
849 850 851
  }


852
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
853 854
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
855
    
856
    // === Create all sets and indices. ===
857
    
Thomas Witkowski's avatar
Thomas Witkowski committed
858
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
859
    
860 861 862 863
    primalDofMap.init(mpiComm, feSpaces, true, true);
    dualDofMap.init(mpiComm, feSpaces, false, false);
    lagrangeMap.init(mpiComm, feSpaces, true, true);
    localDofMap.init(mpiComm, feSpaces, false, false);
864 865
    if (fetiPreconditioner == FETI_DIRICHLET)
      interiorDofMap.init(mpiComm, feSpaces, false, false);
866

867
    updateDofData();
868
    
869
    // === Create matrices for the FETI-DP method. ===
870
    
871 872
    int nRowsRankB = localDofMap.getRankDofs();
    int nRowsOverallB = localDofMap.getOverallDofs();
873 874
    int nRowsRankPrimal = primalDofMap.getRankDofs();
    int nRowsOverallPrimal = primalDofMap.getOverallDofs();
875
    
Thomas Witkowski's avatar
Thomas Witkowski committed
876 877
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
878
    
879
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
880 881
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
882
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
883
    
884
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
885 886
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
887
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
888
    
889
    MatCreateMPIAIJ(PETSC_COMM_WORLD,