ParallelCoarseSpaceMatVec.cc 8.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
//
// 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 "AMDiS.h"
#include "parallel/ParallelCoarseSpaceMatVec.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
117
118
    if (checkMeshChange()) {
      // Mesh has been changed, recompute interior DOF mapping.
      interiorMap->setComputeMatIndex(!localMatrix);
119
      interiorMap->update();
120

121
122
123
      int nMat = uniqueCoarseMap.size() + 1;
      nnz.resize(nMat);
      for (int i = 0; i < nMat; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
124
	nnz[i].resize(nMat);
125
126
127
128
129
130
131
132
133
134
135
136
137
138
	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
139
140
141
142
143
144
145
	  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());
146
147
148
149
150
	}
      }
    }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
170
171
    VecCreateMPI(mpiCommGlobal, nRankInteriorRows, nOverallInteriorRows, &vecSol[0]);
    VecCreateMPI(mpiCommGlobal, nRankInteriorRows, nOverallInteriorRows, &vecRhs[0]);
172

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

    for (int i = 0; i < nCoarseMap + 1; i++) {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
190
      for (int j = 0; j < nCoarseMap + 1; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
191
192
193
194
195
196
197
198
	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());
199

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


Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
215
  void ParallelCoarseSpaceMatVec::matDestroy()
216
  {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
217
    FUNCNAME("ParallelCoarseSpaceMatVec::matDestroy()");
218
219
220
221
222
223
224
225

    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
226
  void ParallelCoarseSpaceMatVec::vecDestroy()
227
  {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
228
229
230
231
232
233
234
235
236
237
238
239
240
    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()");
241
242
243
244
245
246
247
248
249
250

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

251

Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
  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]);
    }
  }


276
  bool ParallelCoarseSpaceMatVec::checkMeshChange()
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
  {
    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
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316

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

317
}