[petsc-users] MatScale returns different results depending on matrix size

[Roland Richter]

Hei,

I would like to scale a given matrix with a fixed scalar value, and
therefore would like to use MatScale(). Nevertheless, I observed an
interesting behavior depending on the size of the matrix, and currently
I am not sure why.

When running the attached code, I intend to divide all elements in the
matrix by a constant factor of 10. If I have three or fewer rows and
1024 columns, I get the expected result. If I have four or more rows
(with the same number of columns), suddenly my scaling factor seems to
be 0.01 instead of 0.1 for the PETSc-matrix. The armadillo-based matrix
still behaves as expected.

I currently do not understand that behavior, but do not see any problems
with the code either. Are there any possible explanations for that behavior?

Thank you very much,

regards,

Roland Richter


#include <iostream>

#include <armadillo>
#include <petscmat.h>
#include <petscvec.h>
#include <fftw3.h>
#include <fftw3-mpi.h>
#include <metis.h>
#include <chrono>

void retrieve_data_from_PETSc(const Mat petsc_mat, arma::cx_mat &out_data,
							  const arma::uword Ntime, const arma::uword Nradius) {
	if(out_data.n_rows != Ntime && out_data.n_cols != Nradius) {
		out_data = arma::zeros<arma::cx_mat>(Ntime, Nradius);
	}
	Mat local_mat;
	arma::Col<int> vector_indices_radius = arma::linspace<arma::Col<int>>(0, Nradius - 1, Nradius);
	arma::Col<int> vector_indices_time = arma::linspace<arma::Col<int>>(0, Ntime - 1, Ntime);
	MatCreateRedundantMatrix(petsc_mat, Ntime * Nradius, MPI_COMM_NULL, MAT_INITIAL_MATRIX, &local_mat);
	MatAssemblyBegin(local_mat, MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(local_mat, MAT_FINAL_ASSEMBLY);
	MatGetValues(local_mat, Nradius, vector_indices_radius.memptr(), Ntime, vector_indices_time.memptr(), out_data.memptr());
	MatDestroy(&local_mat);
	out_data = out_data.st();
}

void store_data_in_PETSc(const arma::cx_mat &in_data, Mat &petsc_mat) {
	const arma::uword Ntime = in_data.n_cols;
	const arma::uword Nradius = in_data.n_rows;
	arma::Col<int> vector_indices_radius = arma::linspace<arma::Col<int>>(0, Nradius - 1, Nradius);
	arma::Col<int> vector_indices_time = arma::linspace<arma::Col<int>>(0, Ntime - 1, Ntime);
	MatZeroEntries(petsc_mat);
	MatAssemblyBegin(petsc_mat, MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(petsc_mat, MAT_FINAL_ASSEMBLY);
	arma::cx_mat local_mat = in_data.st();
	MatSetValues(petsc_mat, Nradius, vector_indices_radius.memptr(), Ntime, vector_indices_time.memptr(), local_mat.memptr(), INSERT_VALUES);
	MatAssemblyBegin(petsc_mat, MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(petsc_mat, MAT_FINAL_ASSEMBLY);
}

void test_scaling_arma(const arma::cx_mat & in_mat,
					   arma::cx_mat &out_mat,
					   const arma::cx_double &scaling_factor) {
	out_mat = in_mat;
	out_mat *= scaling_factor;
}

void test_scaling_petsc(const Mat &in_mat,
						Mat &out_mat,
						const PetscScalar &scaling_factor) {
	MatCopy (in_mat, out_mat, SAME_NONZERO_PATTERN);
	MatScale (out_mat, scaling_factor);
}

void test_scaling(const size_t matrix_size_rows, const size_t matrix_size_cols, const bool print_matrices) {
	arma::cx_mat in_mat = arma::zeros<arma::cx_mat>(matrix_size_rows, matrix_size_cols),
			out_mat = arma::zeros<arma::cx_mat>(matrix_size_rows, matrix_size_cols);
	arma::cx_rowvec matrix_vec = arma::conv_to<arma::cx_rowvec>::from(arma::linspace<arma::rowvec>(0, matrix_size_cols - 1, matrix_size_cols));
	in_mat.each_row([&](arma::cx_rowvec &a){
		a = matrix_vec;
	});

	Mat petsc_in_mat, petsc_out_mat;
	arma::cx_mat petsc_out_comparison_mat = arma::zeros<arma::cx_mat>(matrix_size_rows, matrix_size_cols);
	MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, matrix_size_rows, matrix_size_cols, NULL, &petsc_in_mat);
	MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, matrix_size_rows, matrix_size_cols, NULL, &petsc_out_mat);

	store_data_in_PETSc (in_mat, petsc_in_mat);
	store_data_in_PETSc (petsc_out_comparison_mat, petsc_out_mat);

	const std::complex<double> scaling_factor{1. / 10, 0.};

	test_scaling_arma (in_mat, out_mat, scaling_factor);
	test_scaling_petsc (petsc_in_mat, petsc_out_mat, scaling_factor);

	retrieve_data_from_PETSc (petsc_out_mat, petsc_out_comparison_mat, matrix_size_cols, matrix_size_rows);

	if(print_matrices) {
		std::cout << "In-matrix, ARMA:\n" << in_mat
				  << "\n\nOut-matrix, ARMA:\n" << out_mat
				  << "\n\nComparison-out-matrix, ARMA:\n" << petsc_out_comparison_mat
				  << "\n\nDifference: \n" << arma::abs(petsc_out_comparison_mat - out_mat)
				  <<'\n';
	}
	std::cout << "Arma and PETSc are equal: " << arma::approx_equal(out_mat, petsc_out_comparison_mat, "reldiff", 1e-8) << '\n';

	MatDestroy (&petsc_in_mat);
	MatDestroy (&petsc_out_mat);
}

int main(int argc, char **args) {

	PetscMPIInt rank, size;
	PetscInitialize(&argc, &args, (char*) 0, NULL);

	MPI_Comm_size(PETSC_COMM_WORLD, &size);
	MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

	test_scaling (3, 1024, false);
	test_scaling (4, 1024, false);

	PetscFinalize();
	return 0;
}