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 30
using namespace std;

31
namespace AMDiS { namespace Parallel {
32

33
  PetscSolverGlobalMatrix::PetscSolverGlobalMatrix(string name, bool setOptions)
34 35 36
    : PetscSolver(name),
      zeroStartVector(false),
      printMatInfo(false)
37 38
  { FUNCNAME_DBG("PetscSolverGlobalMatrix()");
  
39
    bool matSolverPackage = false;
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
    if (setOptions) {
      PetscParameters params;
      
      // 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.find(kspSolver) != params.matSolverPackage.end()) {
	// direct solvers
	PetscOptionsInsertString(("-" + kspPrefix + "ksp_type preonly").c_str());
	PetscOptionsInsertString(("-" + kspPrefix + "pc_type lu").c_str());
	PetscOptionsInsertString(("-" + kspPrefix + "pc_factor_mat_solver_package " + kspSolver).c_str());
	setMaxIterations(1);
	zeroStartVector = true;
	matSolverPackage = true;
      } else if (params.emptyParam.find(kspSolver) == params.emptyParam.end() && solverName != "petsc") {    
	// 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.find(precon) == params.emptyParam.end()) {
	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());
82 83
    }
    if (!matSolverPackage) {
84
      Parameters::get(name + "->use zero start vector", zeroStartVector);
85
    }
86
    Parameters::get("parallel->print matrix info", printMatInfo);
87 88 89 90 91 92 93 94 95 96 97 98 99
    
#if DEBUG != 0
    bool printOptionsInfo = false;
    Parameters::get("parallel->debug->print options info", printOptionsInfo);
    if (printOptionsInfo) {
    MSG("PetscOptionsView:\n");
    PetscViewer viewer;
    PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
    PetscViewerSetType(viewer, PETSCVIEWERASCII);
    PetscOptionsView(viewer);
    PetscViewerDestroy(&viewer);
    }
#endif
100 101 102
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
103
  void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
104 105 106
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
107 108 109 110
    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");
    
111
#if (DEBUG != 0)
112
    Timer t;
113
#endif
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
114

115
    createMatVec(*seqMat);
116

Thomas Witkowski's avatar
Thomas Witkowski committed
117
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
118
      fillPetscMatrixWithCoarseSpace(seqMat);
119 120
      return;
    }
121
    
122 123
    // === Create PETSc vector (solution and a temporary vector). ===

Thomas Witkowski's avatar
Thomas Witkowski committed
124
#if (DEBUG != 0)
125 126
    MSG("Fill petsc matrix 1 needed %.5f seconds\n", t.elapsed());
    t.reset();
Thomas Witkowski's avatar
Thomas Witkowski committed
127 128 129 130
#endif

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
131
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
132 133
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
134 135
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
136 137

#if (DEBUG != 0)
138 139
    MSG("Fill petsc matrix 2 needed %.5f seconds\n", t.elapsed());
    t.reset();
Thomas Witkowski's avatar
Thomas Witkowski committed
140 141
#endif

142
    matAssembly();
143

144 145
    if (printMatInfo) {
      MatInfo matInfo;
146
      MatGetInfo(getMatInterior(), MAT_GLOBAL_SUM, &matInfo);
147 148 149 150 151 152 153
      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));
    }
154 155


156
    // === Init PETSc solver and preconditioner objects. ===
157

158 159
    initSolver(kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
160
    initPreconditioner(pcInterior);
161

162

Thomas Witkowski's avatar
Thomas Witkowski committed
163
#if (DEBUG != 0)
164
    MSG("Fill petsc matrix 3 needed %.5f seconds\n", t.elapsed());
Thomas Witkowski's avatar
Thomas Witkowski committed
165
#endif
166 167


168 169 170 171 172 173 174 175 176 177 178
    // === 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);
    }
179 180 181
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
182
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *seqMat)
183
  {
184
    FUNCNAME_DBG("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
185

186
    TEST_EXIT_DBG(interiorMap)("No interiorMap! Should not happen!\n");
187 188
    TEST_EXIT_DBG(coarseSpaceMap.size() == seqMat->getSize())
      ("Wrong sizes %d %d\n", coarseSpaceMap.size(), seqMat->getSize());
189

190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
    // === 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);

206 207 208
    bool localMatrix = 
      (meshDistributor->getMeshLevelData().getMpiComm(meshLevel) == MPI::COMM_SELF);

209 210 211
    // === Traverse all sequentially created matrices and add the values to ===
    // === the global PETSc matrices.                                       ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
212
    int nComponents = seqMat->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
213 214 215 216 217
    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
	DOFMatrix* dofMat = (*seqMat)[rowComponent][colComponent];

	if (!dofMat)
218 219
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
220 221
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
222

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

Thomas Witkowski's avatar
Thomas Witkowski committed
225 226
	traits::col<Matrix>::type col(dofMat->getBaseMatrix());
	traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
227 228
	
	// Traverse all rows.
Thomas Witkowski's avatar
Thomas Witkowski committed
229 230
	for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()), 
	       cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
231

232
	  bool isRowCoarse = isCoarseSpace(rowComponent, cursor.value());
Thomas Witkowski's avatar
Thomas Witkowski committed
233 234 235

	  // For the case, this is a dirichlet row we have to check whether the 
	  // rank is also owner of this row DOF.
236 237 238
	  if (dirichletRows.count(cursor.value())) {
	    if ((!isRowCoarse && !(*interiorMap)[rowComponent].isRankDof(cursor.value())) ||
		(isRowCoarse && !(*rowCoarseSpace)[rowComponent].isRankDof(cursor.value())))
239
	      continue;    
Thomas Witkowski's avatar
Thomas Witkowski committed
240
	  }
241 242 243 244 245 246 247 248 249 250
  
	  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
251
	    bool isColCoarse = isCoarseSpace(colComponent, col(*icursor));
252

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
253
	    if (isColCoarse == false)
254
	      if ((*interiorMap)[colComponent].isSet(col(*icursor)) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
255 256
		continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
257 258 259
	    if (isColCoarse == isRowCoarse) {
	      cols.push_back(col(*icursor));
	      values.push_back(value(*icursor));
260
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
261 262
	      colsOther.push_back(col(*icursor));
	      valuesOther.push_back(value(*icursor));
263 264 265 266 267 268
	    }
	  }  // for each nnz in row


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

Thomas Witkowski's avatar
Thomas Witkowski committed
269
	  if (isRowCoarse) {
Thomas Witkowski's avatar
Thomas Witkowski committed
270 271
	    for (unsigned int i = 0; i < cols.size(); i++)
	      cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]);
272

273
	    int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, cursor.value());
Thomas Witkowski's avatar
Thomas Witkowski committed
274
	    MatSetValues(getMatCoarseByComponent(rowComponent, colComponent),
275 276
	    		 1, &rowIndex, cols.size(),
	    		 &(cols[0]), &(values[0]), ADD_VALUES);
277 278

	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
279 280 281
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  interiorMap->getMatIndex(colComponent, colsOther[i]) + 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
282
		  rStartInterior;	      
283

Thomas Witkowski's avatar
Thomas Witkowski committed
284
	      MatSetValues(getMatCoarseInteriorByComponent(rowComponent), 
285 286
	       		   1, &rowIndex, colsOther.size(),
 	       		   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
287 288
	    }
	  } else {
289
	    if ((*interiorMap)[rowComponent].isSet(cursor.value()) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
290 291
	      continue;

292
	    int localRowIndex = 
293
	      (localMatrix ? 
294 295
	       interiorMap->getLocalMatIndex(rowComponent, cursor.value()) :
	       interiorMap->getMatIndex(rowComponent, cursor.value()));
296

Thomas Witkowski's avatar
Thomas Witkowski committed
297
	    for (unsigned int i = 0; i < cols.size(); i++) {
298
	      if (localMatrix)
Thomas Witkowski's avatar
Thomas Witkowski committed
299
		cols[i] = interiorMap->getLocalMatIndex(colComponent, cols[i]);
300
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
301
		cols[i] = interiorMap->getMatIndex(colComponent, cols[i]);
302 303
	    }
	    
304 305
	    MatSetValues(getMatInterior(), 1, &localRowIndex, cols.size(),
  	     		 &(cols[0]), &(values[0]), ADD_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
306
	    
307
	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
308
	      int globalRowIndex = 
309
		interiorMap->getMatIndex(rowComponent, cursor.value()) + rStartInterior;
Thomas Witkowski's avatar
Thomas Witkowski committed
310
      
Thomas Witkowski's avatar
Thomas Witkowski committed
311 312 313
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  colCoarseSpace->getMatIndex(colComponent, colsOther[i]);
314

315 316 317
	      MatSetValues(getMatInteriorCoarseByComponent(colComponent), 
			   1, &globalRowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
318 319 320 321 322 323
	    }
	  }
	} 
      }
    }

324
    matAssembly();
325

326 327
    // === Create solver for the non primal (thus local) variables. ===

328
    KSPCreate(domainComm, &kspInterior);
329
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(),
330
		    SAME_NONZERO_PATTERN);
331 332 333
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
334 335
    if (isSymmetric) {
      PCSetType(pcInterior, PCCHOLESKY);
336
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
337 338
    } else {
      PCSetType(pcInterior, PCLU);
339
      if (localMatrix)
340 341 342 343
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
      else
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    }
344
    KSPSetFromOptions(kspInterior);
345 346 347
  }


348 349
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
350
    FUNCNAME_DBG("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
351

352
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
353
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
354
    
Thomas Witkowski's avatar
Thomas Witkowski committed
355
    // === Transfer values from DOF vector to the PETSc vector. === 
Thomas Witkowski's avatar
Thomas Witkowski committed
356
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
357
      for (int i = 0; i < vec->getSize(); i++)
358
	setDofVector(getVecRhsInterior(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
359
		     getVecRhsCoarseByComponent(i), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
360 361
    } else {
      for (int i = 0; i < vec->getSize(); i++)
362
	setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
363 364
    }

365
    vecRhsAssembly();
366 367 368 369 370 371 372 373 374 375 376 377

    // === 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
378 379 380
  }


381 382
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
383 384 385 386 387 388 389
  {
    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
390 391
      TEST_EXIT(coarseSpaceMap.size() == 0)("Not yet supported!\n");

392
      VecSet(getVecSolInterior(), 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
393 394
      
      for (int i = 0; i < nComponents; i++)
395
	setDofVector(getVecSolInterior(), vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
396

397
      vecSolAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
398 399
    }

400

Thomas Witkowski's avatar
Thomas Witkowski committed
401 402
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
403 404 405
    if (nullspace.size() > 0 || 
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
406
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
407

408 409 410 411 412 413 414 415 416 417
      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
418
      if (nullspace.size() > 0) {
419
	VecDuplicate(getVecSolInterior(), &nullspaceBasis);
420
	setDofVector(nullspaceBasis, *(nullspace[0]), true);
Thomas Witkowski's avatar
Thomas Witkowski committed
421 422 423
	
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
424 425

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
426
	
427
	MatNullSpaceCreate(domainComm, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
Thomas Witkowski's avatar
Thomas Witkowski committed
428 429
			   1, &nullspaceBasis, &matNullspace);

430
	MatMult(getMatInterior(), nullspaceBasis, getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
431
	PetscReal n;
432
	VecNorm(getVecSolInterior(), NORM_2, &n);
Thomas Witkowski's avatar
Thomas Witkowski committed
433 434
	MSG("NORM IS: %e\n", n);
      } else {
435
	MatNullSpaceCreate(domainComm, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
436
      }
437

Thomas Witkowski's avatar
Thomas Witkowski committed
438 439 440
      MSG("NULLSPACE IS NOT REMOVED!\n");
      // MatSetNullSpace(getMatInterior(), matNullspace);
      // KSPSetNullSpace(kspInterior, matNullspace);
441

Thomas Witkowski's avatar
Thomas Witkowski committed
442
      // === Remove null space, if requested. ===
443

Thomas Witkowski's avatar
Thomas Witkowski committed
444
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
445
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
446 447

	MSG("Remove nullspace from rhs vector.\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
448
	
449
	MatNullSpaceRemove(matNullspace, getVecRhsInterior(), PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
450 451 452 453
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
454 455
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
456
    // PETSc.
457
    solve(getVecRhsInterior(), getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
458

459 460

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
461
      MatNullSpaceDestroy(&matNullspace);
462 463 464 465
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
466 467
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
468
    VecGetArray(getVecSolInterior(), &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
469

470
    int c = 0;
471 472
    for (int component = 0; component < nComponents; component++) {
      DOFVector<double> &dv = *(vec.getDOFVector(component));
473

474
      DofMap& d = (*interiorMap)[component].getMap();
475 476 477
      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
478 479
    }

480
    VecRestoreArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
481 482 483 484

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


488 489 490
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    Vec tmp;
491
    if (domainComm.Get_size() == 1)
492
      createLocalVec(*interiorMap, tmp);
493
    else
494
      createVec(*interiorMap, tmp);
495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520

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


521 522
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
523
    matDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
524

525 526
    exitPreconditioner(pcInterior);
    exitSolver(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
527 528 529
  }


530 531
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
532
    vecDestroy();
533 534 535
  }


536 537 538 539 540 541 542 543 544 545 546 547 548 549 550
  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;
551
      Parameters::get("parallel->solver->is block " + boost::lexical_cast<string>(i),
552 553 554
		      blockComponents);
      int nComponents = static_cast<int>(blockComponents.size());

555
      TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
556 557 558 559 560 561 562 563

      // 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());
      }

564
      createFieldSplit(pc, isNames[i].c_str(), blockComponents);
565 566 567 568
    }
  }


569
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc, 
570
						 const char* splitName, 
571 572 573 574
						 vector<int> &components)
  {
    IS is;
    interiorMap->createIndexSet(is, components[0], components.size());
575
    PCFieldSplitSetIS(pc, splitName, is);
576 577
    ISDestroy(&is);
  }
578 579 580 581 582 583 584 585 586
  
  
  void PetscSolverGlobalMatrix::extractVectorComponent(Vec input, int i, Vec *output, int numberOfComponents)
  {
    IS is;
    interiorMap->createIndexSet(is, i, numberOfComponents);
    VecGetSubVector(input, is, output);
    ISDestroy(&is);
  }
587

588 589 590 591 592 593 594 595 596 597
  void PetscSolverGlobalMatrix::extractMatrixComponent(Mat input, int startRow, int numberOfRows, int startCol, int numberOfCols, Mat *output)
  {
    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);
  }
  
598

599
  void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
600
  {
601
    KSPCreate(domainComm, &ksp);
602 603 604 605 606
    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());
607
    MatSetOptionsPrefix(getMatInterior(), kspPrefix.c_str());
608 609 610 611 612 613 614 615
    KSPSetFromOptions(ksp);

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


616
  void PetscSolverGlobalMatrix::exitSolver(KSP &ksp)
617 618 619 620 621
  {
    KSPDestroy(&ksp);
  }  


622 623 624 625 626 627 628 629
  void PetscSolverGlobalMatrix::initPreconditioner(PC pc)
  {
    PCSetFromOptions(pc);
    createFieldSplit(pc);
  }

 
  void PetscSolverGlobalMatrix::exitPreconditioner(PC pc)
630
  { }
631 632


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

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

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

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

    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;

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

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

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

      int globalRowDof = rowMultiIndex.global;
675

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

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

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

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

	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.
696
	  MultiIndex colMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
697
	  if ((*interiorMap)[colComp].find(col(*icursor), colMultiIndex) == false)
698 699 700
	    continue;

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

	  // 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
712 713
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
714 715 716 717 718
	  } else {
	    // === Row index is not periodic, but column index is. ===

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

Thomas Witkowski's avatar
Thomas Witkowski committed
722 723 724 725 726 727 728 729 730
	    // 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;
731
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
732
	    for (unsigned int i = 0; i < newCols.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
733
	      cols.push_back(interiorMap->getMatIndex(colComp, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
734 735 736 737 738
	      values.push_back(scaledValue);	      
	    }
	  }
	}

739
  	MatSetValues(getMatInterior(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
740
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
741 742 743 744 745 746 747 748 749 750 751 752 753 754
      } 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
755
	  int globalColDof = (*interiorMap)[colComp][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
756 757 758 759 760 761 762 763 764

	  // 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;
765 766 767 768
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
769 770 771 772 773 774 775 776 777 778

	  // 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
779 780
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
781 782 783

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

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

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


	// === 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;
802

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


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

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

822
    ParallelDofMapping *rowCoarseSpace = 
823
      (coarseSpaceMap.size() ? coarseSpaceMap[rowComp] : nullptr);
Thomas Witkowski's avatar
Thomas Witkowski committed
824 825

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

Thomas Witkowski's avatar
Thomas Witkowski committed
827 828 829
    // 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
830 831
      
      DegreeOfFreedom dof = dofIt.getDOFIndex();
832

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

836 837
      bool isCoarseDof = isCoarseSpace(rowComp, dof);

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

845
      if (isCoarseDof) {
846
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("vecCoarse not set! Should not happen!\n");
847

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
871 872 873
	  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
874
	    int mappedIndex = interiorMap->getMatIndex(rowComp, mappedDof);
875

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

886 887 888 889 890 891 892 893

  PetscSolver* PetscSolverGlobalMatrix::createSubSolver(int component,
							string kspPrefix)
  {
    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
  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp)
  {
    MatNullSpace matNullSpace;
947
    MatNullSpaceCreate(domainComm, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
948 949 950
    KSPSetNullSpace(ksp, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
  }
951

952
} } // end namespace Parallel, AMDiS