PetscSolverGlobalMatrix.cc 26.8 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"
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 {

20
21
22
23
24
25
26
27
28
29
30

  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
31
  void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
32
33
34
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");

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

41
    createMatVec(*seqMat);
42

Thomas Witkowski's avatar
Thomas Witkowski committed
43
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
44
      fillPetscMatrixWithCoarseSpace(seqMat);
45
46
47
      return;
    }

48
49
    // === Create PETSc vector (solution and a temporary vector). ===

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

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

66
    matAssembly();
67

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


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

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

86

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


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


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

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

114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
    // === 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
133
    int nComponents = seqMat->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
134
135
136
137
138
    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
	DOFMatrix* dofMat = (*seqMat)[rowComponent][colComponent];

	if (!dofMat)
139
140
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
141
142
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
143

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
153
	  bool isRowCoarse = isCoarseSpace(rowComponent, *cursor);
Thomas Witkowski's avatar
Thomas Witkowski committed
154
155
156
157
158
159

	  // 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)))
160
	      continue;    
Thomas Witkowski's avatar
Thomas Witkowski committed
161
	  }
162
163
164
165
166
167
168
169
170
171
  
	  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
172
	    bool isColCoarse = isCoarseSpace(colComponent, col(*icursor));
173

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
174
	    if (isColCoarse == false)
175
	      if ((*interiorMap)[colComponent].isSet(col(*icursor)) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
176
177
		continue;

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


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

Thomas Witkowski's avatar
Thomas Witkowski committed
190
	  if (isRowCoarse) {
Thomas Witkowski's avatar
Thomas Witkowski committed
191
192
	    for (unsigned int i = 0; i < cols.size(); i++)
	      cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]);
193

Thomas Witkowski's avatar
Thomas Witkowski committed
194
	    int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, *cursor);
Thomas Witkowski's avatar
Thomas Witkowski committed
195
	    MatSetValues(getMatCoarseByComponent(rowComponent, colComponent),
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
196
			 1, &rowIndex, cols.size(),
197
198
199
			 &(cols[0]), &(values[0]), ADD_VALUES);

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

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

213
	    int localRowIndex = 
Thomas Witkowski's avatar
Thomas Witkowski committed
214
215
216
	      (subdomainLevel == 0 ? 
	       interiorMap->getLocalMatIndex(rowComponent, *cursor) :
	       interiorMap->getMatIndex(rowComponent, *cursor));
217

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

Thomas Witkowski's avatar
Thomas Witkowski committed
236
  	      MatSetValues(getMatInteriorCoarseByComponent(colComponent), 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
237
			   1, &globalRowIndex, colsOther.size(),
238
239
240
241
242
243
244
  			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  }
	} 
      }
    }

245
    matAssembly();
246

247
248
249
    // === Create solver for the non primal (thus local) variables. ===

    KSPCreate(mpiCommLocal, &kspInterior);
250
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(),
251
		    SAME_NONZERO_PATTERN);
252
253
254
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
255
256
    if (isSymmetric) {
      PCSetType(pcInterior, PCCHOLESKY);
257
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
258
259
260
261
262
263
264
    } else {
      PCSetType(pcInterior, PCLU);
      if (subdomainLevel == 0)
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
      else
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    }
265
    KSPSetFromOptions(kspInterior);
266
267
268
  }


269
270
271
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
272

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

286
    vecRhsAssembly();
287
288
289
290
291
292
293
294
295
296
297
298

    // === 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
299
300
301
  }


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

313
      VecSet(getVecSolInterior(), 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
314
315
      
      for (int i = 0; i < nComponents; i++)
316
	setDofVector(getVecSolInterior(), vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
317

318
      vecSolAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
319
320
    }

321

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

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

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
349
350
	
	MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
			   1, &nullspaceBasis, &matNullspace);

351
	MatMult(getMatInterior(), nullspaceBasis, getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
352
	PetscReal n;
353
	VecNorm(getVecSolInterior(), NORM_2, &n);
Thomas Witkowski's avatar
Thomas Witkowski committed
354
355
356
357
	MSG("NORM IS: %e\n", n);
      } else {
	MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
      }
358

Thomas Witkowski's avatar
Thomas Witkowski committed
359
360
361
      MSG("NULLSPACE IS NOT REMOVED!\n");
      // MatSetNullSpace(getMatInterior(), matNullspace);
      // KSPSetNullSpace(kspInterior, matNullspace);
362

Thomas Witkowski's avatar
Thomas Witkowski committed
363
      // === Remove null space, if requested. ===
364

Thomas Witkowski's avatar
Thomas Witkowski committed
365
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
366
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
367
368

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

Thomas Witkowski's avatar
Thomas Witkowski committed
377
    // PETSc.
378
    solve(getVecRhsInterior(), getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
379

380
381

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
382
      MatNullSpaceDestroy(&matNullspace);
383
384
385
386
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
387
388
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
389
    VecGetArray(getVecSolInterior(), &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
390

391
    int c = 0;
392
393
    for (int component = 0; component < nComponents; component++) {
      DOFVector<double> &dv = *(vec.getDOFVector(component));
394

395
      DofMap& d = (*interiorMap)[component].getMap();
396
397
398
      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
399
400
    }

401
    VecRestoreArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
402
403
404
405

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


409
410
411
412
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()");

413
414
415
    double wtime = MPI::Wtime();
    double t0 = 0.0, t1 = 0.0;

416
417
    Vec tmp;
    if (mpiCommLocal.Get_size() == 1)
Thomas Witkowski's avatar
Thomas Witkowski committed
418
      interiorMap->createLocalVec(tmp);
419
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
420
      interiorMap->createVec(tmp);
421
422
423
424
425
426
427
428
429
430
431

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

432
433
434
    t0 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
435
    KSPSolve(kspInterior, tmp, tmp);
436
437
438
    t1 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
439
440
441
442
443
444
445
446
447
448
449

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

452
    //    MSG("TIMEING: %.5f %.5f\n", t0, t1);
453
454
455
  }


456
457
458
459
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");

460
    matDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
461

462
463
464
    exitPreconditioner(pcInterior);

    exitSolver(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
465
466
467
  }


468
469
470
471
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()");

472
    vecDestroy();
473
474
475
  }


476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
  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());

495
      TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
496
497
498
499
500
501
502
503

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

504
      createFieldSplit(pc, isNames[i].c_str(), blockComponents);
505
506
507
508
    }
  }


509
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc, 
510
						 const char* splitName, 
511
512
513
514
515
516
						 vector<int> &components)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createFieldSplit()");

    IS is;
    interiorMap->createIndexSet(is, components[0], components.size());
517
    PCFieldSplitSetIS(pc, splitName, is);
518
519
520
521
    ISDestroy(&is);
  }


522
  void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
  {
    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);
  }  


548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
  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
563
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
Thomas Witkowski's avatar
Thomas Witkowski committed
564
					     int rowComp, int colComp)
Thomas Witkowski's avatar
Thomas Witkowski committed
565
566
567
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
568
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
569
570

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
571
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
572
573
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
574
575
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
576
577
578
579
580
581
582
583
584
585
586

    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;

587
588
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
589
590
    std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows();
      
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
591
592
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
593
594

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
597
598
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
599
      // Global index of the current row DOF.
600
      MultiIndex rowMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
601
      if ((*interiorMap)[rowComp].find(*cursor, rowMultiIndex) == false)
602
603
604
	continue;

      int globalRowDof = rowMultiIndex.global;
605

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

Thomas Witkowski's avatar
Thomas Witkowski committed
609
610
611
612
      // 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
613
614
615
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

616
	// Get PETSc's mat row index.
Thomas Witkowski's avatar
Thomas Witkowski committed
617
	int rowIndex = interiorMap->getMatIndex(rowComp, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
618
619
620
621
622
623
624
625

	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.
626
	  MultiIndex colMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
627
	  if ((*interiorMap)[colComp].find(col(*icursor), colMultiIndex) == false)
628
629
630
	    continue;

	  int globalColDof = colMultiIndex.global;
Thomas Witkowski's avatar
Thomas Witkowski committed
631
	  // Test if the current col dof is a periodic dof.
632
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
633
	  // Get PETSc's mat col index.
Thomas Witkowski's avatar
Thomas Witkowski committed
634
	  int colIndex = interiorMap->getMatIndex(colComp, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
635
636
637
638
639
640
641

	  // 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
642
643
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
644
645
646
647
648
	  } else {
	    // === Row index is not periodic, but column index is. ===

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

Thomas Witkowski's avatar
Thomas Witkowski committed
652
653
654
655
656
657
658
659
660
	    // 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;
661
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
662
	    for (unsigned int i = 0; i < newCols.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
663
	      cols.push_back(interiorMap->getMatIndex(colComp, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
666
667
668
	      values.push_back(scaledValue);	      
	    }
	  }
	}

669
  	MatSetValues(getMatInterior(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
670
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
671
672
673
674
675
676
677
678
679
680
681
682
683
684
      } 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
685
	  int globalColDof = (*interiorMap)[colComp][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
686
687
688
689
690
691
692
693
694

	  // 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;
695
696
697
698
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
699
700
701
702
703
704
705
706
707
708

	  // 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
709
710
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
711
712
713

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

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

718
719
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
720
721
722
723
724
725
726
727
728
729
730
731
	  }
	}


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

733
	  MatSetValues(getMatInterior(), 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
734
735
736
737
738
739
740
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
741
742
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
743
					     DOFVector<double>* vec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
744
					     int rowComp, 
745
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
746
747
748
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

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

752
    ParallelDofMapping *rowCoarseSpace = 
Thomas Witkowski's avatar
Thomas Witkowski committed
753
754
755
      (coarseSpaceMap.size() ? coarseSpaceMap[rowComp] : NULL);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
757
758
759
    // 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
760
761
762
763
764
      
      DegreeOfFreedom dof = dofIt.getDOFIndex();
      
      if (rankOnly && !(*interiorMap)[rowComp].isRankDof(dof))
	continue;
765

766
767
      bool isCoarseDof = isCoarseSpace(rowComp, dof);

Thomas Witkowski's avatar
Thomas Witkowski committed
768
      // Dirichlet rows can be set only be the owner ranks.
769
770
771
772
773
      if (dirichletValues.count(dof)) {
	if ((!isCoarseDof && !((*interiorMap)[rowComp].isRankDof(dof))) ||
	    (isCoarseDof && !((*rowCoarseSpace)[rowComp].isRankDof(dof))))
	  continue;
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
774

775
      if (isCoarseDof) {
776
777
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
778
	int index = rowCoarseSpace->getMatIndex(rowComp, dof);
779
780
	VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES);
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
781
	if ((*interiorMap)[rowComp].isSet(dof) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
782
783
	  continue;

784
	// Calculate global row index of the DOF.
Thomas Witkowski's avatar
Thomas Witkowski committed
785
	DegreeOfFreedom globalRowDof = (*interiorMap)[rowComp][dof].global;
786
787
788
789
790
	
	// Get PETSc's mat index of the row DOF.
	int index = 0;
	if (interiorMap->isMatIndexFromGlobal())
	  index = 
Thomas Witkowski's avatar
Thomas Witkowski committed
791
	    interiorMap->getMatIndex(rowComp, globalRowDof) + rStartInterior;
792
793
	else
	  index =
Thomas Witkowski's avatar
Thomas Witkowski committed
794
	    interiorMap->getMatIndex(rowComp, dof) + rStartInterior;
795

796
	if (perMap.isPeriodic(feSpace, globalRowDof)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
797
798
799
	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
	  double value = *dofIt / (perAsc.size() + 1.0);
	  VecSetValue(vecInterior, index, value, ADD_VALUES);
800

Thomas Witkowski's avatar
Thomas Witkowski committed
801
802
803
	  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
804
	    int mappedIndex = interiorMap->getMatIndex(rowComp, mappedDof);
805

Thomas Witkowski's avatar
Thomas Witkowski committed
806
	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
807
808
809
	  }	  
	} else {	  
	  // The DOF index is not periodic.
Thomas Witkowski's avatar
Thomas Witkowski committed
810
811
	  VecSetValue(vecInterior, index, *dofIt, ADD_VALUES);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
812
813
814
815
      }
    }
  }

816
817
818
819
820
821
822
823
824
825

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

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

    PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
826
    subSolver->setKspPrefix(kspPrefix);
827
828
829
830
831
832
833
834
835
836
837
838
    subSolver->setMeshDistributor(meshDistributor,
				  mpiCommGlobal,
				  mpiCommLocal);
    subSolver->init(fe, fe);

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

    return subSolver;
  }

839
840
841
842

  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp,
						     int constFeSpace,
						     bool test)
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
  {
    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);
866
    }
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
867
868
869
870
871
872
873
874
    
    
    MatNullSpace matNullSpace;
    MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, &matNullSpace);
    KSPSetNullSpace(ksp, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
  }

875

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
876
877
878
879
880
881
882
883
884
  void PetscSolverGlobalMatrix::setConstantNullSpace(KSP ksp)
  {
    FUNCNAME("PetscSolverGlobalMatrix::setConstantNullSpace()");

    MatNullSpace matNullSpace;
    MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
    KSPSetNullSpace(ksp, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);
  }
885

Thomas Witkowski's avatar
Thomas Witkowski committed
886
}