// Remove zeros from the matrix Ain and get matrix Aout // ierr = dropZero(Ain, Aout, 1e-13); // Then work with matrix Aout as with usual petsc matrix int dropZero(Mat& Ain, Mat& Aout, double droptol) { int i, j; int numbercols; int nrows; // first count the nonzero elements PetscFunctionBegin; ierr = MatGetSize(Ain, &nrows, &numbercols); CHKERRA(ierr); int* number_of_nbs = new int[nrows]; const int *cols; int ncols; const PetscScalar *vals; cout << "Size of matrix is " << nrows << " rows and " << numbercols << " columns\n"; int counter = 0; int nocounter = 0; for (i = 0; i < nrows; i++) { ierr = MatGetRow(Ain, i, &ncols, &cols, &vals); CHKERRA(ierr); int loccounter = 0; // ncols - number of nonzero el-s in the row with number i // cols - number of the column // vals - nonzero values for (j = 0; j < ncols; j++) { if (fabs(vals[j]) > droptol) { counter++; loccounter++; } else { nocounter++; } } number_of_nbs[i] = loccounter; ierr = MatRestoreRow(Ain, i, &ncols, &cols, &vals); CHKERRA(ierr); } cout << "Found " << counter << " nonzero and " << nocounter << " zero elements ( < 1e-12)\n" << flush; ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nrows, numbercols, 0, number_of_nbs, &Aout); CHKERRA(ierr); delete[] number_of_nbs; Array1dReal resultElMat(50); Array1dInt rowIdx(1); Array1dInt colIdx(50); for (i = 0; i < nrows; i++) { ierr = MatGetRow(Ain, i, &ncols, &cols, &vals); CHKERRA(ierr); int loccounter = 0; rowIdx(0) = i; for (j = 0; j < ncols; j++) { if (fabs(vals[j]) > droptol) { colIdx(loccounter) = cols[j]; resultElMat(loccounter) = vals[j]; loccounter++; } } ierr = MatRestoreRow(Ain, i, &ncols, &cols, &vals); CHKERRA(ierr); ierr = MatSetValues(Aout, 1, &rowIdx(0), loccounter, &colIdx(0), &resultElMat(0), ADD_VALUES); } ierr = MatAssemblyBegin(Aout, MAT_FINAL_ASSEMBLY); CHKERRA(ierr); ierr = MatAssemblyEnd(Aout, MAT_FINAL_ASSEMBLY); CHKERRA(ierr); PetscFunctionReturn(0); return 0; }