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