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


13
#include "AMDiS.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14
15
16
17
18
19
#include "parallel/PetscSolverGlobalMatrix.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"

namespace AMDiS {

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
20
  void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
21
22
23
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
24
25
26
27
28
29
    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();

30
    createMatVec(*seqMat);
31

Thomas Witkowski's avatar
Thomas Witkowski committed
32
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
33
      fillPetscMatrixWithCoarseSpace(seqMat);
34
35
36
      return;
    }

37
38
    // === Create PETSc vector (solution and a temporary vector). ===

Thomas Witkowski's avatar
Thomas Witkowski committed
39
40
41
42
43
44
#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
45
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
46
47
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
48
49
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
50
51
52
53
54

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

55
    matAssembly();
56
57

    removeDirichletRows(seqMat, dofMap, getMatInterior());
58
59
60
    
    if (printMatInfo) {
      MatInfo matInfo;
61
      MatGetInfo(getMatInterior(), MAT_GLOBAL_SUM, &matInfo);
62
63
64
65
66
67
68
      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));
    }
69
70


71
    // === Init PETSc solver and preconditioner objects. ===
72

73
74
    initSolver(kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
75
    initPreconditioner(pcInterior);
76

77

Thomas Witkowski's avatar
Thomas Witkowski committed
78
79
80
#if (DEBUG != 0)
    MSG("Fill petsc matrix 3 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif
81
82


83
84
85
86
87
88
89
90
91
92
93
    // === 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);
    }
94
95
96
  }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
97
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *seqMat)
98
99
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
100
101

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

105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
    // === 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);

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
124
    int nComponents = seqMat->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
125
126
127
128
129
    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
	DOFMatrix* dofMat = (*seqMat)[rowComponent][colComponent];

	if (!dofMat)
130
131
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
132
133
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
134

Thomas Witkowski's avatar
Thomas Witkowski committed
135
136
	traits::col<Matrix>::type col(dofMat->getBaseMatrix());
	traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
137
138
	
	// Traverse all rows.
Thomas Witkowski's avatar
Thomas Witkowski committed
139
140
	for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()), 
	       cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
141

Thomas Witkowski's avatar
Thomas Witkowski committed
142
	  bool isRowCoarse = isCoarseSpace(rowComponent, *cursor);
143
144
145
146
147
148
149
150
151
152
  
	  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
153
	    bool isColCoarse = isCoarseSpace(colComponent, col(*icursor));
154

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
155
	    if (isColCoarse == false)
156
	      if ((*interiorMap)[colComponent].isSet(col(*icursor)) == false)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
157
158
		continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
159
160
161
	    if (isColCoarse == isRowCoarse) {
	      cols.push_back(col(*icursor));
	      values.push_back(value(*icursor));
162
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
163
164
	      colsOther.push_back(col(*icursor));
	      valuesOther.push_back(value(*icursor));
165
166
167
168
169
170
	    }
	  }  // for each nnz in row


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

Thomas Witkowski's avatar
Thomas Witkowski committed
171
	  if (isRowCoarse) {
Thomas Witkowski's avatar
Thomas Witkowski committed
172
173
	    for (unsigned int i = 0; i < cols.size(); i++)
	      cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]);
174

Thomas Witkowski's avatar
Thomas Witkowski committed
175
	    int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, *cursor);
Thomas Witkowski's avatar
Thomas Witkowski committed
176
	    MatSetValues(getMatCoarseByComponent(rowComponent, colComponent),
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
177
			 1, &rowIndex, cols.size(),
178
179
180
			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
181
182
183
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  interiorMap->getMatIndex(colComponent, colsOther[i]) + 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
184
		  rStartInterior;	      
185

Thomas Witkowski's avatar
Thomas Witkowski committed
186
	      MatSetValues(getMatCoarseInteriorByComponent(rowComponent), 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
187
			   1, &rowIndex, colsOther.size(),
188
189
190
 			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  } else {
191
	    if ((*interiorMap)[rowComponent].isSet(*cursor) == false)
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
192
193
	      continue;

194
	    int localRowIndex = 
Thomas Witkowski's avatar
Thomas Witkowski committed
195
196
197
	      (subdomainLevel == 0 ? 
	       interiorMap->getLocalMatIndex(rowComponent, *cursor) :
	       interiorMap->getMatIndex(rowComponent, *cursor));
198

Thomas Witkowski's avatar
Thomas Witkowski committed
199
	    for (unsigned int i = 0; i < cols.size(); i++) {
200
	      if (subdomainLevel == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
201
		cols[i] = interiorMap->getLocalMatIndex(colComponent, cols[i]);
202
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
203
		cols[i] = interiorMap->getMatIndex(colComponent, cols[i]);
204
205
	    }
	    
206
  	    MatSetValues(getMatInterior(), 1, &localRowIndex, cols.size(),
207
  			 &(cols[0]), &(values[0]), ADD_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
208
	    
209
	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
210
211
212
	      int globalRowIndex = 
		interiorMap->getMatIndex(rowComponent, *cursor) + rStartInterior;
      
Thomas Witkowski's avatar
Thomas Witkowski committed
213
214
215
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  colCoarseSpace->getMatIndex(colComponent, colsOther[i]);
216

Thomas Witkowski's avatar
Thomas Witkowski committed
217
  	      MatSetValues(getMatInteriorCoarseByComponent(colComponent), 
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
218
			   1, &globalRowIndex, colsOther.size(),
219
220
221
222
223
224
225
  			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  }
	} 
      }
    }

226
    matAssembly();
227

228
229
230
    // === Create solver for the non primal (thus local) variables. ===

    KSPCreate(mpiCommLocal, &kspInterior);
231
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(),
232
		    SAME_NONZERO_PATTERN);
233
234
235
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
236
237
    if (isSymmetric) {
      PCSetType(pcInterior, PCCHOLESKY);
238
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
239
240
241
242
243
244
245
    } else {
      PCSetType(pcInterior, PCLU);
      if (subdomainLevel == 0)
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
      else
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    }
246
    KSPSetFromOptions(kspInterior);
247
248
249
  }


250
251
252
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
253

254
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
255
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
256
    
Thomas Witkowski's avatar
Thomas Witkowski committed
257
    // === Transfer values from DOF vector to the PETSc vector. === 
Thomas Witkowski's avatar
Thomas Witkowski committed
258
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
259
      for (int i = 0; i < vec->getSize(); i++)
260
	setDofVector(getVecRhsInterior(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
261
		     getVecRhsCoarseByComponent(i), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
262
263
    } else {
      for (int i = 0; i < vec->getSize(); i++)
264
	setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
265
266
    }

267
    vecRhsAssembly();
268

269
    removeDirichletRows(vec, dofMap, getVecRhsInterior());
270
271
272
273
274
275
276
277
278
279
280
281

    // === 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
282
283
284
  }


285
286
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
287
288
289
290
291
292
293
  {
    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
294
295
      TEST_EXIT(coarseSpaceMap.size() == 0)("Not yet supported!\n");

296
      VecSet(getVecSolInterior(), 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
297
298
      
      for (int i = 0; i < nComponents; i++)
299
	setDofVector(getVecSolInterior(), vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
300

301
      vecSolAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
    }

304

Thomas Witkowski's avatar
Thomas Witkowski committed
305
306
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
307
308
309
    if (nullspace.size() > 0 || 
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
310
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
311

312
313
314
315
316
317
318
319
320
321
      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
322
      if (nullspace.size() > 0) {
323
	VecDuplicate(getVecSolInterior(), &nullspaceBasis);
Thomas Witkowski's avatar
Thomas Witkowski committed
324
325
326
327
328
	for (int i = 0; i < nComponents; i++)
	  setDofVector(nullspaceBasis, nullspace[0]->getDOFVector(i), i, true);
	
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
329
330

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
331
332
333
334
	
	MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
			   1, &nullspaceBasis, &matNullspace);

335
	MatMult(getMatInterior(), nullspaceBasis, getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
336
	PetscReal n;
337
	VecNorm(getVecSolInterior(), NORM_2, &n);
Thomas Witkowski's avatar
Thomas Witkowski committed
338
339
340
341
	MSG("NORM IS: %e\n", n);
      } else {
	MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
      }
342

Thomas Witkowski's avatar
Thomas Witkowski committed
343
344
345
      MSG("NULLSPACE IS NOT REMOVED!\n");
      // MatSetNullSpace(getMatInterior(), matNullspace);
      // KSPSetNullSpace(kspInterior, matNullspace);
346

Thomas Witkowski's avatar
Thomas Witkowski committed
347
      // === Remove null space, if requested. ===
348

Thomas Witkowski's avatar
Thomas Witkowski committed
349
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
350
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
351
352

	MSG("Remove nullspace from rhs vector.\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
353
	
354
	MatNullSpaceRemove(matNullspace, getVecRhsInterior(), PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
355
356
357
358
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
359
360
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
361
    // PETSc.
362
    solve(getVecRhsInterior(), getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
363

364
365

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
366
      MatNullSpaceDestroy(&matNullspace);
367
368
369
370
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
371
372
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
373
    VecGetArray(getVecSolInterior(), &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
374

375
    int c = 0;
376
377
    for (int component = 0; component < nComponents; component++) {
      DOFVector<double> &dv = *(vec.getDOFVector(component));
378

379
      DofMap& d = (*interiorMap)[component].getMap();
380
381
382
      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
383
384
    }

385
    VecRestoreArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
386
387
388
389

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


393
394
395
396
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()");

397
398
399
    double wtime = MPI::Wtime();
    double t0 = 0.0, t1 = 0.0;

400
401
    Vec tmp;
    if (mpiCommLocal.Get_size() == 1)
Thomas Witkowski's avatar
Thomas Witkowski committed
402
      interiorMap->createLocalVec(tmp);
403
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
404
      interiorMap->createVec(tmp);
405
406
407
408
409
410
411
412
413
414
415

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

416
417
418
    t0 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
419
    KSPSolve(kspInterior, tmp, tmp);
420
421
422
    t1 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
423
424
425
426
427
428
429
430
431
432
433

    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);
434
435
    t0 += MPI::Wtime() - wtime;

436
    //    MSG("TIMEING: %.5f %.5f\n", t0, t1);
437
438
439
  }


440
441
442
443
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");

444
    matDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
445

446
447
448
    exitPreconditioner(pcInterior);

    exitSolver(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
449
450
451
  }


452
453
454
455
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()");

456
    vecDestroy();
457
458
459
  }


460
461
462
463
464
465
466
467
468
469
  void PetscSolverGlobalMatrix::removeDirichletRows(Matrix<DOFMatrix*> *seqMat, 
						    ParallelDofMapping &dofMap,
						    Mat mpiMat)
  {
    FUNCNAME("PetscSolverGlobalMatrix::removeDirichletRows()");

    if (!handleDirichletRows)
      return;

    int nComponents = seqMat->getSize();
470
471
472
    vector<int> dirichletRows;
    vector<int> dirichletCols;
    vector<double> dirichletValues;
473
474

    for (int i = 0; i < nComponents; i++) {
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
      bool dirichletRow = false;
      int dirichletMainCol = -1;

      for (int j = 0; j < nComponents; j++) {
	DOFMatrix *dofMat = (*seqMat)[i][j];
	if (!dofMat)
	  continue;

	BoundaryIndexMap& bounds = 
	  const_cast<BoundaryIndexMap&>(dofMat->getBoundaryManager()->getBoundaryConditionMap());

	for (BoundaryIndexMap::iterator bIt = bounds.begin(); bIt != bounds.end(); ++bIt) {
	  if (bIt->second && bIt->second->isDirichlet()) {
	    dirichletRow = true;
	    if ((dynamic_cast<DirichletBC*>(bIt->second))->applyBoundaryCondition()) {
	      TEST_EXIT_DBG(dirichletMainCol == -1)("Should not happen!\n");
	      dirichletMainCol = j;
	    } 
	  }
	}
      }

      if (!dirichletRow)
498
499
	continue;

500
      DOFMatrix *dofMat = (*seqMat)[i][dirichletMainCol];
501
      const FiniteElemSpace *feSpace = dofMat->getRowFeSpace();
502
503
      TEST_EXIT(dofMat->getRowFeSpace() == dofMat->getColFeSpace())
	("I have to think about this scenario! Really possible?\n");
504
505
506
507
508
509

      std::set<DegreeOfFreedom> &dRows = dofMat->getDirichletRows();
      for (std::set<DegreeOfFreedom>::iterator dofIt = dRows.begin();
	   dofIt != dRows.end(); ++dofIt) {
	MultiIndex multiIndex;
	if (dofMap[feSpace].find(*dofIt, multiIndex) &&
510
511
512
513
514
515
516
	    dofMap[feSpace].isRankDof(*dofIt)) {
	  int rowIndex = dofMap.getMatIndex(i, multiIndex.global);
	  int colIndex = dofMap.getMatIndex(dirichletMainCol, multiIndex.global);
	  dirichletRows.push_back(rowIndex);
	  dirichletCols.push_back(colIndex);
	  dirichletValues.push_back(1.0);
	}
517
518
519
      }
    }

520
    MatZeroRows(mpiMat, dirichletRows.size(), &(dirichletRows[0]), 0.0,
521
522
		PETSC_NULL, PETSC_NULL);

523
524
    for (int i = 0; i < static_cast<int>(dirichletRows.size()); i++)
      MatSetValue(mpiMat, dirichletRows[i], dirichletCols[i], dirichletValues[i], INSERT_VALUES);    
525
526
527
    
    MatAssemblyBegin(mpiMat, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mpiMat, MAT_FINAL_ASSEMBLY);
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
  }

  
  void PetscSolverGlobalMatrix::removeDirichletRows(SystemVector *seqVec, 
						    ParallelDofMapping &dofMap,
						    Vec mpiVec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::removeDirichletRows()");

    if (!handleDirichletRows)
      return;

    int nComponents = seqVec->getSize();

    for (int i = 0; i < nComponents; i++) {
      DOFVector<double> *dofVec = seqVec->getDOFVector(i);
      const FiniteElemSpace *feSpace = dofVec->getFeSpace();

      map<DegreeOfFreedom, double>& dValues = dofVec->getDirichletValues();
      for (map<DegreeOfFreedom, double>::iterator dofIt = dValues.begin();
	   dofIt != dValues.end(); ++dofIt) {
	MultiIndex multiIndex;
	if (dofMap[feSpace].find(dofIt->first, multiIndex) &&
	    dofMap[feSpace].isRankDof(dofIt->first))
	  VecSetValue(mpiVec, dofMap.getMatIndex(i, multiIndex.global),
		      dofIt->second, INSERT_VALUES);
      }
    }

    VecAssemblyBegin(mpiVec);
    VecAssemblyEnd(mpiVec);
  }


562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
  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());

581
      TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
582
583
584
585
586
587
588
589

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

590
      createFieldSplit(pc, isNames[i], blockComponents);
591
592
593
594
    }
  }


595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc, 
						 string splitName, 
						 vector<int> &components)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createFieldSplit()");

    IS is;
    interiorMap->createIndexSet(is, components[0], components.size());
    PCFieldSplitSetIS(pc, isNames[i].c_str(), is);
    ISDestroy(&is);
  }


  void PetscSolverGlobalMatrix::initSolver(KSP ksp)
  {
    FUNCNAME("PetscSolverGlobalMatrix::initSolver()");

    KSPCreate(mpiCommGlobal, &ksp);
    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);
  }  


  void PetscSolverGlobalMatrix::exitSolver(KSP ksp)
  {
    FUNCNAME("PetscSolverGlobalMatrix::exitSolver()");

    KSPDestroy(&ksp);
  }  


634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
  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
649
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
650
					     int nRowMat, int nColMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
651
652
653
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
654
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
655
656

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
657
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
658
659
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
660
661
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
662
663
664
665
666
667
668
669
670
671
672

    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;

673
674
675
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
676
677
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
678
679

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

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
682
683
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
684
      // Global index of the current row DOF.
685
      MultiIndex rowMultiIndex;
686
      if ((*interiorMap)[nRowMat].find(*cursor, rowMultiIndex) == false)
687
688
689
	continue;

      int globalRowDof = rowMultiIndex.global;
690

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

Thomas Witkowski's avatar
Thomas Witkowski committed
694
695
696
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

697
	// Get PETSc's mat row index.
698
	int rowIndex = interiorMap->getMatIndex(nRowMat, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
699
700
701
702
703
704
705
706

	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.
707
	  MultiIndex colMultiIndex;
708
	  if ((*interiorMap)[nColMat].find(col(*icursor), colMultiIndex) == false)
709
710
711
	    continue;

	  int globalColDof = colMultiIndex.global;
Thomas Witkowski's avatar
Thomas Witkowski committed
712
	  // Test if the current col dof is a periodic dof.
713
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
714
	  // Get PETSc's mat col index.
715
	  int colIndex = interiorMap->getMatIndex(nColMat, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
718
719
720
721
722

	  // 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
723
724
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
725
726
727
728
729
	  } else {
	    // === Row index is not periodic, but column index is. ===

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

Thomas Witkowski's avatar
Thomas Witkowski committed
733
734
735
736
737
738
739
740
741
	    // 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;
742
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
743
	    for (unsigned int i = 0; i < newCols.size(); i++) {
744
	      cols.push_back(interiorMap->getMatIndex(nColMat, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
745
746
747
748
749
	      values.push_back(scaledValue);	      
	    }
	  }
	}

750
  	MatSetValues(getMatInterior(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
751
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
752
753
754
755
756
757
758
759
760
761
762
763
764
765
      } 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.
766
	  int globalColDof = (*interiorMap)[nColMat][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
767
768
769
770
771
772
773
774
775

	  // 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;
776
777
778
779
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
780
781
782
783
784
785
786
787
788
789

	  // 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
790
791
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
792
793
794

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

795
	  for (unsigned int i = 0; i < entry.size(); i++) {
796
797
	    int rowIdx = interiorMap->getMatIndex(nRowMat, entry[i].first);
	    int colIdx = interiorMap->getMatIndex(nColMat, entry[i].second);
Thomas Witkowski's avatar
Thomas Witkowski committed
798

799
800
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
801
802
803
804
805
806
807
808
809
810
811
812
	  }
	}


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

814
	  MatSetValues(getMatInterior(), 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
815
816
817
818
819
820
821
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
822
823
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
824
					     DOFVector<double>* vec, 
825
					     int nRowVec, 
826
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
827
828
829
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

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

833
834
    ParallelDofMapping *rowCoarseSpace = 
      (coarseSpaceMap.size() ? coarseSpaceMap[nRowVec] : NULL);
835

Thomas Witkowski's avatar
Thomas Witkowski committed
836
837
838
    // Traverse all used DOFs in the dof vector.
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
839

840
      if (rankOnly && !(*interiorMap)[nRowVec].isRankDof(dofIt.getDOFIndex()))
Thomas Witkowski's avatar
Thomas Witkowski committed
841
842
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
843
      if (isCoarseSpace(nRowVec, dofIt.getDOFIndex())) {
844
845
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");

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

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

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

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

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

}