source: issm/trunk-jpl/src/wrappers/matlab/io/MatlabMatrixToPetscMat.cpp@ 13749

/* \file MatlabMatrixToPetscMatrix.cpp
 * \brief: convert a sparse or dense matlab matrix to a serial Petsc matrix:
 */

#ifdef HAVE_CONFIG_H
        #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
#include "../../c/shared/shared.h"

/*Petsc includes: */
#include <petscmat.h>
#include <petscvec.h>
#include <petscksp.h>

/*Matlab includes: */
#include "./matlabio.h"

PetscMat* MatlabMatrixToPetscMat(const mxArray* mxmatrix){

        int dummy;
        PetscMat* matrix=new PetscMat();

        MatlabMatrixToPetscMat(&matrix->matrix, &dummy, &dummy, mxmatrix);

        return matrix;
}
int MatlabMatrixToPetscMat(Mat* pmatrix,int* pmatrix_rows,int* pmatrix_cols,const mxArray* mxmatrix){

        int rows, cols;
        double *mxmatrix_ptr = NULL;
        double *tmatrix      = NULL;
        int ierr;
        int i,j;

        /*output: */
        Mat matrix = NULL;

        /*matlab indices: */
        mwIndex *ir = NULL;
        mwIndex *jc = NULL;
        double  *pr = NULL;
        int     count;
        int     nnz;
        int     nz;

        /*petsc indices: */
        int *idxm = NULL;
        int *idxn = NULL;

        /*Ok, first check if we are dealing with a sparse or full matrix: */
        if (mxIsSparse(mxmatrix)){

                /*Dealing with sparse matrix: recover size first: */
                mxmatrix_ptr=(double*)mxGetPr(mxmatrix);
                rows=mxGetM(mxmatrix);
                cols=mxGetN(mxmatrix);
                nnz=mxGetNzmax(mxmatrix);
                if(rows){
                        nz=(int)((double)nnz/(double)rows);
                }
                else{
                        nz=0;
                }

                ierr=MatCreateSeqAIJ(PETSC_COMM_SELF,rows,cols,nz,PETSC_NULL,&matrix);CHKERRQ(ierr);

                /*Now, get ir,jc and pr: */
                pr=mxGetPr(mxmatrix);
                ir=mxGetIr(mxmatrix);
                jc=mxGetJc(mxmatrix);

                /*Now, start inserting data into sparse matrix: */
                count=0;
                for(i=0;i<cols;i++){
                        for(j=0;j<(jc[i+1]-jc[i]);j++){
                                MatSetValue(matrix,ir[count],i,pr[count],INSERT_VALUES);
                                count++;
                        }
                }
        }
        else{
                /*Dealing with dense matrix: recover pointer and size: */
                mxmatrix_ptr=(double*)mxGetPr(mxmatrix);
                rows=mxGetM(mxmatrix);
                cols=mxGetN(mxmatrix);

                /*transpose, as Petsc now does not allows MAT_COLUMN_ORIENTED matrices in MatSetValues: */
                tmatrix=xNew<double>(rows*cols);
                for(i=0;i<cols;i++){
                        for(j=0;j<rows;j++){
                                *(tmatrix+rows*i+j)=*(mxmatrix_ptr+cols*j+i);
                        }
                }

                /*Create serial matrix: */
                ierr=MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&matrix);CHKERRQ(ierr);

                /*Insert mxmatrix_ptr values into petsc matrix: */
                idxm=xNew<int>(rows);
                idxn=xNew<int>(cols);

                for(i=0;i<rows;i++)idxm[i]=i;
                for(i=0;i<cols;i++)idxn[i]=i;

                ierr=MatSetValues(matrix,rows,idxm,cols,idxn,tmatrix,INSERT_VALUES); CHKERRQ(ierr);

                xDelete<double>(tmatrix);
        }

        /*Assemble matrix: */
        MatAssemblyBegin(matrix,MAT_FINAL_ASSEMBLY); 
        MatAssemblyEnd(matrix,MAT_FINAL_ASSEMBLY);

        /*Assign output pointer: */
        *pmatrix=matrix;
        if(pmatrix_rows) *pmatrix_rows=rows;
        if(pmatrix_cols) *pmatrix_cols=cols;

        return 1;
}