Strange behavior with MatLoad

General Discussion
1

Dear FreeFem users,

I would like to use ObjectView and MatLoad for saving the PETSc numbering solution (vectors) in binary, possible to dome post-processing in Matlab, etc. However, when I am loading the results back in FreeFem, I get a very stange behavior: the unknown numbering change per process. The following MWE demonstrates this behavior:

load "PETSc"

{
int n=1000;
real[int] rhs(n+mpirank*n*0.1);

real[int, int] rhsView(rhs.n,1);
rhsView(:,0) = rhs;
ObjectView(rhsView, format = "binary", name = "rhsout");
cout << "writing: rank " << mpirank << "/" << mpisize << ", rhsView.n=" << rhsView.n << ", rhs.n=" << rhs.n << endl;
}

{
real[int, int] rhsView(1, 1);
MatLoad(rhsView, format = "binary", name = "rhsout");
real[int] rhs = rhsView.asarray;  
cout << "reading: rank " << mpirank << "/" << mpisize << ", rhsView.n=" << rhsView.n << ", rhs.n=" << rhs.n << endl;
}

The resulting output is the following:

writing: rank 1/2, rhsView.n=1100, rhs.n=1100
writing: rank 0/2, rhsView.n=1000, rhs.n=1000
reading: rank 1/2, rhsView.n=1050, rhs.n=1050
reading: rank 0/2, rhsView.n=1050, rhs.n=1050

I would expect the numerings to be the same. Is this a bug, or am I misunderstanding something?

2

PETSc does not store the partitioning. This has the big advantage that you can save with, e.g., 1 process, and then load with, e.g., 4 processes. In plain PETSc, you can force the local dimension before loading a Mat so that the partitioning suits your needs, but this is not implemented in the FreeFEM plugin. I could add that if you want, though, I’m not sure when (probably not before the end of this weekend).

3

Thanks for the quick reply, this makes complete sense.

I think there is no need for you implement additional features, I think I can just (i) save the vector, (ii) save the numbering accessed via GlobalNumbering (along with the corresponding ranks), (iii) load the vector and and the numbering on rank 0, and (iv) distribute the vector myself among the processes. Please correct me if I am wrong.

A note: the example FreeFem-sources/diffusion-substructuring-2d-PETSc.edp at master · FreeFem/FreeFem-sources · GitHub throws me an error at line 45 I think due to incorrect indexing of numbering[i]. I am not sure what the example exactly does, so I cannot propose a fix unfortunately.

4

I was trying to implement what I described above, however, it did not work. The error can be reproduced in the MWE I posted initially: I would expect it to work with reading only on a single core:

if(mpirank==0){
real[int, int] rhsView(1, 1);
MatLoad(rhsView, format = "binary", name = "rhsout");
real[int] rhs = rhsView.asarray;  
cout << "reading: rank " << mpirank << "/" << mpisize << ", rhsView.n=" << rhsView.n << ", rhs.n=" << rhs.n << endl;
}

However, at MatLoad the script seems to freeze. Is this behavior intended, or is this a bug?
(Again, I think I can circumvent this, I am just curious and I would like to understand the behavior)

5

You probably need to specify communicator = mpiCommSelf in the MatLoad().

6

Thanks, now it works perfectly!

7

Hi, @aszaboa

I have recently been learning parallel programming with FreeFEM. Could you elaborate on how the method you mentioned actually works? Do you have any script code for the implementation examples of your (i)-(iv)? I would like to learn from it.

Best regards.

8

Hi,

This is the macro I use to save the numbering of the processes:

    macro saveGlobal2LocalNumbering(infespace, inMat, inmeshName){

        /* just a dummy vector */
        real[int] tempPETSc;
        infespace gendef(uTemp);
        ChangeNumbering(inMat, uTemp[], tempPETSc);
        /* for some reason PETSc can only export (n,1) matrices ... */
        real[int, int] numExport(tempPETSc.n, 1);
        real[int, int] numExportFF(uTemp[].n, 1);

        /* mpiranks for the global indices */
        int[int] index2mpirank(tempPETSc.n);
        index2mpirank = mpirank;
        /* exporting ranks */
        for [i,  ranki : index2mpirank] numExport(i,0) = ranki;
        ObjectView(numExport, format = "binary", name = (wdir + "/" + inmeshName + "_mpiranks_PETSc"));

        /* exporting vector space index 2 rank */
        uTemp[] = mpirank;
        numExportFF(:,0)= uTemp[];
        ObjectView(numExportFF, format = "binary", name = (wdir + "/" + inmeshName + "_mpiranks_FE"));

    }//EOM

This is the macro to load the numbering:

    macro loadGlobal2LocalNumbering(vLocal2Global,inName,nGlobal){
        int[int] rankNumbering;
        /* loading the data */
        if(mpirank==0) {
            real[int, int] vecView(1, 1);
            MatLoad(vecView, format = "binary", name =  inName, communicator = mpiCommSelf);
            nGlobal = vecView.n;
            rankNumbering.resize(nGlobal);
            for [i, numi:rankNumbering] numi = vecView(i,0);
        }
        
        /* broadcast and resize numbering */
        broadcast(processor(0),nGlobal);
        if(mpirank!=0) 
            rankNumbering.resize(nGlobal);
        mpiBarrier(mpiCommWorld);
        broadcast(processor(0),rankNumbering);
        
        /* gettig the sizes of the local numberings */
        int nLocal = 0;
        for(int i=0; i<nGlobal; i++) {
            if(rankNumbering[i] == mpirank)
                nLocal++;
        }

        /* saving the actual global2local numberings */
        vLocal2Global.resize(nLocal);
        nLocal=-1;
        for(int i=0; i<nGlobal; i++) {
            if(rankNumbering[i] == mpirank)
                vLocal2Global[++nLocal] = i;
        }
            
        mpiBarrier(mpiCommWorld);
    } //EOM

These are the macros for writing and reading vectors:

    macro exportPETScBinary(vecPETSc, outputFileName) {
        real[int, int] solallView(vecPETSc.n, 1);
        solallView(:, 0) = vecPETSc;
        ObjectView(solallView, format = "binary", name = outputFileName);
    } //EOM

    macro loadPETScBinary(vec,loc2GlobNum,nGlob,fname){

        real[int] vecGlobIn(nGlob);
        if(mpirank==0) {
            real[int, int] vecGlobView(1, 1);
            MatLoad(vecGlobView, format = "binary", name =  fname, communicator = mpiCommSelf);
            vecGlobIn = vecGlobView.asarray;
        }

        broadcast(processor(0), vecGlobIn);
            for(int i=0; i<loc2GlobNum.n; i++)
                vec[i] = vecGlobIn[loc2GlobNum[i]];

    } //EOM

I hope that these are helpful. Let me know if you have questions!

9

Hi @aszaboa , Thank you!
That’s an excellent idea and the logic is very clear, it has been a great help to me.

By the way, I’d like to ask if you know whether there is a way to convert a PETSc Mat to a FreeFEM matrix. I know the reverse conversion is easy—simply use Mat A(B).

10

I am happy to hear that the scripts are helpful. I have no idea if there is a simple method to convert a PETSc Mat into a FreeFEM matrix. I guess you could approach a PETSc mat as a bundle of vectors, and do the conversion that way - I believe PETSc partitions the matrices like that, i.e., each process get a certain number of rows of a matrix. I am not sure whether the sparse matrix structure complicates the approach.

11

I greatly appreciate your detailed reply. I will continue to follow any progress in this regard and will share my feedback here as appropriate.