PetscHelper.cc 9.75 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
22


#include "parallel/PetscHelper.h"
23
#include "parallel/PetscSolver.h"
24
#include "Global.h"
25

26
27
28
29
30
namespace AMDiS
{
  namespace Parallel
  {
    namespace petsc_helper
31
    {
32
33

      using namespace std;
34

35
36
37
38
      void getMatLocalColumn(Mat mat, PetscMatCol &matCol)
      {
	PetscInt firstIndex, lastIndex;
	MatGetOwnershipRange(mat, &firstIndex, &lastIndex);
39

40
41
42
	PetscInt nCols;
	const PetscInt *cols;
	const PetscScalar *values;
43

44
45
	for (int row = firstIndex; row < lastIndex; row++) {
	  MatGetRow(mat, row, &nCols, &cols, &values);
46

47
48
49
50
51
	  for (int i = 0; i < nCols; i++) {
	    if (values[i] != 0.0) {
	      matCol[cols[i]].first.push_back(row - firstIndex);
	      matCol[cols[i]].second.push_back(values[i]);
	    }
52
53
	  }

54
55
	  MatRestoreRow(mat, row, &nCols, &cols, &values);
	}
56
57
58
      }


59
60
61
62
      void setMatLocalColumn(Mat mat, int column, Vec vec)
      {
	PetscInt firstIndex;
	MatGetOwnershipRange(mat, &firstIndex, PETSC_NULL);
63

64
65
	PetscInt vecSize;
	VecGetLocalSize(vec, &vecSize);
66

67
68
69
	PetscScalar *tmpValues;
	VecGetArray(vec, &tmpValues);
	for (int i  = 0; i < vecSize; i++)
70
	  MatSetValue(mat,
71
72
73
74
75
76
		      firstIndex + i,
		      column,
		      tmpValues[i],
		      ADD_VALUES);
	VecRestoreArray(vec, &tmpValues);
      }
77
78


79
80
81
      void getColumnVec(const SparseCol &matCol, Vec vec)
      {
	VecSet(vec, 0.0);
82
	VecSetValues(vec, matCol.first.size(),
83
84
85
86
		    &(matCol.first[0]), &(matCol.second[0]), INSERT_VALUES);
	VecAssemblyBegin(vec);
	VecAssemblyEnd(vec);
      }
87
88


89
90
91
      void blockMatMatSolve(MPI::Intracomm mpiComm, KSP ksp, Mat mat0, Mat &mat1)
      {
	FUNCNAME("blockMatMatSolve()");
92

93
94
95
96
97
	// === We have to  calculate mat1 = ksp mat0:                       ===
	// ===    - get all local column vectors from mat0                  ===
	// ===    - solve with ksp for this column vector as the rhs vector ===
	// ===    - set the result to the corresponding column vector of    ===
	// ===      matrix mat1                                             ===
98

99
100
101
102
103
104
105
106
107
	// Transform matrix mat0 into a sparse column format.
	PetscMatCol mat0_cols;
	getMatLocalColumn(mat0, mat0_cols);

	int nFirstCol, nLastCol;
	MatGetOwnershipRangeColumn(mat0, &nFirstCol, &nLastCol);

	int dnnz = 0;
	int onnz = 0;
108
	for (PetscMatCol::iterator it = mat0_cols.begin();
109
110
111
112
113
114
	    it != mat0_cols.end(); ++it) {
	  if (it->first >= nFirstCol && it->first < nLastCol)
	    dnnz++;
	  else
	    onnz++;
	}
115
116


117
118
119
	PetscInt localRows, localCols, globalRows, globalCols;
	MatGetLocalSize(mat0, &localRows, &localCols);
	MatGetSize(mat0, &globalRows, &globalCols);
120

121
122
123
124
	MatCreateAIJ(mpiComm,
		    localRows, localCols, globalRows, globalCols,
		    dnnz, PETSC_NULL, onnz, PETSC_NULL, &mat1);
	MatSetOption(mat1, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
125
126


127
128
129
130
	Vec tmpVec;
	VecCreateSeq(PETSC_COMM_SELF, localRows, &tmpVec);

	// Solve for all column vectors of mat A_BPi and create matrix matK
131
	for (PetscMatCol::iterator it = mat0_cols.begin();
132
133
134
135
136
	    it != mat0_cols.end(); ++it) {
	  getColumnVec(it->second, tmpVec);
	  KSPSolve(ksp, tmpVec, tmpVec);
	  setMatLocalColumn(mat1, it->first, tmpVec);
	}
137

138
	VecDestroy(&tmpVec);
139

140
141
142
	MatAssemblyBegin(mat1, MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(mat1, MAT_FINAL_ASSEMBLY);
      }
143
144


145
      void matNestConvert(Mat matNest, Mat &mat)
146
147
      {
	FUNCNAME("matNestConvert()");
148

149
150
	PetscInt nestRows, nestCols;
	MatNestGetSize(matNest, &nestRows, &nestCols);
151

152
153
	TEST_EXIT(nestRows == 2 && nestCols == 2)
	  ("This funciton is only implemented for 2x2 nested matrices!\n");
154

155
156
157
158
159
	Mat mat00, mat01, mat10, mat11;
	MatNestGetSubMat(matNest, 0, 0, &mat00);
	MatNestGetSubMat(matNest, 0, 1, &mat01);
	MatNestGetSubMat(matNest, 1, 0, &mat10);
	MatNestGetSubMat(matNest, 1, 1, &mat11);
160

161
162
163
	int nRankRows0, nOverallRows0;
	MatGetLocalSize(mat00, &nRankRows0, PETSC_NULL);
	MatGetSize(mat00, &nOverallRows0, PETSC_NULL);
164

165
166
167
	int nRankRows1, nOverallRows1;
	MatGetLocalSize(mat11, &nRankRows1, PETSC_NULL);
	MatGetSize(mat11, &nOverallRows1, PETSC_NULL);
168

169
170
	int firstRow0;
	MatGetOwnershipRange(mat00, &firstRow0, PETSC_NULL);
171

172
173
	int firstRow1;
	MatGetOwnershipRange(mat11, &firstRow1, PETSC_NULL);
174

175
176
177
	int nRankRows = nRankRows0 + nRankRows1;
	int nOverallRows = nOverallRows0 + nOverallRows1;
	int firstRow = firstRow0 + firstRow1;
178

179
	int mpiSize = MPI::COMM_WORLD.Get_size();
180
#ifndef NDEBUG
181
	int mpiRank = MPI::COMM_WORLD.Get_rank();
182
#endif
183
184
185
186
	vector<int> allFirstRow0(mpiSize + 1, 0);
	vector<int> allFirstRow1(mpiSize + 1, 0);
	MPI::COMM_WORLD.Allgather(&nRankRows0, 1, MPI_INT, &(allFirstRow0[1]), 1, MPI_INT);
	MPI::COMM_WORLD.Allgather(&nRankRows1, 1, MPI_INT, &(allFirstRow1[1]), 1, MPI_INT);
187

188
189
190
191
	for (int i = 1; i < mpiSize + 1; i++) {
	  allFirstRow0[i] += allFirstRow0[i - 1];
	  allFirstRow1[i] += allFirstRow1[i - 1];
	}
192

193
194
195
	TEST_EXIT_DBG(allFirstRow0[mpiRank] == firstRow0)("Should not happen!\n");
	TEST_EXIT_DBG(allFirstRow1[mpiRank] == firstRow1)("Should not happen!\n");

196
	MatCreateAIJ(PETSC_COMM_WORLD,
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
		    nRankRows, nRankRows,
		    nOverallRows, nOverallRows,
		    25, PETSC_NULL, 25, PETSC_NULL, &mat);
	MatSetOption(mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

	PetscInt nCols;
	const PetscInt *cols;
	const PetscScalar *vals;

	for (int i = 0; i < nRankRows0; i++) {
	  PetscInt newRowIndex = firstRow + i;

	  MatGetRow(mat00, firstRow0 + i, &nCols, &cols, &vals);
	  for (int j = 0; j < nCols; j++) {
	    int rank = -1;
	    for (int k = 0; k < mpiSize; k++) {
	      if (cols[j] >= allFirstRow0[k] && cols[j] < allFirstRow0[k + 1]) {
		rank = k;
		break;
	      }
217
	    }
218
	    TEST_EXIT_DBG(rank != -1)("Should not happen!\n");
219

220
221
	    PetscInt newColIndex = cols[j] + allFirstRow1[rank];
	    MatSetValue(mat, newRowIndex, newColIndex, vals[j], INSERT_VALUES);
222
	  }
223
224
225
226
227
228
229
230
231
232
233
234
	  MatRestoreRow(mat00, firstRow0 + i, &nCols, &cols, &vals);

	  MatGetRow(mat01, firstRow0 + i, &nCols, &cols, &vals);
	  for (int j = 0; j < nCols; j++) {
	    int rank = -1;
	    for (int k = 0; k < mpiSize; k++) {
	      if (cols[j] >= allFirstRow1[k] && cols[j] < allFirstRow1[k + 1]) {
		rank = k;
		break;
	      }
	    }
	    TEST_EXIT_DBG(rank != -1)("Should not happen!\n");
235

236
237
238
	    PetscInt newColIndex = cols[j] + allFirstRow0[rank + 1];
	    MatSetValue(mat, newRowIndex, newColIndex, vals[j], INSERT_VALUES);
	  }
239
	  MatRestoreRow(mat01, firstRow0 + i, &nCols, &cols, &vals);
240
	}
241

242
243
244
245
246
247
248
249
250
251
252
	for (int i = 0; i < nRankRows1; i++) {
	  PetscInt newRowIndex = firstRow + nRankRows0 + i;

	  MatGetRow(mat10, firstRow1 + i, &nCols, &cols, &vals);
	  for (int j = 0; j < nCols; j++) {
	    int rank = -1;
	    for (int k = 0; k < mpiSize; k++) {
	      if (cols[j] >= allFirstRow0[k] && cols[j] < allFirstRow0[k + 1]) {
		rank = k;
		break;
	      }
253
	    }
254
	    TEST_EXIT_DBG(rank != -1)("Should not happen!\n");
255

256
257
	    PetscInt newColIndex = cols[j] + allFirstRow1[rank];
	    MatSetValue(mat, newRowIndex, newColIndex, vals[j], INSERT_VALUES);
258
	  }
259
260
261
262
263
264
265
266
267
268
269
270
	  MatRestoreRow(mat10, firstRow1 + i, &nCols, &cols, &vals);

	  MatGetRow(mat11, firstRow1 + i, &nCols, &cols, &vals);
	  for (int j = 0; j < nCols; j++) {
	    int rank = -1;
	    for (int k = 0; k < mpiSize; k++) {
	      if (cols[j] >= allFirstRow1[k] && cols[j] < allFirstRow1[k + 1]) {
		rank = k;
		break;
	      }
	    }
	    TEST_EXIT_DBG(rank != -1)("Should not happen!\n");
271

272
273
274
275
	    PetscInt newColIndex = cols[j] + allFirstRow0[rank + 1];
	    MatSetValue(mat, newRowIndex, newColIndex, vals[j], INSERT_VALUES);
	  }
	  MatRestoreRow(mat11, firstRow1 + i, &nCols, &cols, &vals);
276
	}
277
278
279

	MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
280
281
282
      }


283
      void setSolverWithLu(KSP ksp,
284
			  const char* kspPrefix,
285
286
			  KSPType kspType,
			  PCType pcType,
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
			  const MatSolverPackage matSolverPackage,
			  PetscReal rtol,
			  PetscReal atol,
			  PetscInt maxIt)
      {
	KSPSetType(ksp, kspType);
	KSPSetTolerances(ksp, rtol, atol, PETSC_DEFAULT, maxIt);
	//if (kspPrefix != "")
	if (strcmp(kspPrefix, "") != 0)
	  KSPSetOptionsPrefix(ksp, kspPrefix);
	KSPSetFromOptions(ksp);

	PC pc;
	KSPGetPC(ksp, &pc);
	PCSetType(pc, pcType);
	if (matSolverPackage != PETSC_NULL)
	  PCFactorSetMatSolverPackage(pc, matSolverPackage);
	PCSetFromOptions(pc);
305

306
#ifndef NDEBUG
307
308
309
310
311
312
313
    MSG("PetscOptionsView:\n");
    PetscViewer viewer;
    PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
    PetscViewerSetType(viewer, PETSCVIEWERASCII);
    PetscOptionsView(viewer);
    PetscViewerDestroy(&viewer);
#endif
314
      }
315

316

317
      void setSolver(KSP ksp,
318
		    const char* kspPrefix,
319
320
		    KSPType kspType,
		    PCType pcType,
321
322
323
324
		    PetscReal rtol,
		    PetscReal atol,
		    PetscInt maxIt)
      {
325

326
	setSolverWithLu(ksp, kspPrefix, kspType, pcType, PETSC_NULL,
327
328
			rtol, atol, maxIt);
      }
329
330


331
332
333
334
335
336
337
338
339
340
341
      void createSolver(MPI::Intracomm comm, KSP &ksp, Mat m, std::string kspPrefix, int info)
      {
	KSPCreate(comm, &ksp);
	#if (PETSC_VERSION_MINOR >= 5)
	  KSPSetOperators(ksp, m, m);
	#else
	  KSPSetOperators(ksp, m, m, SAME_NONZERO_PATTERN);
	#endif
	KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
	KSPSetType(ksp, KSPBCGS);
	KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
342

343
344
345
346
347
	if (info >= 10)
	  KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, PETSC_NULL);
	else if (info >= 20)
	  KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
      }
348

349
350
351
    } // end namespace petsc_helper
  } // end namespace Parallel
} // end namespace AMDiS