PetscSolverFeti.cc 49 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
//
// 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"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
16
#include "io/VtkWriter.h"
17 18 19 20 21

namespace AMDiS {

  using namespace std;

22
#ifdef HAVE_PETSC_DEV 
Thomas Witkowski's avatar
Thomas Witkowski committed
23
  FetiStatisticsData PetscSolverFeti::fetiStatistics;
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

Thomas Witkowski's avatar
Thomas Witkowski committed
30
    double first = MPI::Wtime();
31 32
    void *ctx;
    MatShellGetContext(mat, &ctx);
33
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
34 35

    MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
Thomas Witkowski's avatar
Thomas Witkowski committed
36
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
37 38 39 40
    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);

Thomas Witkowski's avatar
Thomas Witkowski committed
41 42 43
    PetscSolverFeti::fetiStatistics.nSchurPrimalApply++;
    PetscSolverFeti::fetiStatistics.timeSchurPrimalApply += MPI::Wtime() - first;

44 45 46 47 48 49 50
    return 0;
  }


  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
51 52
    //    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)
53 54 55

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

58 59
    // tmp_vec_b0 = inv(K_BB) trans(L) x
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b0);
Thomas Witkowski's avatar
Thomas Witkowski committed
60
    data->fetiSolver->solveLocalProblem(data->tmp_vec_b0, data->tmp_vec_b0);
61

62 63
    // tmp_vec_b1 = inv(K_BB) K_BPi  inv(S_PiPi) K_PiB tmp_vec_b0
    MatMult(*(data->mat_primal_b), data->tmp_vec_b0, data->tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
64 65

    double first = MPI::Wtime();
66
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
Thomas Witkowski's avatar
Thomas Witkowski committed
67 68 69 70
    PetscSolverFeti::fetiStatistics.nSchurPrimalSolve++;
    PetscSolverFeti::fetiStatistics.timeSchurPrimalSolve +=
      MPI::Wtime() - first;

71
    MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b1);
72

73 74 75
    // y = L (tmp_vec_b0 + tmp_vec_b1)
    VecAXPBY(data->tmp_vec_b0, 1.0, 1.0, data->tmp_vec_b1);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b0, y);
76

Thomas Witkowski's avatar
Thomas Witkowski committed
77 78
    PetscSolverFeti::fetiStatistics.nFetiApply++;

79 80 81 82
    return 0;
  }


83
  // y = PC * x
84
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
85
  {
86
    // Get data for the preconditioner
87 88
    void *ctx;
    PCShellGetContext(pc, &ctx);
89
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
90

91
    // Multiply with scaled Lagrange constraint matrix.
92 93 94
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


95
    // === Restriction of the B nodes to the boundary nodes. ===
96

97 98 99 100
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
101

102 103 104 105 106 107
    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];
108 109

    VecRestoreArray(data->tmp_vec_b, &local_b);
110
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
111 112


113
    // === K_DD - K_DI inv(K_II) K_ID ===
114

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

117
    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
118
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152
    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);
153 154


155
    // === Restriction of the B nodes to the boundary nodes. ===
156

157 158 159 160
    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
161

162
    PetscScalar *local_b, *local_duals;
163
    VecGetArray(data->tmp_vec_b, &local_b);
164
    VecGetArray(data->tmp_vec_duals0, &local_duals);
165

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

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


    // === K_DD ===

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

177

178
    // === Prolongation from local dual nodes to B nodes.
179

180 181 182 183 184
    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];
185

186 187
    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
188

189 190

    // Multiply with scaled Lagrange constraint matrix.
191 192 193 194 195 196
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    return 0;
  }


197 198
  PetscSolverFeti::PetscSolverFeti()
    : PetscSolver(),
Thomas Witkowski's avatar
Thomas Witkowski committed
199 200
      nComponents(-1),
      schurPrimalSolver(0)
201 202 203 204 205 206
  {
    FUNCNAME("PetscSolverFeti::PetscSolverFeti()");

    string preconditionerName = "";
    Parameters::get("parallel->solver->precon", preconditionerName);
    if (preconditionerName == "" || preconditionerName == "none") {
207
      MSG("Create FETI-DP solver with no preconditioner!\n");
208 209
      fetiPreconditioner = FETI_NONE;
    } else if (preconditionerName == "dirichlet") {
210
      MSG("Create FETI-DP solver with Dirichlet preconditioner!\n");
211 212
      fetiPreconditioner = FETI_DIRICHLET;
    } else if (preconditionerName == "lumped") {
213
      MSG("Create FETI-DP solver with lumped preconditioner!\n");
214 215
      fetiPreconditioner = FETI_LUMPED;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
216 217
      ERROR_EXIT("Preconditioner \"%s\" not available!\n", 
		 preconditionerName.c_str());
218
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
219 220 221 222 223

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


227
  void PetscSolverFeti::updateDofData()
228 229
  {
    FUNCNAME("PetscSolverFeti::updateDofData()");
230 231 232

    TEST_EXIT(meshDistributor->getFeSpace()->getBasisFcts()->getDegree() == 1)
      ("Works for linear basis functions only!\n");
233
   
234 235 236 237 238 239 240
    createPrimals();

    createDuals();

    createLagrange();

    createIndexB();
241 242 243
  }


244
  void PetscSolverFeti::createPrimals()
245
  {
246
    FUNCNAME("PetscSolverFeti::createPrimals()");  
247

248 249 250
    // === Define all vertices on the interior boundaries of the macro mesh ===
    // === to be primal variables.                                          ===

251 252 253
    /// Set of DOF indices that are considered to be primal variables.
    
    DofIndexSet primals;
254 255 256 257 258 259
    DofContainerSet& vertices = 
      meshDistributor->getBoundaryDofInfo().geoDofs[VERTEX];
    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);
260

261 262 263 264

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

265
    globalPrimalIndex.clear();
266 267 268 269
    nRankPrimals = 0;
    for (DofIndexSet::iterator it = primals.begin(); it != primals.end(); ++it)
      if (meshDistributor->getIsRankDof(*it)) {
	globalPrimalIndex[*it] = nRankPrimals;
270 271 272
	nRankPrimals++;
      }

273

274 275 276
    // === Get overall number of primals and rank's displacement in the ===
    // === numbering of the primals.                                    ===

277
    nOverallPrimals = 0;
278
    rStartPrimals = 0;
279 280 281
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankPrimals, rStartPrimals, nOverallPrimals);

282 283 284

    // === Create global primal index for all primals. ===

285 286 287 288
    for (DofMapping::iterator it = globalPrimalIndex.begin();
	 it != globalPrimalIndex.end(); ++it)
      it->second += rStartPrimals;

289 290
    MSG("nRankPrimals = %d   nOverallPrimals = %d\n", 
	nRankPrimals, nOverallPrimals);
291

292 293 294 295 296

    // === Communicate primal's global index from ranks that own the     ===
    // === primals to ranks that contain this primals but are not owning ===
    // === them.                                                         ===

297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344
    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (globalPrimalIndex.count(**dofIt))
	  stdMpi.getSendData(it->first).push_back(globalPrimalIndex[**dofIt]);
    stdMpi.updateSendDataSize();

    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      bool recvFromRank = false;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
	if (primals.count(**dofIt) && 
	    meshDistributor->getIsRankDof(**dofIt) == false) {
	  recvFromRank = true;
	  break;
	}

      if (recvFromRank) 
	stdMpi.recv(it->first);
    }
    stdMpi.startCommunication();

    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      int i = 0;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	if (primals.count(**dofIt) && 
	    meshDistributor->getIsRankDof(**dofIt) == false)
	  globalPrimalIndex[**dofIt] = stdMpi.getRecvData(it->first)[i++];
      }
    }

    TEST_EXIT_DBG(primals.size() == globalPrimalIndex.size())
      ("Number of primals %d, but number of global primals on this rank is %d!\n",
       primals.size(), globalPrimalIndex.size());
  }


  void PetscSolverFeti::createDuals()
  {
    FUNCNAME("PetscSolverFeti::createDuals()");
    
345 346
    // === Create for each dual node that is owned by the rank, the set ===
    // === of ranks that contain this node (denoted by W(x_j)).         ===
347 348 349 350 351 352 353 354 355

    boundaryDofRanks.clear();

    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it) {
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
	// If DOF is not primal, i.e., its a dual node
356
	if (globalPrimalIndex.count(**dofIt) == 0) {
357 358 359 360 361 362
	  boundaryDofRanks[**dofIt].insert(mpiRank);
	  boundaryDofRanks[**dofIt].insert(it->first);
	}
      }
    }

363 364 365 366

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

367 368 369 370 371
    StdMpi<vector<std::set<int> > > stdMpi(meshDistributor->getMpiComm());
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
372
	if (globalPrimalIndex.count(**dofIt) == 0)
373 374 375 376 377 378 379 380 381 382
	  stdMpi.getSendData(it->first).push_back(boundaryDofRanks[**dofIt]);

    stdMpi.updateSendDataSize();

    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      bool recvFromRank = false;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
383
	if (globalPrimalIndex.count(**dofIt) == 0) {
384 385 386 387 388 389 390 391 392 393 394 395 396 397
	  recvFromRank = true;
	  break;
	}

      if (recvFromRank)
	stdMpi.recv(it->first);
    }
    stdMpi.startCommunication();

    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      int i = 0;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)	
398
	if (globalPrimalIndex.count(**dofIt) == 0)
399 400 401 402 403 404 405 406 407 408
	  boundaryDofRanks[**dofIt] = stdMpi.getRecvData(it->first)[i++];	      
    }


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

    duals.clear();
    globalDualIndex.clear();

    int nRankAllDofs = meshDistributor->getFeSpace()->getAdmin()->getUsedDofs();
409
    nRankB = nRankAllDofs - globalPrimalIndex.size();
410 411 412 413 414 415 416 417 418 419 420
    nOverallB = 0;
    rStartB = 0;
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankB, rStartB, nOverallB);
    DofContainer allBoundaryDofs;
    meshDistributor->getAllBoundaryDofs(allBoundaryDofs);
    int nRankInteriorDofs = nRankAllDofs - allBoundaryDofs.size();

    int nRankDuals = 0;
    for (DofContainer::iterator it = allBoundaryDofs.begin();
	 it != allBoundaryDofs.end(); ++it) {
421
      if (globalPrimalIndex.count(**it) == 0) {
422 423 424 425 426 427 428 429 430
	duals.insert(**it);
	globalDualIndex[**it] = rStartB + nRankInteriorDofs + nRankDuals;
	nRankDuals++;
      }
    }

    int nOverallDuals = nRankDuals;
    mpi::globalAdd(nOverallDuals);

431 432
    MSG("nRankDuals = %d   nOverallDuals = %d\n",
	nRankDuals, nOverallDuals);
433 434 435 436 437 438 439
  }

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

440 441 442
    // === Reserve for each dual node, on the rank that owns this node, the ===
    // === appropriate number of Lagrange constraints.                      ===

443 444 445 446 447 448 449 450 451
    nRankLagrange = 0;
    for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) {
      if (meshDistributor->getIsRankDof(*it)) {
	dofFirstLagrange[*it] = nRankLagrange;
	int degree = boundaryDofRanks[*it].size();
	nRankLagrange += (degree * (degree - 1)) / 2;
      }
    }

452 453 454 455 456

    // === Get the overall number of Lagrange constraints and create the ===
    // === mapping dofFirstLagrange, that defines for each dual boundary ===
    // === node the first Lagrange constraint global index.              ===

457
    nOverallLagrange = 0;
458
    rStartLagrange = 0;
459 460 461 462 463 464 465
    mpi::getDofNumbering(meshDistributor->getMpiComm(),
			 nRankLagrange, rStartLagrange, nOverallLagrange);

    for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it)
      if (meshDistributor->getIsRankDof(*it))
	dofFirstLagrange[*it] += rStartLagrange;

466 467
    MSG("nRankLagrange = %d  nOverallLagrange = %d\n",
	nRankLagrange, nOverallLagrange);
468 469


470
    // === Communicate dofFirstLagrange to all other ranks. ===
471 472 473 474 475 476 477

    StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
    RankToDofContainer& sendDofs = meshDistributor->getSendDofs();
    for (RankToDofContainer::iterator it = sendDofs.begin();
	 it != sendDofs.end(); ++it)
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt) {
478
	if (globalPrimalIndex.count(**dofIt) == 0) {
479 480 481 482 483 484 485 486 487 488 489 490
	  TEST_EXIT_DBG(dofFirstLagrange.count(**dofIt))("Should not happen!\n");
	  stdMpi.getSendData(it->first).push_back(dofFirstLagrange[**dofIt]);
	}
      }
    stdMpi.updateSendDataSize();

    RankToDofContainer& recvDofs = meshDistributor->getRecvDofs();
    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      bool recvData = false;
      for (DofContainer::iterator dofIt = it->second.begin();
	   dofIt != it->second.end(); ++dofIt)
491
	if (globalPrimalIndex.count(**dofIt) == 0) {
492 493 494 495 496 497 498 499 500 501 502 503 504 505
	  recvData = true;
	  break;
	}
	  
      if (recvData)
	stdMpi.recv(it->first);
    }

    stdMpi.startCommunication();

    for (RankToDofContainer::iterator it = recvDofs.begin();
	 it != recvDofs.end(); ++it) {
      int counter = 0;
      for (unsigned int i = 0; i < it->second.size(); i++)
506
	if (globalPrimalIndex.count(*(it->second[i])) == 0)
507
	  dofFirstLagrange[*(it->second[i])] = stdMpi.getRecvData(it->first)[counter++];
508
    }     
509 510 511 512 513 514 515
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
516
    localIndexB.clear();
517
    DOFAdmin* admin = meshDistributor->getFeSpace()->getAdmin();
518 519 520 521

    // === 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.                                 ===
522 523
    
    for (int i = 0; i < admin->getUsedSize(); i++)
524 525
      if (admin->isDofFree(i) == false && globalPrimalIndex.count(i) == 0)
	if (duals.count(i) == 0 && globalPrimalIndex.count(i) == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
526
	  localIndexB[i] = -1;
527

528 529 530

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

531
    nLocalInterior = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
532 533 534
    for (DofMapping::iterator it = localIndexB.begin(); 
	 it != localIndexB.end(); ++it) {
      it->second = nLocalInterior;
535
      nLocalInterior++;
536
    }
537
    nLocalDuals = duals.size();
538

539
    TEST_EXIT_DBG(nLocalInterior + globalPrimalIndex.size() + duals.size() == 
540 541 542
		  static_cast<unsigned int>(admin->getUsedDofs()))
      ("Should not happen!\n");

543 544 545

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

546 547
    for (DofIndexSet::iterator it = duals.begin();
	 it != duals.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
548
      localIndexB[*it] = globalDualIndex[*it] - rStartB;
549 550 551
  }


552
  void PetscSolverFeti::createMatLagrange()
553 554 555
  {
    FUNCNAME("PetscSolverFeti::createMatLagrange()");

556 557
    // === Create distributed matrix for Lagrange constraints. ===

558
    MatCreateMPIAIJ(PETSC_COMM_WORLD,
559 560 561 562
		    nRankLagrange * nComponents, 
		    nRankB * nComponents,
		    nOverallLagrange * nComponents, 
		    nOverallB * nComponents,
563 564 565
		    2, PETSC_NULL, 2, PETSC_NULL,
		    &mat_lagrange);

566 567 568 569 570 571 572
    // === 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.                                                     ===

573 574 575 576
    for (DofIndexSet::iterator it = duals.begin(); it != duals.end(); ++it) {
      TEST_EXIT_DBG(dofFirstLagrange.count(*it))("Should not happen!\n");
      TEST_EXIT_DBG(boundaryDofRanks.count(*it))("Should not happen!\n");

577
      // Global index of the first Lagrange constriant for this node.
578
      int index = dofFirstLagrange[*it];
579
      // Copy set of all ranks that contain this dual node.
580
      vector<int> W(boundaryDofRanks[*it].begin(), boundaryDofRanks[*it].end());
581
      // Number of ranks that contain this dual node.
582 583 584 585 586 587 588
      int degree = W.size();

      TEST_EXIT_DBG(globalDualIndex.count(*it))("Should not happen!\n");
      int dualCol = globalDualIndex[*it];

      for (int i = 0; i < degree; i++) {
	for (int j = i + 1; j < degree; j++) {
589 590
	  if (W[i] == mpiRank || W[j] == mpiRank) {
	    // Set the constraint for all components of the system.
591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609
	    for (int k = 0; k < nComponents; k++) {
	      int rowIndex = index * nComponents + k;
	      int colIndex = dualCol * nComponents + k;
	      double value = (W[i] == mpiRank ? 1.0 : -1.0);
	      MatSetValue(mat_lagrange, rowIndex, colIndex, value, 
			  INSERT_VALUES);
	    }
	  }

	  index++;
	}
      }
    }

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


610
  void PetscSolverFeti::createSchurPrimalKsp()
611 612 613
  {
    FUNCNAME("PetscSolverFeti::createSchurPrimal()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641
      schurPrimalData.mat_primal_primal = &mat_primal_primal;
      schurPrimalData.mat_primal_b = &mat_primal_b;
      schurPrimalData.mat_b_primal = &mat_b_primal;
      schurPrimalData.fetiSolver = this;
      
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(schurPrimalData.tmp_vec_b));
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankPrimals * nComponents, nOverallPrimals * nComponents,
		   &(schurPrimalData.tmp_vec_primal));

      MatCreateShell(PETSC_COMM_WORLD,
		     nRankPrimals * nComponents, nRankPrimals * nComponents,
		     nOverallPrimals * nComponents, nOverallPrimals * nComponents,
		     &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 {
642 643
      MSG("Create direct schur primal solver!\n");

644 645 646 647 648 649
      TEST_EXIT(ksp_b)("No solver object for local problem!\n");

      double wtime = MPI::Wtime();

      int nRowsRankPrimal = nRankPrimals * nComponents;
      int nRowsOverallPrimal = nOverallPrimals * nComponents;
Thomas Witkowski's avatar
Thomas Witkowski committed
650 651
      int nRowsRankB = nRankB * nComponents;
      int nRowsOverallB = nOverallB * nComponents;
652

Thomas Witkowski's avatar
Thomas Witkowski committed
653 654 655 656 657 658
      Mat matBPi;
      MatCreateMPIAIJ(PETSC_COMM_WORLD,
		      nRowsRankB, nRowsRankPrimal, 
		      nRowsOverallB, nRowsOverallPrimal,
		      30, PETSC_NULL, 30, PETSC_NULL, &matBPi);
      Mat matPrimal;
659

660 661 662
      PetscInt nCols;
      const PetscInt *cols;
      const PetscScalar *values;
Thomas Witkowski's avatar
Thomas Witkowski committed
663

664 665 666 667 668 669 670 671 672 673 674
      int nLocalPrimals = globalPrimalIndex.size() * nComponents;
      map<int, vector<pair<int, double> > > mat_b_primal_cols;

      for (int i = 0; i < nRankB * nComponents; i++) {
	PetscInt row = rStartB * nComponents + i;
	MatGetRow(mat_b_primal, row, &nCols, &cols, &values);

	for (int j = 0; j < nCols; j++)
	  if (values[j] != 0.0)
	    mat_b_primal_cols[cols[j]].push_back(make_pair(i, values[j]));

Thomas Witkowski's avatar
Thomas Witkowski committed
675
	MatRestoreRow(mat_b_primal, row, &nCols, &cols, &values);
676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696
      }

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

      TEST_EXIT(mat_b_primal_cols.size() ==
		globalPrimalIndex.size() * nComponents)("Should not happen!\n");
      for (map<int, vector<pair<int, double> > >::iterator it = mat_b_primal_cols.begin();
	   it != mat_b_primal_cols.end(); ++it) {
	Vec tmpVec;
	VecCreateSeq(PETSC_COMM_SELF, nRankB * nComponents, &tmpVec);

 	for (unsigned int i = 0; i < it->second.size(); i++) 
 	  VecSetValue(tmpVec, 
 		      it->second[i].first, it->second[i].second, INSERT_VALUES);

	VecAssemblyBegin(tmpVec);
	VecAssemblyEnd(tmpVec);

       	KSPSolve(ksp_b, tmpVec, tmpVec);

Thomas Witkowski's avatar
Thomas Witkowski committed
697
	PetscScalar *tmpValues;
698 699
	VecGetArray(tmpVec, &tmpValues);
	for (int i  = 0; i < nRankB * nComponents; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
700 701 702 703 704
	  MatSetValue(matBPi, 
		      rStartB * nComponents + i,
		      it->first,
		      tmpValues[i],
		      ADD_VALUES);
705 706 707 708 709
	VecRestoreArray(tmpVec, &tmpValues);

	VecDestroy(&tmpVec);
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
710 711 712
      MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
      MatMatMult(mat_primal_b, matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
713 714 715
      MatAXPY(mat_primal_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);

      MatDestroy(&matPrimal);
Thomas Witkowski's avatar
Thomas Witkowski committed
716
      MatDestroy(&matBPi);
717 718 719 720 721 722 723 724 725 726

      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
727

728 729
      MSG("Creating Schur primal matrix needed %.5f seconds.\n",
	  MPI::Wtime() - wtime);
Thomas Witkowski's avatar
Thomas Witkowski committed
730
    }
731 732 733 734 735 736 737
  }


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

738 739 740 741 742
    if (schurPrimalSolver == 0) {
      schurPrimalData.mat_primal_primal = PETSC_NULL;
      schurPrimalData.mat_primal_b = PETSC_NULL;
      schurPrimalData.mat_b_primal = PETSC_NULL;
      schurPrimalData.fetiSolver = NULL;
743

744 745
      VecDestroy(&schurPrimalData.tmp_vec_b);
      VecDestroy(&schurPrimalData.tmp_vec_primal);
746

747 748 749 750 751
      MatDestroy(&mat_schur_primal);
      KSPDestroy(&ksp_schur_primal);
    } else {
      KSPDestroy(&ksp_schur_primal);
    }
752 753 754 755 756 757 758
  }


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

759 760
    // === Create FETI-DP solver object. ===

761 762 763
    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
764
    fetiData.fetiSolver = this;
765
    fetiData.ksp_schur_primal = &ksp_schur_primal;
766

767 768 769 770 771 772 773 774 775
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents,
		 &(fetiData.tmp_vec_b0));
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents,
		 &(fetiData.tmp_vec_b1));
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankPrimals * nComponents, nOverallPrimals * nComponents,
		 &(fetiData.tmp_vec_primal));
776 777

    MatCreateShell(PETSC_COMM_WORLD,
778 779
		   nRankLagrange * nComponents, nRankLagrange * nComponents,
		   nOverallLagrange * nComponents, nOverallLagrange * nComponents,
780
		   &fetiData, &mat_feti);
781 782 783 784 785 786 787
    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);
788 789


790
    // === Create FETI-DP preconditioner object. ===
791

792 793 794 795
    if (fetiPreconditioner != FETI_NONE) {
      MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
      MatScale(mat_lagrange_scaled, 0.5);
    }
796

797 798 799
    switch (fetiPreconditioner) {
    case FETI_DIRICHLET:           
      KSPCreate(PETSC_COMM_SELF, &ksp_interior);
Thomas Witkowski's avatar
Thomas Witkowski committed
800 801
      KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, 
		      SAME_NONZERO_PATTERN);
802 803 804 805 806 807 808 809 810 811
      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;
      
812 813 814
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiDirichletPreconData.tmp_vec_b));      
815 816 817 818 819 820 821 822 823 824 825 826 827 828 829
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_duals1));
      MatGetVecs(mat_interior_interior, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_interior));
      
      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiDirichletPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiDirichletPrecon);
      
      break;

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

830 831 832
      VecCreateMPI(PETSC_COMM_WORLD, 
		   nRankB * nComponents, nOverallB * nComponents,
		   &(fetiLumpedPreconData.tmp_vec_b));
833 834 835 836 837 838 839 840
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals0));
      MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals1));

      KSPGetPC(ksp_feti, &precon_feti);
      PCSetType(precon_feti, PCSHELL);
      PCShellSetContext(precon_feti, static_cast<void*>(&fetiLumpedPreconData));
      PCShellSetApply(precon_feti, petscApplyFetiLumpedPrecon);
      
841 842
      break;
    default:
843 844
      break;
    }
845 846 847 848 849 850 851
  }
  

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

852 853
    // === Destroy FETI-DP solver object. ===

854 855 856
    fetiData.mat_primal_b = PETSC_NULL;
    fetiData.mat_b_primal = PETSC_NULL;
    fetiData.mat_lagrange = PETSC_NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
857
    fetiData.fetiSolver = NULL;
858
    fetiData.ksp_schur_primal = PETSC_NULL;
859

860 861
    VecDestroy(&fetiData.tmp_vec_b0);
    VecDestroy(&fetiData.tmp_vec_b1);
862
    VecDestroy(&fetiData.tmp_vec_primal);
863 864
    MatDestroy(&mat_feti);
    KSPDestroy(&ksp_feti);
865 866


867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894
    // === 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;
895 896
    default:
      break;
897
    }
898 899 900 901 902 903 904 905 906
  }


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

907 908 909 910 911 912
    PetscScalar *ttt;
    VecGetArray(vec_sol_primal, &ttt);
    for (int i = 0; i < nRankPrimals; i++)
      MSG("PRIM VAL = %f\n", ttt[i * 3 + 1]);
    VecRestoreArray(vec_sol_primal, &ttt);

913
    // === Get local part of the solution for B variables. ===
914 915 916 917 918

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


919 920
    // === Create scatter to get solutions of all primal nodes that are ===
    // === contained in rank's domain.                                  ===
921 922 923 924 925 926 927 928 929 930
    
    vector<PetscInt> globalIsIndex, localIsIndex;
    globalIsIndex.reserve(globalPrimalIndex.size() * nComponents);
    localIsIndex.reserve(globalPrimalIndex.size() * nComponents);

    {
      int counter = 0;
      for (DofMapping::iterator it = globalPrimalIndex.begin();
	   it != globalPrimalIndex.end(); ++it) {
	for (int i = 0; i < nComponents; i++) {
931 932
	  MSG("COPY GLOBAL PRIM %d of COMP %d TO LOCAL TMP %d\n",
	      it->second, i, counter);
933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961
	  globalIsIndex.push_back(it->second * nComponents + i);
	  localIsIndex.push_back(counter++);
	}
      }
    }
    
    IS globalIs, localIs;
    ISCreateGeneral(PETSC_COMM_SELF, 
		    globalIsIndex.size(), 
		    &(globalIsIndex[0]),
		    PETSC_USE_POINTER,
		    &globalIs);

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

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

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

962 963 964
    ISDestroy(&globalIs);
    ISDestroy(&localIs);    
    VecScatterDestroy(&primalScatter);    
965 966 967 968

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

969 970 971
    for (int i = 0; i < globalPrimalIndex.size(); i++)
      MSG("TEST VAL %f\n", localSolPrimal[i * 3 + 1]);

972

973
    // === And copy from PETSc local vectors to the DOF vectors. ===
974 975 976 977

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

Thomas Witkowski's avatar
Thomas Witkowski committed
978 979 980
      for (DofMapping::iterator it = localIndexB.begin();
	   it != localIndexB.end(); ++it) {
	int petscIndex = it->second * nComponents + i;
981 982 983 984 985 986
	dofVec[it->first] = localSolB[petscIndex];
      }

      int counter = 0;
      for (DofMapping::iterator it = globalPrimalIndex.begin();
	   it != globalPrimalIndex.end(); ++it) {
987 988 989 990
	if (i == 1)
	  MSG("AND SET THE VAL OF DOF %d TO VAL %f\n", 
	      it->first,
	      localSolPrimal[counter * nComponents + i]);
991 992 993 994 995
	dofVec[it->first] = localSolPrimal[counter * nComponents + i];
	counter++;
      }
    }

996
    MSG("VAL BOUNDS: %f %f\n", vec.getDOFVector(1)->min(), vec.getDOFVector(1)->max());
997 998 999

    VecRestoreArray(vec_sol_b, &localSolB);
    VecRestoreArray(local_sol_primal, &localSolPrimal);
1000
    VecDestroy(&local_sol_primal);
1001 1002 1003
  }


1004
  void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
1005 1006
  {
    FUNCNAME("PetscSolverFeti::fillPetscMatrix()");   
1007

1008
    nComponents = mat->getSize();
1009 1010 1011

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

1012 1013
    updateDofData();

1014 1015 1016 1017 1018 1019 1020

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

    int nRowsRankB = nRankB * nComponents;
    int nRowsOverallB = nOverallB * nComponents;
    int nRowsRankPrimal = nRankPrimals * nComponents;
    int nRowsOverallPrimal = nOverallPrimals * nComponents;
1021
    int nRowsInterior = nLocalInterior * nComponents;
1022
    int nRowsDual = nLocalDuals * nComponents;
1023

Thomas Witkowski's avatar
Thomas Witkowski committed
1024 1025
    MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 30, PETSC_NULL,
		    &mat_b_b);
1026 1027

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
1028 1029
		    nRowsRankPrimal, nRowsRankPrimal, 
		    nRowsOverallPrimal, nRowsOverallPrimal,
1030
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_primal);
1031 1032

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
1033 1034
		    nRowsRankB, nRowsRankPrimal, 
		    nRowsOverallB, nRowsOverallPrimal,
1035
		    30, PETSC_NULL, 30, PETSC_NULL, &mat_b_primal);
1036 1037

    MatCreateMPIAIJ(PETSC_COMM_WORLD,
1038 1039
		    nRowsRankPrimal, nRowsRankB,
		    nRowsOverallPrimal, nRowsOverallB,
Thomas Witkowski's avatar
Thomas Witkowski committed
1040
		    15, PETSC_NULL, 15, PETSC_NULL, &mat_primal_b);
1041

1042

1043 1044 1045 1046
    // === Create matrices for FETI-DP preconditioner. ===

    if (fetiPreconditioner != FETI_NONE)
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1047
		      nRowsDual, nRowsDual, 30, PETSC_NULL,
1048 1049 1050 1051
		      &mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatCreateSeqAIJ(PETSC_COMM_SELF,
1052
		      nRowsInterior, nRowsInterior, 30, PETSC_NULL,
1053