PetscHelper.cc 9.68 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
			  const petsc::MatSolverPackage_t matSolverPackage,
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
			  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)
303
	  petsc::pc_factor_set_mat_solver_package(pc, matSolverPackage);
304
	PCSetFromOptions(pc);
305

306
#ifndef NDEBUG
307
308
309
310
311
312
        MSG("PetscOptionsView:\n");
        PetscViewer viewer;
        PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
        PetscViewerSetType(viewer, PETSCVIEWERASCII);
        petsc::options_view(viewer);
        PetscViewerDestroy(&viewer);
313
#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
      void createSolver(MPI::Intracomm comm, KSP &ksp, Mat m, std::string kspPrefix, int info)
      {
	KSPCreate(comm, &ksp);
333
        petsc::ksp_set_operators(ksp, m, m);
334
335
336
	KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
	KSPSetType(ksp, KSPBCGS);
	KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
337

338
	if (info >= 10)
339
	  petsc::ksp_monitor_set(ksp, KSPMonitorDefault);
340
	else if (info >= 20)
341
	  petsc::ksp_monitor_set(ksp, KSPMonitorTrueResidualNorm);
342
      }
343

344
345
346
    } // end namespace petsc_helper
  } // end namespace Parallel
} // end namespace AMDiS