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

20
  void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
Thomas Witkowski's avatar
Thomas Witkowski committed
21
22
23
24
25
26
27
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");

    TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT_DBG(mat)("No DOF matrix defined!\n");

    double wtime = MPI::Wtime();
28
29
30
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
    int nRankRows = meshDistributor->getNumberRankDofs(feSpaces);
    int nOverallRows = meshDistributor->getNumberOverallDofs(feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
31

32
    // === Create PETSc vector (solution and a temporary vector). ===
Thomas Witkowski's avatar
Thomas Witkowski committed
33

34
35
    VecCreateMPI(mpiComm, nRankRows, nOverallRows, &petscSolVec);
    VecCreateMPI(mpiComm, nRankRows, nOverallRows, &petscTmpVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
36

37
38
39
40
41
42
43
44
    int testddd = 1;
    Parameters::get("block size", testddd);
    if (testddd > 1) {
      VecSetBlockSize(petscSolVec, testddd);
      VecSetBlockSize(petscTmpVec, testddd);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
45
    int recvAllValues = 0;
46
47
48
49
    int sendValue = 
      static_cast<int>(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz);
    meshDistributor->getMpiComm().Allreduce(&sendValue, &recvAllValues, 
					    1, MPI_INT, MPI_SUM);
Thomas Witkowski's avatar
Thomas Witkowski committed
50

51
52
53
    recvAllValues = 1;

    if (!d_nnz || recvAllValues != 0 || alwaysCreateNnzStructure) {
54
55
56
57
58

      // Global DOF mapping must only be recomputed, if the NNZ structure has
      // changed (thus, also the mesh has changed).
      createGlobalDofMapping(feSpaces);

Thomas Witkowski's avatar
Thomas Witkowski committed
59
60
61
62
63
64
65
66
67
68
69
70
71
72
      if (d_nnz) {
	delete [] d_nnz;
	d_nnz = NULL;
	delete [] o_nnz;
	o_nnz = NULL;
      }

      createPetscNnzStructure(mat);
      lastMeshNnz = meshDistributor->getLastMeshChangeIndex();
    }


    // === Create PETSc matrix with the computed nnz data structure. ===

73
    MatCreateMPIAIJ(mpiComm, nRankRows, nRankRows, 
74
 		    nOverallRows, nOverallRows,
75
		    0, d_nnz, 0, o_nnz, &petscMatrix);
76

77
78
79
80
81
    if (testddd > 1) {
      MatSetBlockSize(petscMatrix, testddd);
      MSG("MAT SET BLOCK SIZE: %d\n", testddd);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
82
83
84
85
86
87
88
#if (DEBUG != 0)
    MSG("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif

#if (DEBUG != 0)
    int a, b;
    MatGetOwnershipRange(petscMatrix, &a, &b);
89
    TEST_EXIT(a == meshDistributor->getStartDofs(feSpaces))
Thomas Witkowski's avatar
Thomas Witkowski committed
90
      ("Wrong matrix ownership range!\n");
91
    TEST_EXIT(b == meshDistributor->getStartDofs(feSpaces) + nRankRows)
Thomas Witkowski's avatar
Thomas Witkowski committed
92
93
94
95
96
97
      ("Wrong matrix ownership range!\n");
#endif


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

98
    int nComponents = mat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
99
100
101
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
	if ((*mat)[i][j])
102
	  setDofMatrix((*mat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
103
104
105
106
107
108
109
110

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

    MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);

111
    // === Init PETSc solver. ===
112
    KSPCreate(mpiComm, &solver);
113
114
115
116
    KSPGetPC(solver, &pc);
    KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); 
    KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(solver, KSPBCGS);
117
    KSPSetOptionsPrefix(solver, kspPrefix.c_str());
118
119
120
121
122
123
124
    KSPSetFromOptions(solver);
    PCSetFromOptions(pc);

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

125
126
127
128
129
130
131
    MSG("Fill petsc matrix needed %.5f seconds\n", MPI::Wtime() - wtime);
  }


  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
132

133
134
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
    TEST_EXIT_DBG(meshDistributor)("No mesh distributor defined!\n");
135

136
137
138
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(vec);
    int nRankRows = meshDistributor->getNumberRankDofs(feSpaces);
    int nOverallRows = meshDistributor->getNumberOverallDofs(feSpaces);
139

140
    VecCreateMPI(mpiComm, nRankRows, nOverallRows, &petscRhsVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
141

142
143
144
145
146
147
    int testddd = 1;
    Parameters::get("block size", testddd);
    if (testddd > 1)
      VecSetBlockSize(petscRhsVec, testddd);


Thomas Witkowski's avatar
Thomas Witkowski committed
148
    // === Transfer values from DOF vector to the PETSc vector. === 
149
    for (int i = 0; i < vec->getSize(); i++)
150
      setDofVector(petscRhsVec, vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
151
152
153

    VecAssemblyBegin(petscRhsVec);
    VecAssemblyEnd(petscRhsVec);
154
155

    if (removeRhsNullSpace) {
156
      MSG("Remove constant null space from the RHS!\n");
157
      MatNullSpace sp;
158
      MatNullSpaceCreate(mpiComm, PETSC_TRUE, 0, PETSC_NULL, &sp);
159
160
161
      MatNullSpaceRemove(sp, petscRhsVec, PETSC_NULL);
      MatNullSpaceDestroy(&sp);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
162
163
164
  }


165
166
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
167
168
169
170
171
172
173
174
175
176
  {
    FUNCNAME("PetscSolverGlobalMatrix::solvePetscMatrix()");

    int nComponents = vec.getSize();

    // === Set old solution to be initiual guess for PETSc solver. ===
    if (!zeroStartVector) {
      VecSet(petscSolVec, 0.0);
      
      for (int i = 0; i < nComponents; i++)
177
	setDofVector(petscSolVec, vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
180
181
182
183
      
      VecAssemblyBegin(petscSolVec);
      VecAssemblyEnd(petscSolVec);
    }

    // PETSc.
184
185
186
187
188
    PetscErrorCode solverError = KSPSolve(solver, petscRhsVec, petscSolVec);
    if (solverError != 0) {
      AMDiS::finalize();
      exit(-1);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
189
190
191
192
193

    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
    VecGetArray(petscSolVec, &vecPointer);

194
    int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
195
    for (int i = 0; i < nComponents; i++) {
196
      DOFVector<double> &dv = *(vec.getDOFVector(i));
197
198
      const FiniteElemSpace *feSpace = dv.getFeSpace();
      int nRankDofs = meshDistributor->getNumberRankDofs(feSpace);
199

Thomas Witkowski's avatar
Thomas Witkowski committed
200
      for (int j = 0; j < nRankDofs; j++)
201
	dv[meshDistributor->mapLocalToDof(feSpace, j)] = vecPointer[c++]; 
Thomas Witkowski's avatar
Thomas Witkowski committed
202
203
204
205
206
207
208
209
210
211
212
213
    }

    VecRestoreArray(petscSolVec, &vecPointer);


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


    // Print iteration counter and residual norm of the solution.
    printSolutionInfo(adaptInfo);
214
215
216
217
218
219
220
221
222
  }


  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");

    MatDestroy(&petscMatrix);
    KSPDestroy(&solver);
223
224
    VecDestroy(&petscSolVec);
    VecDestroy(&petscTmpVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
225
226
227
  }


228
229
230
231
232
233
234
235
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()");

    VecDestroy(&petscRhsVec);
  }


236
237
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* mat,
					     int nRowMat, int nColMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
238
239
240
241
242
243
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

    TEST_EXIT(mat)("No DOFMatrix!\n");

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
244
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
    typedef DOFMatrix::base_matrix_type Matrix;

    traits::col<Matrix>::type col(mat->getBaseMatrix());
    traits::const_value<Matrix>::type value(mat->getBaseMatrix());

    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;

260
261
262
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();

Thomas Witkowski's avatar
Thomas Witkowski committed
263
264
265
266
267
268
    // === Traverse all rows of the dof matrix and insert row wise the values ===
    // === to the PETSc matrix.                                               ===

    for (cursor_type cursor = begin<row>(mat->getBaseMatrix()), 
	   cend = end<row>(mat->getBaseMatrix()); cursor != cend; ++cursor) {

269
270
271
      const FiniteElemSpace *rowFe = mat->getRowFeSpace();
      const FiniteElemSpace *colFe = mat->getColFeSpace();

Thomas Witkowski's avatar
Thomas Witkowski committed
272
      // Global index of the current row DOF.
273
      int globalRowDof = 
274
	meshDistributor->mapDofToGlobal(rowFe, *cursor);
Thomas Witkowski's avatar
Thomas Witkowski committed
275
      // Test if the current row DOF is a periodic DOF.
276
      bool periodicRow = perMap.isPeriodic(rowFe, globalRowDof);
277

Thomas Witkowski's avatar
Thomas Witkowski committed
278
279
280
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

281
	// Get PETSc's mat row index.
282
	int rowIndex = dofToMatIndex.get(nRowMat, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
283
284
285
286
287
288
289
290

	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.
291
	  int globalColDof = 
292
	    meshDistributor->mapDofToGlobal(colFe, col(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
293
	  // Test if the current col dof is a periodic dof.
294
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
295
	  // Get PETSc's mat col index.
296
	  int colIndex = dofToMatIndex.get(nColMat, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
297
298
299
300
301
302
303
304
305
306
307
308
309
310

	  // 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.
	    cols.push_back(colIndex);
	    values.push_back(value(*icursor));
	  } else {
	    // === Row index is not periodic, but column index is. ===

	    // Create set of all periodic associations of the column index.
	    std::set<int> perAsc;
311
	    std::set<int>& perColAsc = perMap.getAssociations(colFe, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
312
313
	    for (std::set<int>::iterator it = perColAsc.begin(); 
		 it != perColAsc.end(); ++it)
314
	      if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it))
Thomas Witkowski's avatar
Thomas Witkowski committed
315
316
317
318
319
320
321
322
323
324
325
326
327
328
		perAsc.insert(*it);
    
	    // 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;

	    // First, add the original matrix index.
	    newCols.push_back(globalColDof);
329

Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
332
333
334
335
	    // And add all periodic matrix indices.
	    for (std::set<int>::iterator it = perAsc.begin(); 
		 it != perAsc.end(); ++it) {
	      int nCols = static_cast<int>(newCols.size());

	      for (int i = 0; i < nCols; i++) {
336
 		TEST_EXIT_DBG(perMap.isPeriodic(colFe, *it, newCols[i]))
Thomas Witkowski's avatar
Thomas Witkowski committed
337
338
339
 		  ("Wrong periodic DOF associations at boundary %d with DOF %d!\n",
		   *it, newCols[i]);

340
		newCols.push_back(perMap.map(colFe, *it, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
341
342
343
344
	      }
	    }

	    for (unsigned int i = 0; i < newCols.size(); i++) {
345
	      cols.push_back(dofToMatIndex.get(nColMat, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
	      values.push_back(scaledValue);	      
	    }
	  }
	}

	MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), 
		     &(cols[0]), &(values[0]), ADD_VALUES);	
      } 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.
367
	  int globalColDof = 
368
	    meshDistributor->mapDofToGlobal(colFe, col(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
369
370
371
372
373
374
375
376
377
378
379

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

380
381
	  if (perMap.isPeriodic(colFe, globalColDof)) {
	    std::set<int>& perColAsc = perMap.getAssociations(colFe, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
382
383
	    for (std::set<int>::iterator it = perColAsc.begin(); 
		 it != perColAsc.end(); ++it)
384
	      if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it))
Thomas Witkowski's avatar
Thomas Witkowski committed
385
386
387
		perAsc.insert(*it);
	  }

388
	  std::set<int>& perRowAsc = perMap.getAssociations(rowFe, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
	  for (std::set<int>::iterator it = perRowAsc.begin(); 
	       it != perRowAsc.end(); ++it)
391
	    if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it))
Thomas Witkowski's avatar
Thomas Witkowski committed
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
	      perAsc.insert(*it);

	  // 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;
	  
	  // First, add the original entry.
	  entry.push_back(make_pair(globalRowDof, globalColDof));

407
408
	  // Then, traverse the periodic associations of the row and column
	  // indices and create the corresponding entries.
Thomas Witkowski's avatar
Thomas Witkowski committed
409
410
411
412
	  for (std::set<int>::iterator it = perAsc.begin(); it != perAsc.end(); ++it) {
	    int nEntry = static_cast<int>(entry.size());
	    for (int i = 0; i < nEntry; i++) {
	      int perRowDof = 0;
413

414
415
	      if (perMap.isPeriodic(rowFe, *it, entry[i].first))
		perRowDof = perMap.map(rowFe, *it, entry[i].first);
Thomas Witkowski's avatar
Thomas Witkowski committed
416
417
418
419
	      else
		perRowDof = entry[i].first;

	      int perColDof;
420
421
	      if (perMap.isPeriodic(colFe, *it, entry[i].second))
		perColDof = perMap.map(colFe, *it, entry[i].second);
Thomas Witkowski's avatar
Thomas Witkowski committed
422
423
424
425
426
427
428
429
430
431
432
	      else
		perColDof = entry[i].second;	      	      
	      

	      entry.push_back(make_pair(perRowDof, perColDof));
	    }
	  }


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

433
	  for (unsigned int i = 0; i < entry.size(); i++) {
434
435
	    int rowIdx = dofToMatIndex.get(nRowMat, entry[i].first);
	    int colIdx = dofToMatIndex.get(nColMat, entry[i].second);
Thomas Witkowski's avatar
Thomas Witkowski committed
436

437
438
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
439
440
441
442
443
444
445
446
447
448
449
450
	  }
	}


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

Thomas Witkowski's avatar
Thomas Witkowski committed
452
453
454
455
456
457
458
459
	  MatSetValues(petscMatrix, 1, &rowIndex, rowIt->second.size(),
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


460
461
  void PetscSolverGlobalMatrix::setDofVector(Vec& petscVec, 
					     DOFVector<double>* vec, 
462
					     int nRowVec, 
463
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
464
465
466
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

467
    const FiniteElemSpace *feSpace = vec->getFeSpace();
468
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
469

Thomas Witkowski's avatar
Thomas Witkowski committed
470
471
472
    // Traverse all used DOFs in the dof vector.
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
473
      if (rankOnly && !meshDistributor->getIsRankDof(feSpace, dofIt.getDOFIndex()))
Thomas Witkowski's avatar
Thomas Witkowski committed
474
475
476
477
	continue;

      // Calculate global row index of the DOF.
      DegreeOfFreedom globalRowDof = 
478
	meshDistributor->mapDofToGlobal(feSpace, dofIt.getDOFIndex());
479
      // Get PETSc's mat index of the row DOF.
480
      int index = dofToMatIndex.get(nRowVec, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
481

482
483
      if (perMap.isPeriodic(feSpace, globalRowDof)) {
	std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
484
485
486
	double value = *dofIt / (perAsc.size() + 1.0);
	VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);

487
488
489
	for (std::set<int>::iterator perIt = perAsc.begin(); 
	     perIt != perAsc.end(); ++perIt) {
	  int mappedDof = perMap.map(feSpace, *perIt, globalRowDof);
490
	  int mappedIndex = dofToMatIndex.get(nRowVec, mappedDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
	  VecSetValues(petscVec, 1, &mappedIndex, &value, ADD_VALUES);
	}
      } else {
	// The DOF index is not periodic.
	double value = *dofIt;
	VecSetValues(petscVec, 1, &index, &value, ADD_VALUES);
      }
    }
  }


  void PetscSolverGlobalMatrix::createPetscNnzStructure(Matrix<DOFMatrix*> *mat)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createPetscNnzStructure()");

    TEST_EXIT_DBG(!d_nnz)("There is something wrong!\n");
    TEST_EXIT_DBG(!o_nnz)("There is something wrong!\n");

509
510
511
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
    int nRankRows = meshDistributor->getNumberRankDofs(feSpaces);
    int rankStartIndex = meshDistributor->getStartDofs(feSpaces);
Thomas Witkowski's avatar
Thomas Witkowski committed
512
513
514
515
516
517
518
519
520
521
522
523
524
    d_nnz = new int[nRankRows];
    o_nnz = new int[nRankRows];
    for (int i = 0; i < nRankRows; i++) {
      d_nnz[i] = 0;
      o_nnz[i] = 0;
    }

    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 vector<pair<int, int> > MatrixNnzEntry;
    typedef map<int, DofContainer> RankToDofContainer;

525
526
527
528
    // Stores to each rank a list of nnz entries (i.e. pairs of row and column
    // index) that this rank will send to. These nnz entries will be assembled
    // on this rank, but because the row DOFs are not DOFs of this rank they 
    // will be send to the owner of the row DOFs.
Thomas Witkowski's avatar
Thomas Witkowski committed
529
530
531
532
    map<int, MatrixNnzEntry> sendMatrixEntry;

    // Maps to each DOF that must be send to another rank the rank number of the
    // receiving rank.
533
    map<pair<DegreeOfFreedom, int>, int> sendDofToRank;
Thomas Witkowski's avatar
Thomas Witkowski committed
534

535

536
537
538
539
540
541
542
543
544
545
    // First, create for all ranks, to which we send data to, MatrixNnzEntry 
    // object with 0 entries.
    for (unsigned int i = 0; i < feSpaces.size(); i++) {
      for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpaces[i]);
	   !it.end(); it.nextRank()) {
	sendMatrixEntry[it.getRank()].resize(0);
	
	for (; !it.endDofIter(); it.nextDof())
	  sendDofToRank[make_pair(it.getDofIndex(), i)] = it.getRank();
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
546
547
    }

548
    // Create list of ranks from which we receive data from.
Thomas Witkowski's avatar
Thomas Witkowski committed
549
    std::set<int> recvFromRank;
550
551
552
553
554
555
556
    for (unsigned int i = 0; i < feSpaces.size(); i++) 
      for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpaces[i]);
	   !it.end(); it.nextRank())
	recvFromRank.insert(it.getRank());


    // === Traverse matrices to create nnz data. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
557

558
    int nComponents = mat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
559
560
561
562
563
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
 	if (!(*mat)[i][j])
	  continue;

564
565
566
567
568
	TEST_EXIT_DBG((*mat)[i][j]->getRowFeSpace() == feSpaces[i])
	  ("Should not happen!\n");
	TEST_EXIT_DBG((*mat)[i][j]->getColFeSpace() == feSpaces[j])
	  ("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
569
570
571
572
573
574
575
576
577
578
579
	Matrix bmat = (*mat)[i][j]->getBaseMatrix();

	traits::col<Matrix>::type col(bmat);
	traits::const_value<Matrix>::type value(bmat);
	  
	typedef traits::range_generator<row, Matrix>::type cursor_type;
	typedef traits::range_generator<nz, cursor_type>::type icursor_type;
	
	for (cursor_type cursor = begin<row>(bmat), 
	       cend = end<row>(bmat); cursor != cend; ++cursor) {

580
	  int globalRowDof =
581
	    meshDistributor->mapDofToGlobal(feSpaces[i], *cursor);
582

583
	  // The corresponding global matrix row index of the current row DOF.
584
	  int petscRowIdx = dofToMatIndex.get(i, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
585
	  
586
	  if (meshDistributor->getIsRankDof(feSpaces[i], *cursor)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
587
	    	    
588
589
	    // === The current row DOF is a rank DOF, so create the       ===
	    // === corresponding nnz values directly on rank's nnz data.  ===
Thomas Witkowski's avatar
Thomas Witkowski committed
590
591
	    
	    // This is the local row index of the local PETSc matrix.
592
	    int localPetscRowIdx = petscRowIdx - rankStartIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
593
594
595
596
	    
	    TEST_EXIT_DBG(localPetscRowIdx >= 0 && localPetscRowIdx < nRankRows)
	      ("Should not happen! \n Debug info: localRowIdx = %d   globalRowIndx = %d   petscRowIdx = %d   localPetscRowIdx = %d   rStart = %d   nCompontens = %d   nRankRows = %d\n",
	       *cursor, 
597
	       meshDistributor->mapDofToGlobal(feSpaces[i], *cursor), 
Thomas Witkowski's avatar
Thomas Witkowski committed
598
599
	       petscRowIdx, 
	       localPetscRowIdx, 
600
	       rankStartIndex,
Thomas Witkowski's avatar
Thomas Witkowski committed
601
602
603
604
605
606
607
	       nComponents, 
	       nRankRows);
	    
	    
	    // Traverse all non zero entries in this row.
	    for (icursor_type icursor = begin<nz>(cursor), 
		   icend = end<nz>(cursor); icursor != icend; ++icursor) {
608
	      int globalColDof =
609
		meshDistributor->mapDofToGlobal(feSpaces[j], col(*icursor));
610
	      int petscColIdx = dofToMatIndex.get(j, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
611
612
613
	      
	      if (value(*icursor) != 0.0 || petscRowIdx == petscColIdx) {
		// The row DOF is a rank DOF, if also the column is a rank DOF, 
614
615
616
617
		// increment the d_nnz values for this row, otherwise the 
		// o_nnz value.
		if (petscColIdx >= rankStartIndex && 
		    petscColIdx < rankStartIndex + nRankRows)
Thomas Witkowski's avatar
Thomas Witkowski committed
618
619
620
621
622
623
		  d_nnz[localPetscRowIdx]++;
		else
		  o_nnz[localPetscRowIdx]++;
	      }    
	    }
	  } else {
624
625
626
627
628
	    // === The current row DOF is not a rank DOF, i.e., its values   ===
	    // === are also created on this rank, but afterthere they will   ===
	    // === be send to another rank. So we need to send also the      ===
	    // === corresponding nnz structure of this row to the corres-    ===
	    // === ponding rank.                                             ===
Thomas Witkowski's avatar
Thomas Witkowski committed
629
630
	    
	    // Send all non zero entries to the member of the row DOF.
631
	    int sendToRank = sendDofToRank[make_pair(*cursor, i)];
Thomas Witkowski's avatar
Thomas Witkowski committed
632
633
634
635
	    
	    for (icursor_type icursor = begin<nz>(cursor), 
		   icend = end<nz>(cursor); icursor != icend; ++icursor) {
	      if (value(*icursor) != 0.0) {
636
		int globalColDof =
637
		  meshDistributor->mapDofToGlobal(feSpaces[j], col(*icursor));
638
		int petscColIdx = dofToMatIndex.get(j, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
		
		sendMatrixEntry[sendToRank].
		  push_back(make_pair(petscRowIdx, petscColIdx));
	      }
	    }
	    
	  } // if (isRankDof[*cursor]) ... else ...
	} // for each row in mat[i][j]
      } 
    }

    // === Send and recv the nnz row structure to/from other ranks. ===

    StdMpi<MatrixNnzEntry> stdMpi(meshDistributor->getMpiComm(), true);
    stdMpi.send(sendMatrixEntry);
    for (std::set<int>::iterator it = recvFromRank.begin(); 
	 it != recvFromRank.end(); ++it)
      stdMpi.recv(*it);
    stdMpi.startCommunication();


660
661
    // === Evaluate the nnz structure this rank got from other ranks and add ===
    // === it to the PETSc nnz data structure.                               ===
Thomas Witkowski's avatar
Thomas Witkowski committed
662
663
664
665
666
667
668
669

    for (map<int, MatrixNnzEntry>::iterator it = stdMpi.getRecvData().begin();
	 it != stdMpi.getRecvData().end(); ++it) {
      if (it->second.size() > 0) {
	for (unsigned int i = 0; i < it->second.size(); i++) {
	  int r = it->second[i].first;
	  int c = it->second[i].second;

670
	  int localRowIdx = r - rankStartIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
671
672
673
674
675

	  TEST_EXIT_DBG(localRowIdx >= 0 && localRowIdx < nRankRows)
	    ("Got row index %d/%d (nRankRows = %d) from rank %d. Should not happen!\n",
	     r, localRowIdx, nRankRows, it->first);
	  
676
	  if (c < rankStartIndex || c >= rankStartIndex + nRankRows)
Thomas Witkowski's avatar
Thomas Witkowski committed
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
	    o_nnz[localRowIdx]++;
	  else
	    d_nnz[localRowIdx]++;
	}
      }
    }

    // The above algorithm for calculating the number of nnz per row over-
    // approximates the value, i.e., the number is always equal or larger to 
    // the real number of nnz values in the global parallel matrix. For small
    // matrices, the problem may arise, that the result is larger than the
    // number of elements in a row. This is fixed in the following.

    if (nRankRows < 100) 
      for (int i = 0; i < nRankRows; i++)
	d_nnz[i] = std::min(d_nnz[i], nRankRows);
  }

695
696
697
698
699
700
701
702

  void PetscSolverGlobalMatrix::createGlobalDofMapping(vector<const FiniteElemSpace*> &feSpaces)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createGlobalDofMapping()");

    int offset = meshDistributor->getStartDofs(feSpaces);
    Mesh *mesh = meshDistributor->getMesh();

703
    dofToMatIndex.clear();
704
  
705
706
707
708
709
710
711
712
    for (unsigned int i = 0; i < feSpaces.size(); i++) {

      // === Create indices for all DOFs in rank' domain. ===
      std::set<const DegreeOfFreedom*> rankDofSet;
      mesh->getAllDofs(feSpaces[i], rankDofSet);
      for (std::set<const DegreeOfFreedom*>::iterator it = rankDofSet.begin();
	   it != rankDofSet.end(); ++it)
	if (meshDistributor->getIsRankDof(feSpaces[i], **it)) {
713
	  int globalIndex = 
714
	    meshDistributor->mapDofToGlobal(feSpaces[i], **it);
715
716
717
718

 	  int globalMatIndex =  
 	    globalIndex - meshDistributor->getStartDofs(feSpaces[i]) + offset;

719
	  dofToMatIndex.add(i, globalIndex, globalMatIndex);
720
721
722
	}
	  

723
      // === Communicate interior boundary DOFs between domains. ===
724
725
726
727
728
729
730
731

      StdMpi<vector<int> > stdMpi(meshDistributor->getMpiComm());
    
      for (DofComm::Iterator it(meshDistributor->getSendDofs(), feSpaces[i]);
	   !it.end(); it.nextRank()) {
	vector<DegreeOfFreedom> sendGlobalDofs;

	for (; !it.endDofIter(); it.nextDof()) {
732
	  int globalIndex = 
733
	    meshDistributor->mapDofToGlobal(feSpaces[i], it.getDofIndex());
734
	  int globalMatIndex = dofToMatIndex.get(i, globalIndex);
735
	  sendGlobalDofs.push_back(globalMatIndex);
736
737
738
739
740
741
742
743
744
745
746
747
748
749
	}

	stdMpi.send(it.getRank(), sendGlobalDofs);
      }

      for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpaces[i]);
	   !it.end(); it.nextRank())
	stdMpi.recv(it.getRank());

      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getRecvDofs(), feSpaces[i]);
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof()) {
750
751
752
753
	  int globalIndex = 
	    meshDistributor->mapDofToGlobal(feSpaces[i], it.getDofIndex());
	  int globalMatIndex = 
	    stdMpi.getRecvData(it.getRank())[it.getDofCounter()];
754

755
	  dofToMatIndex.add(i, globalIndex, globalMatIndex);
756
757
758
759
760
761
762
763
764
765
766
767
	}


      // === Communicate DOFs on periodic boundaries. ===

      stdMpi.clear();

      for (DofComm::Iterator it(meshDistributor->getPeriodicDofs(), feSpaces[i]);
	   !it.end(); it.nextRank()) {
	vector<DegreeOfFreedom> sendGlobalDofs;
	
	for (; !it.endDofIter(); it.nextDof()) {
768
769
	  int ind0 = 
	    meshDistributor->mapDofToGlobal(feSpaces[i], it.getDofIndex());
770
	  int ind1 = dofToMatIndex.get(i, ind0);
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787

	  sendGlobalDofs.push_back(ind0);
	  sendGlobalDofs.push_back(ind1);
	}

	stdMpi.send(it.getRank(), sendGlobalDofs);
	stdMpi.recv(it.getRank());
      }

      stdMpi.startCommunication();

      for (DofComm::Iterator it(meshDistributor->getPeriodicDofs(), feSpaces[i]);
	   !it.end(); it.nextRank())
	for (; !it.endDofIter(); it.nextDof()) {
	  int ind0 = stdMpi.getRecvData(it.getRank())[it.getDofCounter() * 2];
	  int ind1 = stdMpi.getRecvData(it.getRank())[it.getDofCounter() * 2 + 1];

788
	  dofToMatIndex.add(i, ind0, ind1);
789
790
791
792
793
794
795
	}
      
      // === Update offset. ===

      offset += meshDistributor->getNumberRankDofs(feSpaces[i]);
    }
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
796
}