ParallelCoarseSpaceMatVec.cc 8.68 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
//
// 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.


#include "parallel/ParallelCoarseSpaceMatVec.h"
14
#include "parallel/ParallelDofMapping.h"
15
#include "parallel/MatrixNnzStructure.h"
16
17
18
19
20

namespace AMDiS {

  using namespace std;

21
  ParallelCoarseSpaceMatVec::ParallelCoarseSpaceMatVec()
22
    : interiorMap(NULL),
23
      lastMeshNnz(-1),
24
25
26
27
28
      alwaysCreateNnzStructure(false),
      meshDistributor(NULL),
      subdomainLevel(0),
      rStartInterior(0),
      nGlobalOverallInterior(0)
29
  {
30
31
32
33
34
    Parameters::get("parallel->always create nnz structure", 
		    alwaysCreateNnzStructure);
  }


35
36
  void ParallelCoarseSpaceMatVec::setCoarseSpaceDofMapping(ParallelDofMapping *coarseDofs, 
							   int component)
37
  {
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
    FUNCNAME("ParallelCoarseSpaceMatVec::setCoarseSpaceDofMapping()");

    TEST_EXIT_DBG(coarseDofs)("Should not happen!\n");

    if (component == -1) {
      // === Set coarse space for all components. ===

      coarseSpaceMap.clear();

      int nComponents = coarseDofs->getNumberOfComponents();
      for (int i = 0; i < nComponents; i++)
	coarseSpaceMap[i] = coarseDofs;
    } else {
      // === Set coarse space for just one component. ===

      coarseSpaceMap[component] = coarseDofs;
    }
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
55
56
57
58

    if (find(uniqueCoarseMap.begin(), uniqueCoarseMap.end(), coarseDofs) == 
	uniqueCoarseMap.end())
      uniqueCoarseMap.push_back(coarseDofs);
59
60
61
62
63
64
65
  }


  void ParallelCoarseSpaceMatVec::prepare()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec:prepare()");

66
67
    TEST_EXIT(uniqueCoarseMap.size() <= 2)
      ("Not yet implemented for more than two coarse spaces!\n");
68
69
70

    // === Create pointers to PETSc matrix and vector objects. ===

71
72
73
74
75
    int nCoarseMap = uniqueCoarseMap.size();
    mat.resize(nCoarseMap + 1);
    for (int i = 0; i < nCoarseMap + 1; i++)
      mat[i].resize(nCoarseMap + 1);

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
76
77
    vecSol.resize(nCoarseMap + 1);
    vecRhs.resize(nCoarseMap + 1);
78
   
79
80
    // === Create map from component number to its coarse space map. ===

81
    componentIthCoarseMap.resize(coarseSpaceMap.size());
Thomas Witkowski's avatar
bla    
Thomas Witkowski committed
82
    for (unsigned int i = 0; i < coarseSpaceMap.size(); i++) {
83
84
85
86
87
88
89
90
91
92
93
      bool found = false;
      for (int j = 0; j < nCoarseMap; j++) {
	if (coarseSpaceMap[i] == uniqueCoarseMap[j]) {
	  componentIthCoarseMap[i] = j;
	  found = true;
	  break;
	}
      }
      
      TEST_EXIT_DBG(found)("Should not happen!\n");
    }
94
95
96
  }


97
  void ParallelCoarseSpaceMatVec::createMatVec(Matrix<DOFMatrix*>& seqMat)
98
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
99
    FUNCNAME("ParallelCoarseSpaceMatVec::createMatVec()");
100

101
102
103
104
105
106
107
108
    // === Prepare coarse space information and generate the correct number ===
    // === of empty PETSc matrix and vector objects.                        ===

    prepare();

    
    // === Update subdomain data (mostly required for multi-level methods) ===

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
109
    updateSubdomainData();
110

111

112
113
114
    // === If required, recompute non zero structure of the matrix. ===

    bool localMatrix = (coarseSpaceMap.size() && subdomainLevel == 0);
Thomas Witkowski's avatar
Thomas Witkowski committed
115

116
    if (checkMeshChange()) {
Thomas Witkowski's avatar
bo eh    
Thomas Witkowski committed
117
       int nMat = uniqueCoarseMap.size() + 1;
118
119
      nnz.resize(nMat);
      for (int i = 0; i < nMat; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
120
	nnz[i].resize(nMat);
121
122
123
124
125
126
127
128
129
130
131
132
133
134
	for (int j = 0; j < nMat; j++)
	  nnz[i][j].clear();
      }

      nnz[0][0].create(seqMat, mpiCommGlobal, *interiorMap,
		       (coarseSpaceMap.size() == 0 ? &(meshDistributor->getPeriodicMap()) : NULL),
		       meshDistributor->getElementObjectDb(),
		       localMatrix);

      for (int i = 0; i < nMat; i++) {
	for (int j = 0; j < nMat; j++) {
	  if (i == 0 && j == 0)
	    continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
135
136
137
138
139
140
141
	  ParallelDofMapping &rowMap = 
	    (i == 0 ? *interiorMap : *(uniqueCoarseMap[i - 1]));
	  ParallelDofMapping &colMap =
	    (j == 0 ? *interiorMap : *(uniqueCoarseMap[j - 1]));

	  nnz[i][j].create(seqMat, mpiCommGlobal, rowMap, colMap, NULL,
			   meshDistributor->getElementObjectDb());
142
143
144
145
146
	}
      }
    }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
147
148
    // === Create PETSc matrices and PETSc vectors with the computed  ===
    // === nnz data structure.                                        ===
149
    
Thomas Witkowski's avatar
Thomas Witkowski committed
150
151
    int nRankInteriorRows = interiorMap->getRankDofs();
    int nOverallInteriorRows = interiorMap->getOverallDofs();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
152
    
153
    if (localMatrix) {
Thomas Witkowski's avatar
Thomas Witkowski committed
154
      MatCreateSeqAIJ(mpiCommLocal, nRankInteriorRows, nRankInteriorRows,
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
155
		      0, nnz[0][0].dnnz,
156
157
158
		      &mat[0][0]);
      MatSetOption(mat[0][0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
159
160
      MatCreateAIJ(mpiCommGlobal, nRankInteriorRows, nRankInteriorRows, 
		   nOverallInteriorRows, nOverallInteriorRows,
Thomas Witkowski's avatar
Thomas Witkowski committed
161
		   0, nnz[0][0].dnnz, 0, nnz[0][0].onnz,		   
162
163
164
165
		   &mat[0][0]);
      MatSetOption(mat[0][0], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
166
167
    VecCreateMPI(mpiCommGlobal, nRankInteriorRows, nOverallInteriorRows, &vecSol[0]);
    VecCreateMPI(mpiCommGlobal, nRankInteriorRows, nOverallInteriorRows, &vecRhs[0]);
168

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
169
170
171
172
    int nCoarseMap = uniqueCoarseMap.size();
    for (int i = 0; i < nCoarseMap; i++) {
      ParallelDofMapping* cMap = uniqueCoarseMap[i];
      
Thomas Witkowski's avatar
Thomas Witkowski committed
173
174
      int nRankCoarseRows = cMap->getRankDofs();
      int nOverallCoarseRows = cMap->getOverallDofs();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
175
176
      
      MatCreateAIJ(mpiCommGlobal,
Thomas Witkowski's avatar
Thomas Witkowski committed
177
178
		   nRankCoarseRows, nRankCoarseRows,
		   nOverallCoarseRows, nOverallCoarseRows,
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
179
180
		   0, nnz[i + 1][i + 1].dnnz, 0, nnz[i + 1][i + 1].onnz,
		   &mat[i + 1][i + 1]);
181
182
      cMap->createVec(vecSol[i + 1]);
      cMap->createVec(vecRhs[i + 1]);
Thomas Witkowski's avatar
Thomas Witkowski committed
183
184
185
    }

    for (int i = 0; i < nCoarseMap + 1; i++) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
186
      for (int j = 0; j < nCoarseMap + 1; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
189
190
191
192
193
194
	if (i == j)
	  continue;

	int nRowsRankMat = (i == 0 ? nRankInteriorRows : uniqueCoarseMap[i - 1]->getRankDofs());
	int nRowsOverallMat = (i == 0 ? nOverallInteriorRows : uniqueCoarseMap[i - 1]->getOverallDofs());

	int nColsRankMat = (j == 0 ? nRankInteriorRows : uniqueCoarseMap[j - 1]->getRankDofs());
	int nColsOverallMat = (j == 0 ? nOverallInteriorRows : uniqueCoarseMap[j - 1]->getOverallDofs());
195

196
	MatCreateAIJ(mpiCommGlobal,
Thomas Witkowski's avatar
Thomas Witkowski committed
197
198
199
		     nRowsRankMat, nColsRankMat,
		     nRowsOverallMat, nColsOverallMat,
		     0, nnz[i][j].dnnz, 0, nnz[i][j].onnz,
200
		     &mat[i][j]);	  
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
201
	MatCreateAIJ(mpiCommGlobal,
Thomas Witkowski's avatar
Thomas Witkowski committed
202
203
204
205
		     nColsRankMat, nRowsRankMat,
		     nColsOverallMat, nRowsOverallMat,
		     0, nnz[j][i].dnnz, 0, nnz[j][i].onnz,
		     &mat[j][i]);
206
      }
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
207
    }    
208
209
210
  }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
211
  void ParallelCoarseSpaceMatVec::matDestroy()
212
  {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
213
    FUNCNAME("ParallelCoarseSpaceMatVec::matDestroy()");
214
215
216
217
218
219
220
221

    int nMatrix = mat.size();
    for (int i = 0; i < nMatrix; i++)
      for (int j = 0; j < nMatrix; j++)
	MatDestroy(&mat[i][j]);
  }


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
222
  void ParallelCoarseSpaceMatVec::vecDestroy()
223
  {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
224
225
226
227
228
229
230
231
232
233
234
235
236
    FUNCNAME("ParallelCoarseSpaceMatVec::vecDestroy()");

    int nVec = vecSol.size();
    for (int i = 0; i < nVec; i++) {
      VecDestroy(&vecSol[i]);
      VecDestroy(&vecRhs[i]);
    }
  }


  void ParallelCoarseSpaceMatVec::matAssembly()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::matAssembly()");
237
238
239
240
241
242
243
244
245
246

    int nMatrix = mat.size();
    for (int i = 0; i < nMatrix; i++) {
      for (int j = 0; j < nMatrix; j++) {
	MatAssemblyBegin(mat[i][j], MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(mat[i][j], MAT_FINAL_ASSEMBLY);  
      }
    }
  }

247

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
  void ParallelCoarseSpaceMatVec::vecRhsAssembly()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::vecRhsAssembly()");

    int nVec = vecRhs.size();
    for (int i = 0; i < nVec; i++) {
      VecAssemblyBegin(vecRhs[i]);
      VecAssemblyEnd(vecRhs[i]);
    }
  }


  void ParallelCoarseSpaceMatVec::vecSolAssembly()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::vecSolAssembly()");

    int nVec = vecRhs.size();
    for (int i = 0; i < nVec; i++) {
      VecAssemblyBegin(vecSol[i]);
      VecAssemblyEnd(vecSol[i]);
    }
  }


272
  bool ParallelCoarseSpaceMatVec::checkMeshChange()
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::checkMeshChange()");

    int recvAllValues = 0;
    int sendValue = 
      static_cast<int>(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz);
    mpiCommGlobal.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);

    if (recvAllValues != 0 || alwaysCreateNnzStructure) {
      lastMeshNnz = meshDistributor->getLastMeshChangeIndex();
      return true;
    }

    return false;
  }

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312

  void ParallelCoarseSpaceMatVec::updateSubdomainData()
  {
    FUNCNAME("ParallelCoarseSpaceMatVec::updateSubdomainData()");

    if (mpiCommLocal.Get_size() == 1) {
      rStartInterior = 0;
      nGlobalOverallInterior = interiorMap->getOverallDofs();
    } else {
      int groupRowsInterior = 0;
      if (mpiCommLocal.Get_rank() == 0)
	groupRowsInterior = interiorMap->getOverallDofs();

      mpi::getDofNumbering(mpiCommGlobal, groupRowsInterior,
			   rStartInterior, nGlobalOverallInterior);

      int tmp = 0;
      if (mpiCommLocal.Get_rank() == 0)
	tmp = rStartInterior;

      mpiCommLocal.Allreduce(&tmp, &rStartInterior, 1, MPI_INT, MPI_SUM);
    }
  }

313
}