1
$\begingroup$

I'm using Trilinos for a problem and I have an Epetra_MultiVector object with 1 vector of length m*n*nFrames. I need to turn it into an Epetra_MultiVector object with nFrames of length m*n. The first m*n entries need to go into the first vector, the second into the second vector and so on.

I have this working fine with just one process. What I'm doing is extracting all the data to a standard C++ array on process 0, reshaping that into a 2D standard C++ array (an array of pointers actually), and then creating a new Epetra_MultiVector with my 2D array.

With more than one process, however, because everything is on process 0, when I create the new Epetra_MultiVector, all the data for the other process is nonsense (accessing undefined memory, I guess).

So my question ultimately is: how do I reshape a Trilinos Epetra_MultiVector correctly with more than one process?

$\endgroup$
2
$\begingroup$

Reshaping arrays is something that's a bit Matlab-specific -- most other environments (including Trilinos, PETSc, ...) allow you to create a vector or matrix but not let you re-interpret it as something else. You need to copy the elements by hand into a matrix of different size.

$\endgroup$
1
$\begingroup$

I've written a function that does it. It does not do it in place, and is not particularly well written, but for posterity I'll add my solution. I realize that I really don't need to reshape the standard C++ array, but I'll leave it there until it becomes a problem for me.

int ReshapeMultiVector(Epetra_MultiVector *&reshapedMultiVec, const int length, const int numVecs, const Epetra_MultiVector &originalMultiVec, Epetra_MpiComm &comm){

    /* Function to reshape a multivector. Does it by copying the values, so additional memory will be required.
     * Note that this function does NOT delete the original vector.
     * INPUTS
     *  reshapedMultiVec - This is a reference to a null pointer that will point to the
     *             new (reshaped) multivector.
     *  length       - Global length of the new multivector
     *  numVecs      - The number of vectors in the new multivector.
     *  originalMultiVec - The original multivector that is being reshaped.
     *  comm         - MPI communicator
     * OUTPUT
     *  int      - Error flag
     */

    int err;

    int numVecsOrig = originalMultiVec.NumVectors();
    int lengthOrig = originalMultiVec.GlobalLength();

    //extract phase_gradx and phase_grady as standard arrays so that we can reshape them.
    double** originalMultiVec_std = new double*[numVecsOrig];
    for(int i=0;i<numVecsOrig;++i){
        originalMultiVec_std[i] = new double[lengthOrig];
    }
    err = MultiVectorToArray(originalMultiVec_std, originalMultiVec, false);
    if(err) return err;

    // Copy the values to all the processors after extracting them.
    for(int i=0;i<numVecsOrig;++i){
        comm.Broadcast(originalMultiVec_std[i], lengthOrig, 0);
    }


    double** reshaped_std = new double*[numVecs];
    for(int i=0;i<numVecs;++i){
        reshaped_std[i] = new double[length];
    }
    err = Reshape2DArray(originalMultiVec_std, lengthOrig, numVecsOrig, reshaped_std, length, numVecs);
    if(err) return err;


    Epetra_Map reshapedMap(length,0,comm);

    //Initialize the vector
    int numMyElements = reshapedMap.NumMyElements();
    int* myGlobalElements = reshapedMap.MyGlobalElements(); 


    reshapedMultiVec = new Epetra_MultiVector(reshapedMap,numVecs);
    double ** Ap = reshapedMultiVec->Pointers();

    for (int j=0; j<numVecs; ++j)
    {
        double * v = Ap[j];

        // Fill it
        for (int i=0; i<numMyElements; ++i)
        {
            v[i] = reshaped_std[j][myGlobalElements[i]];
        }
    }


    // free memory
    for(int i=0;i<numVecsOrig;++i){
        delete[] originalMultiVec_std[i];
    }
    delete[] originalMultiVec_std;
    for(int i=0;i<numVecs;++i){
        delete[] reshaped_std[i];
    }
    delete[] reshaped_std;


    return(EXIT_SUCCESS);
}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.