PetscSolverGlobalMatrix.cc 27.4 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
    ISDestroy(&is);
  }
512 513 514 515 516 517 518 519 520 521
  
  
  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);
  }
522

523 524 525 526 527 528 529 530 531 532 533
  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);
  }
  
534

535
  void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
536 537 538
  {
    FUNCNAME("PetscSolverGlobalMatrix::initSolver()");

539
    KSPCreate(domainComm, &ksp);
540 541 542 543 544 545 546 547 548 549 550 551 552
    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);
  }  


553
  void PetscSolverGlobalMatrix::exitSolver(KSP &ksp)
554 555 556 557 558 559 560
  {
    FUNCNAME("PetscSolverGlobalMatrix::exitSolver()");

    KSPDestroy(&ksp);
  }  


561 562 563 564 565 566 567 568 569 570 571 572 573 574 575
  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
576
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
Thomas Witkowski's avatar
Thomas Witkowski committed
577
					     int rowComp, int colComp)
Thomas Witkowski's avatar
Thomas Witkowski committed
578 579 580
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
581
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
582 583

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
584
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
585 586
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
587 588
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
589 590 591 592 593 594 595 596 597 598 599

    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;

600 601
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
602 603
    std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows();
      
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
604 605
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
606 607

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
610 611
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
612
      // Global index of the current row DOF.
613
      MultiIndex rowMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
614
      if ((*interiorMap)[rowComp].find(*cursor, rowMultiIndex) == false)
615 616 617
	continue;

      int globalRowDof = rowMultiIndex.global;
618

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

Thomas Witkowski's avatar
Thomas Witkowski committed
622 623 624 625
      // 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
626 627 628
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

629
	// Get PETSc's mat row index.
Thomas Witkowski's avatar
Thomas Witkowski committed
630
	int rowIndex = interiorMap->getMatIndex(rowComp, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
631 632 633 634 635 636 637 638

	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.
639
	  MultiIndex colMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
640
	  if ((*interiorMap)[colComp].find(col(*icursor), colMultiIndex) == false)
641 642 643
	    continue;

	  int globalColDof = colMultiIndex.global;
Thomas Witkowski's avatar
Thomas Witkowski committed
644
	  // Test if the current col dof is a periodic dof.
645
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
646
	  // Get PETSc's mat col index.
Thomas Witkowski's avatar
Thomas Witkowski committed
647
	  int colIndex = interiorMap->getMatIndex(colComp, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
648 649 650 651 652 653 654

	  // 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
655 656
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
657 658 659 660 661
	  } else {
	    // === Row index is not periodic, but column index is. ===

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

Thomas Witkowski's avatar
Thomas Witkowski committed
665 666 667 668 669 670 671 672 673
	    // 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;
674
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
675
	    for (unsigned int i = 0; i < newCols.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
676
	      cols.push_back(interiorMap->getMatIndex(colComp, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
677 678 679 680 681
	      values.push_back(scaledValue);	      
	    }
	  }
	}

682
  	MatSetValues(getMatInterior(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
683
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
684 685 686 687 688 689 690 691 692 693 694 695 696 697
      } 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
698
	  int globalColDof = (*interiorMap)[colComp][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
699 700 701 702 703 704 705 706 707

	  // 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;
708 709 710 711
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
712 713 714 715 716 717 718 719 720 721

	  // 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
722 723
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
724 725 726

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

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

731 732
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
733 734 735 736 737 738 739 740 741 742 743 744
	  }
	}


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

746
	  MatSetValues(getMatInterior(), 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
747 748 749 750 751 752 753
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
754 755
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
756
					     DOFVector<double>* vec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
757
					     int rowComp, 
758
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
759 760 761
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

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

765
    ParallelDofMapping *rowCoarseSpace = 
Thomas Witkowski's avatar
Thomas Witkowski committed
766 767 768
      (coarseSpaceMap.size() ? coarseSpaceMap[rowComp] : NULL);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
770 771 772
    // 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
773 774
      
      DegreeOfFreedom dof = dofIt.getDOFIndex();
775

Thomas Witkowski's avatar
Thomas Witkowski committed
776 777
      if (rankOnly && !(*interiorMap)[rowComp].isRankDof(dof))
	continue;
778

779 780
      bool isCoarseDof = isCoarseSpace(rowComp, dof);

Thomas Witkowski's avatar
Thomas Witkowski committed
781
      // Dirichlet rows can be set only be the owner ranks.
782 783 784 785 786
      if (dirichletValues.count(dof)) {
	if ((!isCoarseDof && !((*interiorMap)[rowComp].isRankDof(dof))) ||
	    (isCoarseDof && !((*rowCoarseSpace)[rowComp].isRankDof(dof))))
	  continue;
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
787

788
      if (isCoarseDof) {
789 790
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
791
	int index = rowCoarseSpace->getMatIndex(rowComp, dof);
792 793
	VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES);
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
794
	if ((*interiorMap)[rowComp].isSet(dof) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
795 796
	  continue;

797
	// Calculate global row index of the DOF.
Thomas Witkowski's avatar
Thomas Witkowski committed
798
	DegreeOfFreedom globalRowDof = (*interiorMap)[rowComp][dof].global;
799 800 801 802 803
	
	// Get PETSc's mat index of the row DOF.
	int index = 0;
	if (interiorMap->isMatIndexFromGlobal())
	  index = 
Thomas Witkowski's avatar
Thomas Witkowski committed
804
	    interiorMap->getMatIndex(rowComp, globalRowDof) + rStartInterior;
805 806
	else
	  index =
Thomas Witkowski's avatar
Thomas Witkowski committed
807
	    interiorMap->getMatIndex(rowComp, dof) + rStartInterior;
808

809
	if (perMap.isPeriodic(feSpace, globalRowDof)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
810 811 812
	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
	  double value = *dofIt / (perAsc.size() + 1.0);
	  VecSetValue(vecInterior, index, value, ADD_VALUES);
813

Thomas Witkowski's avatar
Thomas Witkowski committed
814 815 816
	  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
817
	    int mappedIndex = interiorMap->getMatIndex(rowComp, mappedDof);
818

Thomas Witkowski's avatar
Thomas Witkowski committed
819
	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
820 821 822
	  }	  
	} else {	  
	  // The DOF index is not periodic.
Thomas Witkowski's avatar
Thomas Witkowski committed
823 824
	  VecSetValue(vecInterior, index, *dofIt, ADD_VALUES);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
825 826 827 828
      }
    }
  }

829 830 831 832 833 834 835 836 837 838

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

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

    PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
839
    subSolver->setKspPrefix(kspPrefix);
840
    subSolver->setMeshDistributor(meshDistributor, 0);
841 842 843 844 845 846 847 848 849
    subSolver->init(fe, fe);

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

    return subSolver;
  }

850 851 852 853

  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp,
						     int constFeSpace,
						     bool test)
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876
  {
    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);
877
    }
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
878 879 880
    
    
    MatNullSpace matNullSpace;
881
    MatNullSpaceCreate(domainComm, PETSC_FALSE, 1, &nullSpaceBasis, &matNullSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
882 883
    KSPSetNullSpace(ksp, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
884 885

    VecDestroy(&nullSpaceBasis);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
886 887
  }

888

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
889 890 891 892 893
  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp)
  {
    FUNCNAME("PetscSolverGlobalMatrix::setConstantNullSpace()");

    MatNullSpace matNullSpace;
894
    MatNullSpaceCreate(domainComm, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
895 896 897
    KSPSetNullSpace(ksp, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
  }
898

Thomas Witkowski's avatar
Thomas Witkowski committed
899
}