PetscSolverGlobalMatrix.cc 26.6 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1 2 3 4 5 6 7 8 9 10 11
//
// 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.

12 13
#include <mpi.h>
#include "DirichletBC.h"
14
#include "DOFVector.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
15 16 17 18 19 20
#include "parallel/PetscSolverGlobalMatrix.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"

namespace AMDiS {

21 22 23 24 25 26 27 28 29 30 31

  PetscSolverGlobalMatrix::PetscSolverGlobalMatrix(string name)
    : PetscSolver(name),
      zeroStartVector(false),
      printMatInfo(false)
  {
    Parameters::get("parallel->use zero start vector", zeroStartVector);
    Parameters::get("parallel->print matrix info", printMatInfo);
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
32
  void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
33 34 35
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
36 37 38 39 40 41
    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();

42
    createMatVec(*seqMat);
43

Thomas Witkowski's avatar
Thomas Witkowski committed
44
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
45
      fillPetscMatrixWithCoarseSpace(seqMat);
46 47
      return;
    }
48
    
49 50
    // === Create PETSc vector (solution and a temporary vector). ===

Thomas Witkowski's avatar
Thomas Witkowski committed
51 52 53 54 55 56
#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
57
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
58 59
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
60 61
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
62 63 64 65 66

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

67
    matAssembly();
68

69 70
    if (printMatInfo) {
      MatInfo matInfo;
71
      MatGetInfo(getMatInterior(), MAT_GLOBAL_SUM, &matInfo);
72 73 74 75 76 77 78
      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));
    }
79 80


81
    // === Init PETSc solver and preconditioner objects. ===
82

83 84
    initSolver(kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
85
    initPreconditioner(pcInterior);
86

87

Thomas Witkowski's avatar
Thomas Witkowski committed
88 89 90
#if (DEBUG != 0)
    MSG("Fill petsc matrix 3 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif
91 92


93 94 95 96 97 98 99 100 101 102 103
    // === 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);
    }
104 105 106
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
107
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *seqMat)
108 109
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
110 111

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

115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
    // === 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);

131 132 133
    bool localMatrix = 
      (meshDistributor->getMeshLevelData().getMpiComm(meshLevel) == MPI::COMM_SELF);

134 135 136
    // === Traverse all sequentially created matrices and add the values to ===
    // === the global PETSc matrices.                                       ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
137
    int nComponents = seqMat->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
138 139 140 141 142
    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
	DOFMatrix* dofMat = (*seqMat)[rowComponent][colComponent];

	if (!dofMat)
143 144
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
145 146
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
147

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

Thomas Witkowski's avatar
Thomas Witkowski committed
150 151
	traits::col<Matrix>::type col(dofMat->getBaseMatrix());
	traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
152 153
	
	// Traverse all rows.
Thomas Witkowski's avatar
Thomas Witkowski committed
154 155
	for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()), 
	       cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
156

Thomas Witkowski's avatar
Thomas Witkowski committed
157
	  bool isRowCoarse = isCoarseSpace(rowComponent, *cursor);
Thomas Witkowski's avatar
Thomas Witkowski committed
158 159 160 161 162 163

	  // For the case, this is a dirichlet row we have to check whether the 
	  // rank is also owner of this row DOF.
	  if (dirichletRows.count(*cursor)) {
	    if ((!isRowCoarse && !(*interiorMap)[rowComponent].isRankDof(*cursor)) ||
		(isRowCoarse && !(*rowCoarseSpace)[rowComponent].isRankDof(*cursor)))
164
	      continue;    
Thomas Witkowski's avatar
Thomas Witkowski committed
165
	  }
166 167 168 169 170 171 172 173 174 175
  
	  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
176
	    bool isColCoarse = isCoarseSpace(colComponent, col(*icursor));
177

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
178
	    if (isColCoarse == false)
179
	      if ((*interiorMap)[colComponent].isSet(col(*icursor)) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
180 181
		continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
182 183 184
	    if (isColCoarse == isRowCoarse) {
	      cols.push_back(col(*icursor));
	      values.push_back(value(*icursor));
185
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
186 187
	      colsOther.push_back(col(*icursor));
	      valuesOther.push_back(value(*icursor));
188 189 190 191 192 193
	    }
	  }  // for each nnz in row


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

Thomas Witkowski's avatar
Thomas Witkowski committed
194
	  if (isRowCoarse) {
Thomas Witkowski's avatar
Thomas Witkowski committed
195 196
	    for (unsigned int i = 0; i < cols.size(); i++)
	      cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]);
197

Thomas Witkowski's avatar
Thomas Witkowski committed
198
	    int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, *cursor);
Thomas Witkowski's avatar
Thomas Witkowski committed
199
	    MatSetValues(getMatCoarseByComponent(rowComponent, colComponent),
200 201
	    		 1, &rowIndex, cols.size(),
	    		 &(cols[0]), &(values[0]), ADD_VALUES);
202 203

	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
204 205 206
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  interiorMap->getMatIndex(colComponent, colsOther[i]) + 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
207
		  rStartInterior;	      
208

Thomas Witkowski's avatar
Thomas Witkowski committed
209
	      MatSetValues(getMatCoarseInteriorByComponent(rowComponent), 
210 211
	       		   1, &rowIndex, colsOther.size(),
 	       		   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
212 213
	    }
	  } else {
214
	    if ((*interiorMap)[rowComponent].isSet(*cursor) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
215 216
	      continue;

217
	    int localRowIndex = 
218
	      (localMatrix ? 
Thomas Witkowski's avatar
Thomas Witkowski committed
219 220
	       interiorMap->getLocalMatIndex(rowComponent, *cursor) :
	       interiorMap->getMatIndex(rowComponent, *cursor));
221

Thomas Witkowski's avatar
Thomas Witkowski committed
222
	    for (unsigned int i = 0; i < cols.size(); i++) {
223
	      if (localMatrix)
Thomas Witkowski's avatar
Thomas Witkowski committed
224
		cols[i] = interiorMap->getLocalMatIndex(colComponent, cols[i]);
225
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
226
		cols[i] = interiorMap->getMatIndex(colComponent, cols[i]);
227 228
	    }
	    
229 230
	    MatSetValues(getMatInterior(), 1, &localRowIndex, cols.size(),
  	     		 &(cols[0]), &(values[0]), ADD_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
231
	    
232
	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
233 234 235
	      int globalRowIndex = 
		interiorMap->getMatIndex(rowComponent, *cursor) + rStartInterior;
      
Thomas Witkowski's avatar
Thomas Witkowski committed
236 237 238
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  colCoarseSpace->getMatIndex(colComponent, colsOther[i]);
239

240 241 242
	      MatSetValues(getMatInteriorCoarseByComponent(colComponent), 
			   1, &globalRowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
243 244 245 246 247 248
	    }
	  }
	} 
      }
    }

249
    matAssembly();
250

251 252
    // === Create solver for the non primal (thus local) variables. ===

253
    KSPCreate(domainComm, &kspInterior);
254
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(),
255
		    SAME_NONZERO_PATTERN);
256 257 258
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
259 260
    if (isSymmetric) {
      PCSetType(pcInterior, PCCHOLESKY);
261
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
262 263
    } else {
      PCSetType(pcInterior, PCLU);
264
      if (localMatrix)
265 266 267 268
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
      else
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    }
269
    KSPSetFromOptions(kspInterior);
270 271 272
  }


273 274 275
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
276

277
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
278
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
279
    
Thomas Witkowski's avatar
Thomas Witkowski committed
280
    // === Transfer values from DOF vector to the PETSc vector. === 
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
282
      for (int i = 0; i < vec->getSize(); i++)
283
	setDofVector(getVecRhsInterior(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
284
		     getVecRhsCoarseByComponent(i), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
285 286
    } else {
      for (int i = 0; i < vec->getSize(); i++)
287
	setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
288 289
    }

290
    vecRhsAssembly();
291 292 293 294 295 296 297 298 299 300 301 302

    // === 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
303 304 305
  }


306 307
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
308 309 310 311 312 313 314
  {
    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
315 316
      TEST_EXIT(coarseSpaceMap.size() == 0)("Not yet supported!\n");

317
      VecSet(getVecSolInterior(), 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
318 319
      
      for (int i = 0; i < nComponents; i++)
320
	setDofVector(getVecSolInterior(), vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
321

322
      vecSolAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
323 324
    }

325

Thomas Witkowski's avatar
Thomas Witkowski committed
326 327
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
328 329 330
    if (nullspace.size() > 0 || 
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
331
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
332

333 334 335 336 337 338 339 340 341 342
      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
343
      if (nullspace.size() > 0) {
344
	VecDuplicate(getVecSolInterior(), &nullspaceBasis);
345
	setDofVector(nullspaceBasis, *(nullspace[0]), true);
Thomas Witkowski's avatar
Thomas Witkowski committed
346 347 348
	
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
349 350

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
351
	
352
	MatNullSpaceCreate(domainComm, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
Thomas Witkowski's avatar
Thomas Witkowski committed
353 354
			   1, &nullspaceBasis, &matNullspace);

355
	MatMult(getMatInterior(), nullspaceBasis, getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
356
	PetscReal n;
357
	VecNorm(getVecSolInterior(), NORM_2, &n);
Thomas Witkowski's avatar
Thomas Witkowski committed
358 359
	MSG("NORM IS: %e\n", n);
      } else {
360
	MatNullSpaceCreate(domainComm, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
361
      }
362

Thomas Witkowski's avatar
Thomas Witkowski committed
363 364 365
      MSG("NULLSPACE IS NOT REMOVED!\n");
      // MatSetNullSpace(getMatInterior(), matNullspace);
      // KSPSetNullSpace(kspInterior, matNullspace);
366

Thomas Witkowski's avatar
Thomas Witkowski committed
367
      // === Remove null space, if requested. ===
368

Thomas Witkowski's avatar
Thomas Witkowski committed
369
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
370
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
371 372

	MSG("Remove nullspace from rhs vector.\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
373
	
374
	MatNullSpaceRemove(matNullspace, getVecRhsInterior(), PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
375 376 377 378
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
379 380
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
381
    // PETSc.
382
    solve(getVecRhsInterior(), getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
383

384 385

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
386
      MatNullSpaceDestroy(&matNullspace);
387 388 389 390
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
391 392
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
393
    VecGetArray(getVecSolInterior(), &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
394

395
    int c = 0;
396 397
    for (int component = 0; component < nComponents; component++) {
      DOFVector<double> &dv = *(vec.getDOFVector(component));
398

399
      DofMap& d = (*interiorMap)[component].getMap();
400 401 402
      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
403 404
    }

405
    VecRestoreArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
406 407 408 409

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


413 414 415 416 417
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()");

    Vec tmp;
418
    if (domainComm.Get_size() == 1)
Thomas Witkowski's avatar
Thomas Witkowski committed
419
      interiorMap->createLocalVec(tmp);
420
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
421
      interiorMap->createVec(tmp);
422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447

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


448 449 450 451
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");

452
    matDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
453

454 455 456
    exitPreconditioner(pcInterior);

    exitSolver(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
457 458 459
  }


460 461 462 463
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()");

464
    vecDestroy();
465 466 467
  }


468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486
  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());

487
      TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
488 489 490 491 492 493 494 495

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

496
      createFieldSplit(pc, isNames[i].c_str(), blockComponents);
497 498 499 500
    }
  }


501
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc, 
502
						 const char* splitName, 
503 504 505 506 507 508
						 vector<int> &components)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createFieldSplit()");

    IS is;
    interiorMap->createIndexSet(is, components[0], components.size());
509
    PCFieldSplitSetIS(pc, splitName, is);
510 511 512 513
    ISDestroy(&is);
  }


514
  void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
515 516 517
  {
    FUNCNAME("PetscSolverGlobalMatrix::initSolver()");

518
    KSPCreate(domainComm, &ksp);
519 520 521 522 523 524 525 526 527 528 529 530 531
    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);
  }  


532
  void PetscSolverGlobalMatrix::exitSolver(KSP &ksp)
533 534 535 536 537 538 539
  {
    FUNCNAME("PetscSolverGlobalMatrix::exitSolver()");

    KSPDestroy(&ksp);
  }  


540 541 542 543 544 545 546 547 548 549 550 551 552 553 554
  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
555
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
Thomas Witkowski's avatar
Thomas Witkowski committed
556
					     int rowComp, int colComp)
Thomas Witkowski's avatar
Thomas Witkowski committed
557 558 559
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
560
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
561 562

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
563
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
564 565
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
566 567
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
568 569 570 571 572 573 574 575 576 577 578

    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;

579 580
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
581 582
    std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows();
      
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
583 584
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
585 586

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
589 590
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
591
      // Global index of the current row DOF.
592
      MultiIndex rowMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
593
      if ((*interiorMap)[rowComp].find(*cursor, rowMultiIndex) == false)
594 595 596
	continue;

      int globalRowDof = rowMultiIndex.global;
597

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

Thomas Witkowski's avatar
Thomas Witkowski committed
601 602 603 604
      // Dirichlet rows can be set only be the owner ranks.
      if (dirichletRows.count(*cursor) && !((*interiorMap)[rowComp].isRankDof(*cursor)))
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
605 606 607
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

608
	// Get PETSc's mat row index.
Thomas Witkowski's avatar
Thomas Witkowski committed
609
	int rowIndex = interiorMap->getMatIndex(rowComp, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
610 611 612 613 614 615 616 617

	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.
618
	  MultiIndex colMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
619
	  if ((*interiorMap)[colComp].find(col(*icursor), colMultiIndex) == false)
620 621 622
	    continue;

	  int globalColDof = colMultiIndex.global;
Thomas Witkowski's avatar
Thomas Witkowski committed
623
	  // Test if the current col dof is a periodic dof.
624
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
625
	  // Get PETSc's mat col index.
Thomas Witkowski's avatar
Thomas Witkowski committed
626
	  int colIndex = interiorMap->getMatIndex(colComp, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
627 628 629 630 631 632 633

	  // 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
634 635
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
636 637 638 639 640
	  } else {
	    // === Row index is not periodic, but column index is. ===

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

Thomas Witkowski's avatar
Thomas Witkowski committed
644 645 646 647 648 649 650 651 652
	    // 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;
653
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
654
	    for (unsigned int i = 0; i < newCols.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
655
	      cols.push_back(interiorMap->getMatIndex(colComp, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
656 657 658 659 660
	      values.push_back(scaledValue);	      
	    }
	  }
	}

661
  	MatSetValues(getMatInterior(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
662
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
663 664 665 666 667 668 669 670 671 672 673 674 675 676
      } 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
677
	  int globalColDof = (*interiorMap)[colComp][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
678 679 680 681 682 683 684 685 686

	  // 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;
687 688 689 690
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
691 692 693 694 695 696 697 698 699 700

	  // 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
701 702
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
703 704 705

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

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

710 711
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
712 713 714 715 716 717 718 719 720 721 722 723
	  }
	}


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

725
	  MatSetValues(getMatInterior(), 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
726 727 728 729 730 731 732
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
733 734
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
735
					     DOFVector<double>* vec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
736
					     int rowComp, 
737
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
738 739 740
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

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

744
    ParallelDofMapping *rowCoarseSpace = 
Thomas Witkowski's avatar
Thomas Witkowski committed
745 746 747
      (coarseSpaceMap.size() ? coarseSpaceMap[rowComp] : NULL);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
749 750 751
    // 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
752 753
      
      DegreeOfFreedom dof = dofIt.getDOFIndex();
754

Thomas Witkowski's avatar
Thomas Witkowski committed
755 756
      if (rankOnly && !(*interiorMap)[rowComp].isRankDof(dof))
	continue;
757

758 759
      bool isCoarseDof = isCoarseSpace(rowComp, dof);

Thomas Witkowski's avatar
Thomas Witkowski committed
760
      // Dirichlet rows can be set only be the owner ranks.
761 762 763 764 765
      if (dirichletValues.count(dof)) {
	if ((!isCoarseDof && !((*interiorMap)[rowComp].isRankDof(dof))) ||
	    (isCoarseDof && !((*rowCoarseSpace)[rowComp].isRankDof(dof))))
	  continue;
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
766

767
      if (isCoarseDof) {
768 769
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
770
	int index = rowCoarseSpace->getMatIndex(rowComp, dof);
771 772
	VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES);
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
773
	if ((*interiorMap)[rowComp].isSet(dof) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
774 775
	  continue;

776
	// Calculate global row index of the DOF.
Thomas Witkowski's avatar
Thomas Witkowski committed
777
	DegreeOfFreedom globalRowDof = (*interiorMap)[rowComp][dof].global;
778 779 780 781 782
	
	// Get PETSc's mat index of the row DOF.
	int index = 0;
	if (interiorMap->isMatIndexFromGlobal())
	  index = 
Thomas Witkowski's avatar
Thomas Witkowski committed
783
	    interiorMap->getMatIndex(rowComp, globalRowDof) + rStartInterior;
784 785
	else
	  index =
Thomas Witkowski's avatar
Thomas Witkowski committed
786
	    interiorMap->getMatIndex(rowComp, dof) + rStartInterior;
787

788
	if (perMap.isPeriodic(feSpace, globalRowDof)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
789 790 791
	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
	  double value = *dofIt / (perAsc.size() + 1.0);
	  VecSetValue(vecInterior, index, value, ADD_VALUES);
792

Thomas Witkowski's avatar
Thomas Witkowski committed
793 794 795
	  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
796
	    int mappedIndex = interiorMap->getMatIndex(rowComp, mappedDof);
797

Thomas Witkowski's avatar
Thomas Witkowski committed
798
	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
799 800 801
	  }	  
	} else {	  
	  // The DOF index is not periodic.
Thomas Witkowski's avatar
Thomas Witkowski committed
802 803
	  VecSetValue(vecInterior, index, *dofIt, ADD_VALUES);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
804 805 806 807
      }
    }
  }

808 809 810 811 812 813 814 815 816 817

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

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

    PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
818
    subSolver->setKspPrefix(kspPrefix);
819
    subSolver->setMeshDistributor(meshDistributor, 0);
820 821 822 823 824 825 826 827 828
    subSolver->init(fe, fe);

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

    return subSolver;
  }

829 830 831 832

  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp,
						     int constFeSpace,
						     bool test)
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855
  {
    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);
856
    }
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
857 858 859
    
    
    MatNullSpace matNullSpace;
860
    MatNullSpaceCreate(domainComm, PETSC_FALSE, 1, &nullSpaceBasis, &matNullSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
861 862
    KSPSetNullSpace(ksp, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
863 864

    VecDestroy(&nullSpaceBasis);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
865 866
  }

867

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
868 869 870 871 872
  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp)
  {
    FUNCNAME("PetscSolverGlobalMatrix::setConstantNullSpace()");

    MatNullSpace matNullSpace;
873
    MatNullSpaceCreate(domainComm, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
874 875 876
    KSPSetNullSpace(ksp, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
  }
877

Thomas Witkowski's avatar
Thomas Witkowski committed
878
}