PetscSolverGlobalMatrix.cc 30 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
Thomas Witkowski's avatar
Thomas Witkowski committed
20

21
22
#include <mpi.h>
#include "DirichletBC.h"
23
#include "DOFVector.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
24
25
26
#include "parallel/PetscSolverGlobalMatrix.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"
27
#include "solver/PetscTypes.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
28

29
namespace AMDiS { namespace Parallel {
30
31
32
33
34
35

  PetscSolverGlobalMatrix::PetscSolverGlobalMatrix(string name)
    : PetscSolver(name),
      zeroStartVector(false),
      printMatInfo(false)
  {
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
    PetscParameters params;
    bool matSolverPackage = false;
    
    // set the solver
    std::string solverName = "petsc";
    Parameters::get(name, solverName);
    if (solverName == "petsc") 
      Parameters::get(name + "->ksp_type", solverName);
    
    std::string kspSolver = params.solverMap[solverName];
    
    if (params.matSolverPackage[kspSolver]) {
      // direct solvers
      PetscOptionsInsertString(("-" + kspPrefix + "ksp_type preonly").c_str());
      PetscOptionsInsertString(("-" + kspPrefix + "pc_type lu").c_str());
      PetscOptionsInsertString(("-" + kspPrefix + "pc_factor_mat_solver_package " + (kspSolver != "direct" ? kspSolver : "mumps")).c_str());
      setMaxIterations(1);
      zeroStartVector = true;
      matSolverPackage = true;
    } else if (!params.emptyParam[kspSolver]) {    
      // other solvers
      PetscOptionsInsertString(("-" + kspPrefix + "ksp_type " + kspSolver).c_str());
    }
    
    
    // set the preconditioner
    string precon = "";
    Parameters::get(name + "->pc_type", precon);
    if (!precon.size())
      Parameters::get(name + "->left precon", precon);
    if (!matSolverPackage && !params.emptyParam[precon]) {
      precon = (params.preconMap.find(precon) != params.preconMap.end() ? params.preconMap[precon] : precon);
      PetscOptionsInsertString(("-" + kspPrefix + "pc_type " + precon).c_str());
    }
    
    PetscOptionsInsertString(("-" + kspPrefix + "ksp_max_it " + boost::lexical_cast<std::string>(getMaxIterations())).c_str());
    PetscOptionsInsertString(("-" + kspPrefix + "ksp_rtol " + boost::lexical_cast<std::string>(getRelative())).c_str());    
    PetscOptionsInsertString(("-" + kspPrefix + "ksp_atol " + boost::lexical_cast<std::string>(getTolerance())).c_str());   
    
    if (getInfo() >= 20)
      PetscOptionsInsertString(("-" + kspPrefix + "ksp_monitor_true_residual").c_str());
    else if (getInfo() >= 10)
      PetscOptionsInsertString(("-" + kspPrefix + "ksp_monitor").c_str());
    
    if (!matSolverPackage) {
      Parameters::get("parallel->use zero start vector", zeroStartVector);
    }
83
84
85
86
    Parameters::get("parallel->print matrix info", printMatInfo);
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
87
  void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
88
89
90
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
91
92
93
94
95
96
    TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT_DBG(interiorMap)("No parallel mapping object defined!\n");
    TEST_EXIT_DBG(seqMat)("No DOF matrix defined!\n");
    
    double wtime = MPI::Wtime();

97
    createMatVec(*seqMat);
98

Thomas Witkowski's avatar
Thomas Witkowski committed
99
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
100
      fillPetscMatrixWithCoarseSpace(seqMat);
101
102
      return;
    }
103
    
104
105
    // === Create PETSc vector (solution and a temporary vector). ===

Thomas Witkowski's avatar
Thomas Witkowski committed
106
107
108
109
110
111
#if (DEBUG != 0)
    MSG("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif

    // === Transfer values from DOF matrices to the PETSc matrix. === 

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
112
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
113
114
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
115
116
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
117
118
119
120
121

#if (DEBUG != 0)
    MSG("Fill petsc matrix 2 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif

122
    matAssembly();
123

124
125
    if (printMatInfo) {
      MatInfo matInfo;
126
      MatGetInfo(getMatInterior(), MAT_GLOBAL_SUM, &matInfo);
127
128
129
130
131
132
133
      MSG("Matrix info:\n");
      MSG("  memory usage: %e MB\n", matInfo.memory / (1024.0 * 1024.0));
      MSG("  mallocs: %d\n", static_cast<int>(matInfo.mallocs));
      MSG("  nz allocated: %d\n", static_cast<int>(matInfo.nz_allocated));
      MSG("  nz used: %d\n", static_cast<int>(matInfo.nz_used));
      MSG("  nz unneeded: %d\n", static_cast<int>(matInfo.nz_unneeded));
    }
134
135


136
    // === Init PETSc solver and preconditioner objects. ===
137

138
139
    initSolver(kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
140
    initPreconditioner(pcInterior);
141

142

Thomas Witkowski's avatar
Thomas Witkowski committed
143
144
145
#if (DEBUG != 0)
    MSG("Fill petsc matrix 3 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif
146
147


148
149
150
151
152
153
154
155
156
157
158
    // === For debugging allow to write the matrix to a file. ===

    bool dbgWriteMatrix = false;
    Parameters::get("parallel->debug->write matrix", dbgWriteMatrix);
    if (dbgWriteMatrix) {
      PetscViewer matView;
      PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mpi.mat",
			    FILE_MODE_WRITE, &matView);
      MatView(getMatInterior(), matView);
      PetscViewerDestroy(&matView);
    }
159
160
161
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
162
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *seqMat)
163
164
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
165
166

    TEST_EXIT_DBG(interiorMap)("Should not happen!\n");
167
168
    TEST_EXIT_DBG(coarseSpaceMap.size() == seqMat->getSize())
      ("Wrong sizes %d %d\n", coarseSpaceMap.size(), seqMat->getSize());
169

170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
    // === Prepare traverse of sequentially created matrices. ===

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

186
187
188
    bool localMatrix = 
      (meshDistributor->getMeshLevelData().getMpiComm(meshLevel) == MPI::COMM_SELF);

189
190
191
    // === Traverse all sequentially created matrices and add the values to ===
    // === the global PETSc matrices.                                       ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
192
    int nComponents = seqMat->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
193
194
195
196
197
    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
	DOFMatrix* dofMat = (*seqMat)[rowComponent][colComponent];

	if (!dofMat)
198
199
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
200
201
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
202

Thomas Witkowski's avatar
Thomas Witkowski committed
203
204
	std::set<DegreeOfFreedom> &dirichletRows = dofMat->getDirichletRows();

Thomas Witkowski's avatar
Thomas Witkowski committed
205
206
	traits::col<Matrix>::type col(dofMat->getBaseMatrix());
	traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
207
208
	
	// Traverse all rows.
Thomas Witkowski's avatar
Thomas Witkowski committed
209
210
	for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()), 
	       cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
211

212
	  bool isRowCoarse = isCoarseSpace(rowComponent, cursor.value());
Thomas Witkowski's avatar
Thomas Witkowski committed
213
214
215

	  // For the case, this is a dirichlet row we have to check whether the 
	  // rank is also owner of this row DOF.
216
217
218
	  if (dirichletRows.count(cursor.value())) {
	    if ((!isRowCoarse && !(*interiorMap)[rowComponent].isRankDof(cursor.value())) ||
		(isRowCoarse && !(*rowCoarseSpace)[rowComponent].isRankDof(cursor.value())))
219
	      continue;    
Thomas Witkowski's avatar
Thomas Witkowski committed
220
	  }
221
222
223
224
225
226
227
228
229
230
  
	  cols.clear();
	  colsOther.clear();
	  values.clear();	  
	  valuesOther.clear();

	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

Thomas Witkowski's avatar
Thomas Witkowski committed
231
	    bool isColCoarse = isCoarseSpace(colComponent, col(*icursor));
232

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
233
	    if (isColCoarse == false)
234
	      if ((*interiorMap)[colComponent].isSet(col(*icursor)) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
235
236
		continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
237
238
239
	    if (isColCoarse == isRowCoarse) {
	      cols.push_back(col(*icursor));
	      values.push_back(value(*icursor));
240
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
241
242
	      colsOther.push_back(col(*icursor));
	      valuesOther.push_back(value(*icursor));
243
244
245
246
247
248
	    }
	  }  // for each nnz in row


	  // === Set matrix values. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
249
	  if (isRowCoarse) {
Thomas Witkowski's avatar
Thomas Witkowski committed
250
251
	    for (unsigned int i = 0; i < cols.size(); i++)
	      cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]);
252

253
	    int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, cursor.value());
Thomas Witkowski's avatar
Thomas Witkowski committed
254
	    MatSetValues(getMatCoarseByComponent(rowComponent, colComponent),
255
256
	    		 1, &rowIndex, cols.size(),
	    		 &(cols[0]), &(values[0]), ADD_VALUES);
257
258

	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
259
260
261
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  interiorMap->getMatIndex(colComponent, colsOther[i]) + 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
262
		  rStartInterior;	      
263

Thomas Witkowski's avatar
Thomas Witkowski committed
264
	      MatSetValues(getMatCoarseInteriorByComponent(rowComponent), 
265
266
	       		   1, &rowIndex, colsOther.size(),
 	       		   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
267
268
	    }
	  } else {
269
	    if ((*interiorMap)[rowComponent].isSet(cursor.value()) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
270
271
	      continue;

272
	    int localRowIndex = 
273
	      (localMatrix ? 
274
275
	       interiorMap->getLocalMatIndex(rowComponent, cursor.value()) :
	       interiorMap->getMatIndex(rowComponent, cursor.value()));
276

Thomas Witkowski's avatar
Thomas Witkowski committed
277
	    for (unsigned int i = 0; i < cols.size(); i++) {
278
	      if (localMatrix)
Thomas Witkowski's avatar
Thomas Witkowski committed
279
		cols[i] = interiorMap->getLocalMatIndex(colComponent, cols[i]);
280
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
281
		cols[i] = interiorMap->getMatIndex(colComponent, cols[i]);
282
283
	    }
	    
284
285
	    MatSetValues(getMatInterior(), 1, &localRowIndex, cols.size(),
  	     		 &(cols[0]), &(values[0]), ADD_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
286
	    
287
	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
288
	      int globalRowIndex = 
289
		interiorMap->getMatIndex(rowComponent, cursor.value()) + rStartInterior;
Thomas Witkowski's avatar
Thomas Witkowski committed
290
      
Thomas Witkowski's avatar
Thomas Witkowski committed
291
292
293
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  colCoarseSpace->getMatIndex(colComponent, colsOther[i]);
294

295
296
297
	      MatSetValues(getMatInteriorCoarseByComponent(colComponent), 
			   1, &globalRowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
298
299
300
301
302
303
	    }
	  }
	} 
      }
    }

304
    matAssembly();
305

306
307
    // === Create solver for the non primal (thus local) variables. ===

308
    KSPCreate(domainComm, &kspInterior);
309
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(),
310
		    SAME_NONZERO_PATTERN);
311
312
313
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
314
315
    if (isSymmetric) {
      PCSetType(pcInterior, PCCHOLESKY);
316
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
317
318
    } else {
      PCSetType(pcInterior, PCLU);
319
      if (localMatrix)
320
321
322
323
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
      else
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    }
324
    KSPSetFromOptions(kspInterior);
325
326
327
  }


328
329
330
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
331

332
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
333
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
334
    
Thomas Witkowski's avatar
Thomas Witkowski committed
335
    // === Transfer values from DOF vector to the PETSc vector. === 
Thomas Witkowski's avatar
Thomas Witkowski committed
336
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
337
      for (int i = 0; i < vec->getSize(); i++)
338
	setDofVector(getVecRhsInterior(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
339
		     getVecRhsCoarseByComponent(i), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
340
341
    } else {
      for (int i = 0; i < vec->getSize(); i++)
342
	setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
343
344
    }

345
    vecRhsAssembly();
346
347
348
349
350
351
352
353
354
355
356
357

    // === For debugging allow to write the rhs vector to a file. ===

    bool dbgWriteRhs = false;
    Parameters::get("parallel->debug->write rhs", dbgWriteRhs);
    if (dbgWriteRhs) {
      PetscViewer vecView;
      PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mpi.vec",
			    FILE_MODE_WRITE, &vecView);
      VecView(getVecRhsInterior(), vecView);
      PetscViewerDestroy(&vecView);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
358
359
360
  }


361
362
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
363
364
365
366
367
368
369
  {
    FUNCNAME("PetscSolverGlobalMatrix::solvePetscMatrix()");

    int nComponents = vec.getSize();

    // === Set old solution to be initiual guess for PETSc solver. ===
    if (!zeroStartVector) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
370
371
      TEST_EXIT(coarseSpaceMap.size() == 0)("Not yet supported!\n");

372
      VecSet(getVecSolInterior(), 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
373
374
      
      for (int i = 0; i < nComponents; i++)
375
	setDofVector(getVecSolInterior(), vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
376

377
      vecSolAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
378
379
    }

380

Thomas Witkowski's avatar
Thomas Witkowski committed
381
382
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
383
384
385
    if (nullspace.size() > 0 || 
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
386
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
387

388
389
390
391
392
393
394
395
396
397
      if (constNullspaceComponent.size() > 0) {
	nullspace.clear();
	SystemVector *basisVec = new SystemVector(vec);
	basisVec->set(0.0);
	for (unsigned int i = 0; i < constNullspaceComponent.size(); i++)
	  basisVec->getDOFVector(constNullspaceComponent[i])->set(1.0);

	nullspace.push_back(basisVec);
      } 

Thomas Witkowski's avatar
Thomas Witkowski committed
398
      if (nullspace.size() > 0) {
399
	VecDuplicate(getVecSolInterior(), &nullspaceBasis);
400
	setDofVector(nullspaceBasis, *(nullspace[0]), true);
Thomas Witkowski's avatar
Thomas Witkowski committed
401
402
403
	
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
404
405

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
406
	
407
	MatNullSpaceCreate(domainComm, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
Thomas Witkowski's avatar
Thomas Witkowski committed
408
409
			   1, &nullspaceBasis, &matNullspace);

410
	MatMult(getMatInterior(), nullspaceBasis, getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
411
	PetscReal n;
412
	VecNorm(getVecSolInterior(), NORM_2, &n);
Thomas Witkowski's avatar
Thomas Witkowski committed
413
414
	MSG("NORM IS: %e\n", n);
      } else {
415
	MatNullSpaceCreate(domainComm, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
416
      }
417

Thomas Witkowski's avatar
Thomas Witkowski committed
418
419
420
      MSG("NULLSPACE IS NOT REMOVED!\n");
      // MatSetNullSpace(getMatInterior(), matNullspace);
      // KSPSetNullSpace(kspInterior, matNullspace);
421

Thomas Witkowski's avatar
Thomas Witkowski committed
422
      // === Remove null space, if requested. ===
423

Thomas Witkowski's avatar
Thomas Witkowski committed
424
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
425
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
426
427

	MSG("Remove nullspace from rhs vector.\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
428
	
429
	MatNullSpaceRemove(matNullspace, getVecRhsInterior(), PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
430
431
432
433
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
434
435
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
436
    // PETSc.
437
    solve(getVecRhsInterior(), getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
438

439
440

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
441
      MatNullSpaceDestroy(&matNullspace);
442
443
444
445
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
446
447
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
448
    VecGetArray(getVecSolInterior(), &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
449

450
    int c = 0;
451
452
    for (int component = 0; component < nComponents; component++) {
      DOFVector<double> &dv = *(vec.getDOFVector(component));
453

454
      DofMap& d = (*interiorMap)[component].getMap();
455
456
457
      for (DofMap::iterator it = d.begin(); it != d.end(); ++it)
	if (it->second.local != -1)
	  dv[it->first] = vecPointer[c++];
Thomas Witkowski's avatar
Thomas Witkowski committed
458
459
    }

460
    VecRestoreArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
461
462
463
464

    // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to ===  
    // === more than one partition.                                       ===  
    meshDistributor->synchVector(vec);
465
466
467
  }


468
469
470
471
472
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()");

    Vec tmp;
473
    if (domainComm.Get_size() == 1)
474
      createLocalVec(*interiorMap, tmp);
475
    else
476
      createVec(*interiorMap, tmp);
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502

    PetscScalar *tmpValues, *rhsValues;
    VecGetArray(tmp, &tmpValues);
    VecGetArray(rhs, &rhsValues);

    for (int i = 0; i < interiorMap->getRankDofs(); i++)
      tmpValues[i] = rhsValues[i];

    VecRestoreArray(rhs, &rhsValues);
    VecRestoreArray(tmp, &tmpValues);

    KSPSolve(kspInterior, tmp, tmp);

    VecGetArray(tmp, &tmpValues);
    VecGetArray(sol, &rhsValues);

    for (int i = 0; i < interiorMap->getRankDofs(); i++) 
      rhsValues[i] = tmpValues[i];

    VecRestoreArray(sol, &rhsValues);
    VecRestoreArray(tmp, &tmpValues);

    VecDestroy(&tmp);
  }


503
504
505
506
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");

507
    matDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
508

509
510
511
    exitPreconditioner(pcInterior);

    exitSolver(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
512
513
514
  }


515
516
517
518
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()");

519
    vecDestroy();
520
521
522
  }


523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createFieldSplit()");

    vector<string> isNames;
    Parameters::get("parallel->solver->is blocks", isNames);

    int nBlocks = isNames.size();
    if (nBlocks == 0)
      return;

    for (int i = 0; i < nBlocks; i++) {
      MSG("Create for block %s\n", isNames[i].c_str());

      vector<int> blockComponents;
      Parameters::get("parallel->solver->is block " + lexical_cast<string>(i),
		      blockComponents);
      int nComponents = static_cast<int>(blockComponents.size());

542
      TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
543
544
545
546
547
548
549
550

      // Check if blocks are continous
      for (int j = 0; j < nComponents; j++) {
	TEST_EXIT(blockComponents[j] == blockComponents[0] + j)
	  ("Does not yet support not continous IS blocks! Block %s\n", 
	   isNames[i].c_str());
      }

551
      createFieldSplit(pc, isNames[i].c_str(), blockComponents);
552
553
554
555
    }
  }


556
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc, 
557
						 const char* splitName, 
558
559
560
561
562
563
						 vector<int> &components)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createFieldSplit()");

    IS is;
    interiorMap->createIndexSet(is, components[0], components.size());
564
    PCFieldSplitSetIS(pc, splitName, is);
565
566
    ISDestroy(&is);
  }
567
568
569
570
571
572
573
574
575
576
  
  
  void PetscSolverGlobalMatrix::extractVectorComponent(Vec input, int i, Vec *output, int numberOfComponents)
  {
    FUNCNAME("PetscSolverGlobalMatrix::extractVectorComponent()");
    IS is;
    interiorMap->createIndexSet(is, i, numberOfComponents);
    VecGetSubVector(input, is, output);
    ISDestroy(&is);
  }
577

578
579
580
581
582
583
584
585
586
587
588
  void PetscSolverGlobalMatrix::extractMatrixComponent(Mat input, int startRow, int numberOfRows, int startCol, int numberOfCols, Mat *output)
  {
    FUNCNAME("PetscSolverGlobalMatrix::extractMatrixComponent()");
    IS isrow, iscol;
    interiorMap->createIndexSet(isrow, startRow, numberOfRows);
    interiorMap->createIndexSet(iscol, startCol, numberOfCols);
    MatGetSubMatrix(input, isrow, iscol, MAT_INITIAL_MATRIX, output);
    ISDestroy(&iscol);
    ISDestroy(&isrow);
  }
  
589

590
  void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
591
592
593
  {
    FUNCNAME("PetscSolverGlobalMatrix::initSolver()");

594
    KSPCreate(domainComm, &ksp);
595
596
597
598
599
600
601
602
603
604
605
606
607
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), 
		    SAME_NONZERO_PATTERN); 
    KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(ksp, KSPBCGS);
    KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
    KSPSetFromOptions(ksp);

    // Do not delete the solution vector, use it for the initial guess.
    if (!zeroStartVector)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
  }  


608
  void PetscSolverGlobalMatrix::exitSolver(KSP &ksp)
609
610
611
612
613
614
615
  {
    FUNCNAME("PetscSolverGlobalMatrix::exitSolver()");

    KSPDestroy(&ksp);
  }  


616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
  void PetscSolverGlobalMatrix::initPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::initPreconditioner()");

    PCSetFromOptions(pc);
    createFieldSplit(pc);
  }

 
  void PetscSolverGlobalMatrix::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::exitPreconditioner()");
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
631
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
Thomas Witkowski's avatar
Thomas Witkowski committed
632
					     int rowComp, int colComp)
Thomas Witkowski's avatar
Thomas Witkowski committed
633
634
635
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
636
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
637
638

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
639
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
640
641
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
642
643
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
644
645
646
647
648
649
650
651
652
653
654

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

    vector<int> cols;
    vector<double> values;
    cols.reserve(300);
    values.reserve(300);
    
    vector<int> globalCols;

655
656
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
657
658
    std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows();
      
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
659
660
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
661
662

    // === Traverse all rows of the DOF matrix and insert row wise the values ===
Thomas Witkowski's avatar
Thomas Witkowski committed
663
664
    // === to the PETSc matrix.                                               ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
665
666
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
667
      // Global index of the current row DOF.
668
      MultiIndex rowMultiIndex;
669
      if ((*interiorMap)[rowComp].find(cursor.value(), rowMultiIndex) == false)
670
671
672
	continue;

      int globalRowDof = rowMultiIndex.global;
673

Thomas Witkowski's avatar
Thomas Witkowski committed
674
      // Test if the current row DOF is a periodic DOF.
675
      bool periodicRow = perMap.isPeriodic(rowFe, globalRowDof);
676

Thomas Witkowski's avatar
Thomas Witkowski committed
677
      // Dirichlet rows can be set only be the owner ranks.
678
      if (dirichletRows.count(cursor.value()) && !((*interiorMap)[rowComp].isRankDof(cursor.value())))
Thomas Witkowski's avatar
Thomas Witkowski committed
679
680
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
681
682
683
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

684
	// Get PETSc's mat row index.
Thomas Witkowski's avatar
Thomas Witkowski committed
685
	int rowIndex = interiorMap->getMatIndex(rowComp, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
686
687
688
689
690
691
692
693

	cols.clear();
	values.clear();

	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	     icursor != icend; ++icursor) {

	  // Global index of the current column index.
694
	  MultiIndex colMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
695
	  if ((*interiorMap)[colComp].find(col(*icursor), colMultiIndex) == false)
696
697
698
	    continue;

	  int globalColDof = colMultiIndex.global;
Thomas Witkowski's avatar
Thomas Witkowski committed
699
	  // Test if the current col dof is a periodic dof.
700
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
701
	  // Get PETSc's mat col index.
Thomas Witkowski's avatar
Thomas Witkowski committed
702
	  int colIndex = interiorMap->getMatIndex(colComp, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
703
704
705
706
707
708
709

	  // Ignore all zero entries, expect it is a diagonal entry.
 	  if (value(*icursor) == 0.0 && rowIndex != colIndex)
 	    continue;

	  if (!periodicCol) {
	    // Calculate the exact position of the column index in the PETSc matrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
710
711
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
712
713
714
715
716
	  } else {
	    // === Row index is not periodic, but column index is. ===

	    // Create set of all periodic associations of the column index.
	    std::set<int> perAsc;
717
718
719
	    perMap.fillAssociations(colFe, globalColDof, 
				    meshDistributor->getElementObjectDb(), perAsc);

Thomas Witkowski's avatar
Thomas Witkowski committed
720
721
722
723
724
725
726
727
728
	    // Scale value to the number of periodic associations of the column index.
	    double scaledValue = 
	      value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));

	    
	    // === Create set of all matrix column indices due to the periodic ===
	    // === associations of the column DOF index.                       ===

	    vector<int> newCols;
729
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
730
	    for (unsigned int i = 0; i < newCols.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
731
	      cols.push_back(interiorMap->getMatIndex(colComp, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
732
733
734
735
736
	      values.push_back(scaledValue);	      
	    }
	  }
	}

737
  	MatSetValues(getMatInterior(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
738
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
739
740
741
742
743
744
745
746
747
748
749
750
751
752
      } else {
	// === Row DOF index is periodic. ===

	// Because this row is periodic, we will have to add the entries of this 
	// matrix row to multiple rows. The following maps store to each row an
	// array of column indices and values of the entries that must be added to
	// the PETSc matrix.
	map<int, vector<int> > colsMap;
	map<int, vector<double> > valsMap;

	// Traverse all column entries.
	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	     icursor != icend; ++icursor) {
	  // Global index of the current column index.
Thomas Witkowski's avatar
Thomas Witkowski committed
753
	  int globalColDof = (*interiorMap)[colComp][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
754
755
756
757
758
759
760
761
762

	  // Ignore all zero entries, expect it is a diagonal entry.
 	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
 	    continue;

	  // === Add all periodic associations of both, the row and the column ===
	  // === indices to the set perAsc.                                    ===

	  std::set<int> perAsc;
763
764
765
766
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
767
768
769
770
771
772
773
774
775
776

	  // Scale the value with respect to the number of periodic associations.
	  double scaledValue = 
	    value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));


	  // === Create all matrix entries with respect to the periodic  ===
	  // === associations of the row and column indices.             ===

	  vector<pair<int, int> > entry;
Thomas Witkowski's avatar
Blbu  
Thomas Witkowski committed
777
778
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
779
780
781

	  // === Translate the matrix entries to PETSc's matrix.

782
	  for (unsigned int i = 0; i < entry.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
783
784
	    int rowIdx = interiorMap->getMatIndex(rowComp, entry[i].first);
	    int colIdx = interiorMap->getMatIndex(colComp, entry[i].second);
Thomas Witkowski's avatar
Thomas Witkowski committed
785

786
787
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
788
789
790
791
792
793
794
795
796
797
798
799
	  }
	}


	// === Finally, add all periodic rows to the PETSc matrix. ===

	for (map<int, vector<int> >::iterator rowIt = colsMap.begin();
	     rowIt != colsMap.end(); ++rowIt) {
	  TEST_EXIT_DBG(rowIt->second.size() == valsMap[rowIt->first].size())
	    ("Should not happen!\n");

	  int rowIndex = rowIt->first;
800

801
	  MatSetValues(getMatInterior(), 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
802
803
804
805
806
807
808
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
809
810
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
811
					     DOFVector<double>* vec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
812
					     int rowComp, 
813
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
814
815
816
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

817
    const FiniteElemSpace *feSpace = vec->getFeSpace();
818
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
819

820
    ParallelDofMapping *rowCoarseSpace = 
Thomas Witkowski's avatar
Thomas Witkowski committed
821
822
823
      (coarseSpaceMap.size() ? coarseSpaceMap[rowComp] : NULL);

    map<DegreeOfFreedom, double> &dirichletValues = vec->getDirichletValues();
824

Thomas Witkowski's avatar
Thomas Witkowski committed
825
826
827
    // Traverse all used DOFs in the dof vector.
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
828
829
      
      DegreeOfFreedom dof = dofIt.getDOFIndex();
830

Thomas Witkowski's avatar
Thomas Witkowski committed
831
832
      if (rankOnly && !(*interiorMap)[rowComp].isRankDof(dof))
	continue;
833

834
835
      bool isCoarseDof = isCoarseSpace(rowComp, dof);

Thomas Witkowski's avatar
Thomas Witkowski committed
836
      // Dirichlet rows can be set only be the owner ranks.
837
838
839
840
841
      if (dirichletValues.count(dof)) {
	if ((!isCoarseDof && !((*interiorMap)[rowComp].isRankDof(dof))) ||
	    (isCoarseDof && !((*rowCoarseSpace)[rowComp].isRankDof(dof))))
	  continue;
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
842

843
      if (isCoarseDof) {
844
845
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
846
	int index = rowCoarseSpace->getMatIndex(rowComp, dof);
847
848
	VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES);
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
849
	if ((*interiorMap)[rowComp].isSet(dof) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
850
851
	  continue;

852
	// Calculate global row index of the DOF.
Thomas Witkowski's avatar
Thomas Witkowski committed
853
	DegreeOfFreedom globalRowDof = (*interiorMap)[rowComp][dof].global;
854
855
856
857
858
	
	// Get PETSc's mat index of the row DOF.
	int index = 0;
	if (interiorMap->isMatIndexFromGlobal())
	  index = 
Thomas Witkowski's avatar
Thomas Witkowski committed
859
	    interiorMap->getMatIndex(rowComp, globalRowDof) + rStartInterior;
860
861
	else
	  index =
Thomas Witkowski's avatar
Thomas Witkowski committed
862
	    interiorMap->getMatIndex(rowComp, dof) + rStartInterior;
863

864
	if (perMap.isPeriodic(feSpace, globalRowDof)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
865
866
867
	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
	  double value = *dofIt / (perAsc.size() + 1.0);
	  VecSetValue(vecInterior, index, value, ADD_VALUES);
868

Thomas Witkowski's avatar
Thomas Witkowski committed
869
870
871
	  for (std::set<int>::iterator perIt = perAsc.begin(); 
	       perIt != perAsc.end(); ++perIt) {
	    int mappedDof = perMap.map(feSpace, *perIt, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
872
	    int mappedIndex = interiorMap->getMatIndex(rowComp, mappedDof);
873

Thomas Witkowski's avatar
Thomas Witkowski committed
874
	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
875
876
877
	  }	  
	} else {	  
	  // The DOF index is not periodic.
Thomas Witkowski's avatar
Thomas Witkowski committed
878
879
	  VecSetValue(vecInterior, index, *dofIt, ADD_VALUES);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
880
881
882
883
      }
    }
  }

884
885
886
887
888
889
890
891
892
893

  PetscSolver* PetscSolverGlobalMatrix::createSubSolver(int component,
							string kspPrefix)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createSubSolver()");

    vector<const FiniteElemSpace*> fe;
    fe.push_back(componentSpaces[component]);

    PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
894
    subSolver->setKspPrefix(kspPrefix);
895
    subSolver->setMeshDistributor(meshDistributor, 0);
896
897
898
899
900
901
902
903
904
    subSolver->init(fe, fe);

    ParallelDofMapping &subDofMap = subSolver->getDofMap();
    subDofMap[0] = dofMap[component];
    subDofMap.update();

    return subSolver;
  }

905
906
907
908

  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp,
						     int constFeSpace,
						     bool test)
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
  {
    FUNCNAME("PetscSolverGlobalMatrix::setConstantNullSpace()");

    Vec nullSpaceBasis;
    VecDuplicate(getVecSolInterior(), &nullSpaceBasis);
    
    SystemVector basisVec("tmp", componentSpaces, componentSpaces.size(), true);
    basisVec.set(0.0);
    basisVec.getDOFVector(constFeSpace)->set(1.0);
    
    setDofVector(nullSpaceBasis, basisVec, true);
    VecAssemblyBegin(nullSpaceBasis);
    VecAssemblyEnd(nullSpaceBasis);
    VecNormalize(nullSpaceBasis, PETSC_NULL);
    
    if (test) {
      Vec tmp;
      MatGetVecs(getMatInterior(), &tmp, PETSC_NULL);
      MatMult(getMatInterior(), nullSpaceBasis, tmp);
      PetscReal n;
      VecNorm(tmp, NORM_2, &n);
      MSG("NORM IS: %e\n", n);
      VecDestroy(&tmp);
932
    }
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
933
934
935
    
    
    MatNullSpace matNullSpace;
936
    MatNullSpaceCreate(domainComm, PETSC_FALSE, 1, &nullSpaceBasis, &matNullSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
937
938
    KSPSetNullSpace(ksp, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
939
940

    VecDestroy(&nullSpaceBasis);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
941
942
  }

943

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
944
945
946
947
948
  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp)
  {
    FUNCNAME("PetscSolverGlobalMatrix::setConstantNullSpace()");

    MatNullSpace matNullSpace;
949
    MatNullSpaceCreate(domainComm, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
950
951
952
    KSPSetNullSpace(ksp, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
  }
953

954
} } // end namespace Parallel, AMDiS