PetscHelper.cc 9.66 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
        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