How to compute RARt with A and R as distributed (MPI) matrices ?
This works with sequential matrices.
The doc say "currently only implemented for pairs of AIJ matrices and classes
which inherit from AIJ": I supposed that MPIAIJ was someway inheriting from
AIJ, seems that it doesn't.
Is this kind of matrix product possible with distributed matrices in PETSc ? Or
is this a known limitation ?
Do I go the wrong way to do that (= should use another method) ? If yes, what
is the correct one ?
Franck
PS: running debian/testing + gcc-6.3 + bitbucket petsc.
>> mpirun -n 2 matRARt.exe seq
Mat Object: 1 MPI processes
type: seqaij
row 0: (0, 1.) (1, 0.)
row 1: (0, 0.) (1, 1.)
>> mpirun -n 2 matRARt.exe mpi
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: No support for this operation for this object type
[0]PETSC ERROR: Matrix of type <mpiaij> does not support RARt
// How to compute RARt with A and R as distributed (MPI) matrices ?
//
// ~> g++ -o matRARt.exe matRARt.cpp -lpetsc -lm; mpirun -n 2 matRARt.exe
#include <iostream>
#include <string>
#include "petsc.h"
using namespace std;
int main(int argc,char **argv) {
if (argc != 2) {cout << "error: need arg = seq or mpi" << endl; return 1;}
string matType = argv[1]; if (matType != "seq" && matType != "mpi") {cout << "error: need arg = seq or mpi" << endl; return 1;}
PetscInitialize(&argc, &argv, NULL, NULL);
int size = 0; MPI_Comm_size(MPI_COMM_WORLD, &size); if (size != 2) {cout << "error: np != 2" << endl; return 1;}
int rank = 0; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
PetscInt idx[2] = {0, 1}; PetscScalar val[4] = {1., 0., 0., 1.};
Mat R;
if(matType == "mpi") MatCreateAIJ(PETSC_COMM_WORLD, 1, 1, 2, 2, 1, NULL, 1, NULL, &R);
if(matType == "seq") MatCreateSeqAIJ(PETSC_COMM_SELF, 2, 2, 2, NULL, &R);
MatSetValues(R, 2, idx, 2, idx, val, ADD_VALUES);
MatAssemblyBegin(R, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(R, MAT_FINAL_ASSEMBLY);
Mat A;
if(matType == "mpi") MatCreateAIJ(PETSC_COMM_WORLD, 1, 1, 2, 2, 1, NULL, 1, NULL, &A);
if(matType == "seq") MatCreateSeqAIJ(PETSC_COMM_SELF, 2, 2, 2, NULL, &A);
MatSetValues(A, 2, idx, 2, idx, val, ADD_VALUES);
MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
Mat C; MatRARt(A, R, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C);
if(matType == "mpi") MatView(C, PETSC_VIEWER_STDOUT_WORLD);
if(matType == "seq" && rank == 0) MatView(C, PETSC_VIEWER_STDOUT_SELF);
PetscFinalize();
return 0;
}
