MatrixNnzStructure.cc 11.9 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20
21


Thomas Witkowski's avatar
Merge.  
Thomas Witkowski committed
22
#include <algorithm>
23
24
25
26
27
#include "Global.h"
#include "DOFMatrix.h"
#include "parallel/MatrixNnzStructure.h"
#include "parallel/ParallelDofMapping.h"

28
29
using namespace std;

30
namespace AMDiS { namespace Parallel {
31

Thomas Witkowski's avatar
Thomas Witkowski committed
32
33
34
35
36
37
  MatrixNnzStructure::~MatrixNnzStructure()
  {
    clear();
  }


38
39
40
41
  void MatrixNnzStructure::clear()
  {
    if (dnnz) {
      delete [] dnnz;
42
      dnnz = NULL;
43
    }
44

45
46
    if (onnz) {
      delete [] onnz;
47
      onnz = NULL;
48
49
50
    }
  }

51
  void MatrixNnzStructure::create(Matrix<DOFMatrix*> &mat,
52
53
				  ParallelDofMapping &rowDofMap,
				  ParallelDofMapping &colDofMap,
Thomas Witkowski's avatar
Thomas Witkowski committed
54
				  PeriodicMap *perMap,
55
				  const ElementObjectDatabase& elObjDb,
56
57
58
59
				  bool localMatrix)
  {
    FUNCNAME("MatrixNnzStructure::create()");

60
    int nRankRows = rowDofMap.getRankDofs();	// Number of DOFs owned by rank.
61

62
    int rankStartRowIndex = rowDofMap.getStartDofs();		// Smallest global index of a DOF owned by the rank.
63
64

    int nRankCols = colDofMap.getRankDofs();
65

66
    int nOverallCols = colDofMap.getOverallDofs();	// Number of global DOFs (this value is thus the same on all ranks).
67

68
69
70
71
72
73
74
75
76
77
78
    int rankStartColIndex = colDofMap.getStartDofs();

    create(nRankRows, (!localMatrix ? nRankRows : -1));

    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;

    typedef vector<pair<int, int> > MatrixNnzEntry;
79
//     typedef map<int, DofContainer> RankToDofContainer;
80
81
82

    // 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
83
    // on this rank, but because the row DOFs are not DOFs of this rank they
84
85
86
87
88
89
90
    // will be send to the owner of the row DOFs.
    map<int, MatrixNnzEntry> sendMatrixEntry;

    // Maps to each DOF that must be send to another rank the rank number of the
    // receiving rank.
    map<pair<DegreeOfFreedom, int>, int> sendDofToRank;

91
    int nComponents = mat.getNumRows();
92

93
    // First, create for all ranks, to which we send data to, MatrixNnzEntry
94
95
    // object with 0 entries.
    for (int i = 0; i < nComponents; i++) {
96
      const FiniteElemSpace* feSpace = NULL;
97
      for (int j = 0; j < nComponents; j++)
98
99
	if (mat[i][j])
	  feSpace = mat[i][j]->getRowFeSpace();
100

101
102
      TEST_EXIT_DBG(feSpace)("No FE space found!\n");

103
      for (DofComm::Iterator it(rowDofMap.getDofComm(feSpace).getRecvDofs(), feSpace);
104
105
	   !it.end(); it.nextRank()) {
	sendMatrixEntry[it.getRank()].resize(0);
106

107
108
109
110
111
112
113
114
	for (; !it.endDofIter(); it.nextDof())
	  sendDofToRank[make_pair(it.getDofIndex(), i)] = it.getRank();
      }
    }

    // Create list of ranks from which we receive data from.
    std::set<int> recvFromRank;
    for (int i = 0; i < nComponents; i++)  {
115
      const FiniteElemSpace* feSpace = NULL;
116
      for (int j = 0; j < nComponents; j++)
117
118
	if (mat[i][j])
	  feSpace = mat[i][j]->getRowFeSpace();
119

120
      for (DofComm::Iterator it(rowDofMap.getDofComm(feSpace).getSendDofs(), feSpace);
121
122
123
124
125
126
127
128
129
	   !it.end(); it.nextRank())
	recvFromRank.insert(it.getRank());
    }


    // === Traverse matrices to create nnz data. ===

    for (int rowComp = 0; rowComp < nComponents; rowComp++) {
      for (int colComp = 0; colComp < nComponents; colComp++) {
130
	DOFMatrix *dofMat = mat[rowComp][colComp];
131
132
133
134
135
136

 	if (!dofMat)
	  continue;

	const FiniteElemSpace *rowFeSpace = dofMat->getRowFeSpace();
	const FiniteElemSpace *colFeSpace = dofMat->getColFeSpace();
137

138
139
	if (rowDofMap.isDefinedFor(rowComp) == false ||
	    colDofMap.isDefinedFor(colComp) == false)
Thomas Witkowski's avatar
Thomas Witkowski committed
140
	  continue;
141
142
143

	// === Prepare MTL4 iterators for the current DOF matrix. ===

144
	Matrix baseMat = dofMat->getBaseMatrix(); // TBD: warum hier keine Referenz?
145

146
147
	traits::col<Matrix>::type getCol(baseMat);
	traits::row<Matrix>::type getRow(baseMat);
148
149

	traits::const_value<Matrix>::type value(baseMat);
150
151
152
153
154
155
	cursor_type cursor = begin<row>(baseMat);
	cursor_type cend = end<row>(baseMat);


	// === Iterate on all DOFs of the row mapping. ===

156
	DofMap::iterator rowIt = rowDofMap[rowComp].begin(); // row dofmap DOFIt
157
	DofMap::iterator rowEndIt = rowDofMap[rowComp].end();
158
159
160
161
	for (; rowIt != rowEndIt; ++rowIt) {

	  // Go to the corresponding matrix row (note, both the mapping and the
	  // matrix rows are stored in increasing order).
162
	  while (cursor.value() != rowIt->first)
163
	    ++cursor;
164
	  size_t _row = cursor.value(); // DegreeOfFreedom
165
166
167
168

	  // The corresponding global matrix row index of the current row DOF.
	  int petscRowIdx = 0;
	  if (localMatrix) {
169
	    petscRowIdx = rowDofMap.getLocalMatIndex(rowComp, _row);
170
171
172
173
	  } else {
	    if (rowDofMap.isMatIndexFromGlobal())
	      petscRowIdx = rowDofMap.getMatIndex(rowComp, rowIt->second.global);
	    else
174
	      petscRowIdx = rowDofMap.getMatIndex(rowComp, _row);
175
176
	  }

177

178
	  // Set of periodic row associations (is empty, if row DOF is not
179
180
181
182
	  // periodic.
	  std::set<int> perRowAsc;
	  if (perMap)
	    perMap->fillAssociations(rowFeSpace, rowIt->second.global, elObjDb, perRowAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
183

184

185
	  if (localMatrix || rowDofMap[rowComp].isRankDof(_row)) {
186
187
	    // === The current row DOF is a rank DOF, so create the       ===
	    // === corresponding nnz values directly on rank's nnz data.  ===
188

189
190
191
192
193
	    // This is the local row index of the local PETSc matrix.
	    int localPetscRowIdx = petscRowIdx;

	    if (localMatrix == false)
	      localPetscRowIdx -= rankStartRowIndex;
194

195
	    TEST_EXIT(localPetscRowIdx >= 0 && localPetscRowIdx < nRankRows)
196
	      ("Should not happen! \n Debug info: DOF = %d   globalRowIndx = %d   petscRowIdx = %d   localPetscRowIdx = %d   rStart = %d   compontens = %d from %d   nRankRows = %d\n",
197
198
	       _row,
	       rowDofMap[rowComp][_row].global,
199
200
	       petscRowIdx,
	       localPetscRowIdx,
201
202
	       rankStartRowIndex,
	       rowComp,
203
	       nComponents,
204
	       nRankRows);
205

206
207

	    if (localMatrix) {
208
	      for (icursor_type icursor = begin<nz>(cursor),
209
		     icend = end<nz>(cursor); icursor != icend; ++icursor)
210
		if (colDofMap[colComp].isSet(getCol(*icursor)))
211
		  dnnz[localPetscRowIdx]++;
212
213
214
215
	    } else {
	      MultiIndex colDofIndex;

	      // Traverse all non zero entries in this row.
216
	      for (icursor_type icursor = begin<nz>(cursor),
217
		     icend = end<nz>(cursor); icursor != icend; ++icursor) {
218
		if (colDofMap[colComp].find(getCol(*icursor), colDofIndex) == false)
219
		  continue;
220
221

		// Set of periodic row associations (is empty, if row DOF is not
222
		// periodic.
Thomas Witkowski's avatar
Blbu  
Thomas Witkowski committed
223
		std::set<int> perColAsc = perRowAsc;
224
225
		if (perMap)
		  perMap->fillAssociations(colFeSpace, colDofIndex.global, elObjDb, perColAsc);
226

227
		if (perColAsc.empty()) {
228
		  if (colDofMap[colComp].isRankDof(getCol(*icursor)))
229
230
231
232
233
234
		    dnnz[localPetscRowIdx]++;
		  else
		    onnz[localPetscRowIdx]++;
		} else {
		  vector<int> newCols;
		  perMap->mapDof(colFeSpace, colDofIndex.global, perColAsc, newCols);
235

236
		  for (size_t aa = 0; aa < newCols.size(); aa++) {
237
		    int petscColIdx = colDofMap.getMatIndex(colComp, newCols[aa]);
238
239
240

		    // The row DOF is a rank DOF, if also the column is a rank DOF,
		    // increment the diagNnz values for this row,
241
		    // otherwise the offdiagNnz value.
242
		    if (petscColIdx >= rankStartColIndex &&
243
244
245
246
247
248
			petscColIdx < rankStartColIndex + nRankCols)
		      dnnz[localPetscRowIdx]++;
		    else
		      onnz[localPetscRowIdx]++;
		  }
		}
Thomas Witkowski's avatar
Thomas Witkowski committed
249
	      }
250

Thomas Witkowski's avatar
Blbu  
Thomas Witkowski committed
251
252
253
 	      if (!perRowAsc.empty()) {
		vector<int> newRows;
		perMap->mapDof(rowFeSpace, rowIt->second.global, perRowAsc, newRows);
254
255

		dnnz[localPetscRowIdx] +=
Thomas Witkowski's avatar
Blbu  
Thomas Witkowski committed
256
257
		  (newRows.size() - 1) * (onnz[localPetscRowIdx] + dnnz[localPetscRowIdx]);
 	      }
258
259
260
261
262
263
264
	    }
	  } else {
	    // === 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.                                             ===
265

266
	    // Send all non zero entries to the member of the row DOF.
267
	    int sendToRank = sendDofToRank[make_pair(cursor.value(), rowComp)];
268
269
	    MultiIndex colDofIndex;

270
	    for (icursor_type icursor = begin<nz>(cursor),
271
		   icend = end<nz>(cursor); icursor != icend; ++icursor) {
272
	      if (colDofMap[colComp].find(getCol(*icursor), colDofIndex) == false)
273
274
		continue;

275
	      int petscColIdx = (colDofMap.isMatIndexFromGlobal() ?
276
				 colDofMap.getMatIndex(colComp, colDofIndex.global) :
277
				 colDofMap.getMatIndex(colComp, getCol(*icursor)));
278

279
280
281
	      sendMatrixEntry[sendToRank].
		push_back(make_pair(petscRowIdx, petscColIdx));
	    }
282

283
	  } // if (isRankDof[cursor.value()]) ... else ...
284
	} // for each row in mat[i][j]
285
      }
286
287
288
289
290
    }


    if (localMatrix == false) {
      // === Send and recv the nnz row structure to/from other ranks. ===
291

292
      StdMpi<MatrixNnzEntry> stdMpi(rowDofMap.getMpiComm());
293
      stdMpi.send(sendMatrixEntry);
294

295
      for (std::set<int>::iterator it = recvFromRank.begin();
296
297
	   it != recvFromRank.end(); ++it)
	stdMpi.recv(*it);
298
299
      stdMpi.startCommunication();

300
301
      // === Evaluate the nnz structure this rank got from other ranks and add ===
      // === it to the PETSc nnz data structure.                               ===
302

303
304
305
306
307
308
      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;
309

310
	    int localRowIdx = r - rankStartRowIndex;
311

312
313
314
	    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);
315

316
317
318
319
320
321
322
323
324
325
	    if (c < rankStartColIndex || c >= rankStartColIndex + nRankCols)
	      onnz[localRowIdx]++;
	    else
	      dnnz[localRowIdx]++;
	  }
	}
      }
    }

    // The above algorithm for calculating the number of nnz per row over-
326
    // approximates the value, i.e., the number is always equal or larger to
327
328
329
330
    // 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.

Thomas Witkowski's avatar
Thomas Witkowski committed
331
332
    for (int i = 0; i < nRankRows; i++) {
      dnnz[i] = std::min(dnnz[i], nRankCols);
Thomas Witkowski's avatar
Bugfix  
Thomas Witkowski committed
333
334
      if (onnz)
	onnz[i] = std::min(onnz[i], nOverallCols - nRankCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
335
    }
336

337
#ifndef NDEBUG
338
    int nMax = 0;
339
340
341
342
343
344
345
346
347
348
349
350
351
352
    int nSum = 0;
    for (int i = 0; i < nRankRows; i++) {
      nMax = std::max(nMax, dnnz[i]);
      nSum += dnnz[i];
    }

    MSG("NNZ in diag block: max = %d, avrg = %.0f\n",
	nMax, (nSum > 0 ? (static_cast<double>(nSum) / nRankRows) : 0));

    if (!localMatrix) {
      nMax = 0;
      nSum = 0;

      for (int i = 0; i < nRankRows; i++) {
Thomas Witkowski's avatar
Merge.  
Thomas Witkowski committed
353
354
	nMax = std::max(nMax, onnz[i]);
	nSum += onnz[i];
355
356
357
      }

      MSG("NNZ in offdiag block: max = %d, avrg = %.0f\n",
358
	  nMax, (nSum > 0 ? (static_cast<double>(nSum) / nRankRows) : 0));
359
360
361
362
363
364
365
366
367
    }
#endif
  }


  void MatrixNnzStructure::create(int nRows0, int nRows1)
  {
    if (nRows0 == 0)
      return;
368

369
    TEST_EXIT_DBG(nRows0 > 0)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
370
371

    clear();
372

373
374
375
    dnnz = new int[nRows0];
    for (int i = 0; i < nRows0; i++)
      dnnz[i] = 0;
376

377
378
379
    if (nRows1 > 0) {
      onnz = new int[nRows1];
      for (int i = 0; i < nRows1; i++)
380
	onnz[i] = 0;
381
382
383
    }
  }

384
} }