Actual source code: readinvecs.c

      2: /*    Reads in PETSc vectors from a PETSc binary file into matlab

      4:   Since this is called from Matlab it cannot be compiled with C++.
      5:   Modified Sept 28, 2003 RFK: updated obsolete mx functions.
      6: */

     9:  #include petscsys.h
     10:  #include petscvec.h
     11: #include "mex.h"
     12: #include <fcntl.h>
     13: #if defined(PETSC_HAVE_UNISTD_H)
     14: #include <unistd.h>
     15: #endif
     16: #if defined (PETSC_HAVE_IO_H)
     17: #include <io.h>
     18: #endif
     19: #if defined(PETSC_HAVE_STRINGS_H)
     20: #include <strings.h>
     21: #endif
     22: #if defined(PETSC_HAVE_STRING_H)
     23: #include <string.h>
     24: #endif
     25: #if defined(PETSC_HAVE_STROPTS_H)
     26: #include <stropts.h>
     27: #endif

     29: #define PETSC_MEX_ERROR(a) {fprintf(stdout,"ReadInVecs %s \n",a); return -1;}
     30: /*-----------------------------------------------------------------*/
     31: /*
     32:        Reads in a single vector
     33: */
     36: PetscErrorCode ReadInVecs(mxArray *plhs[],int t,int dim,int *dims)
     37: {
     38:   int    cookie = 0,M,compx = 0,i;
     39:
     40:   /* get size of matrix */
     41:   if (PetscBinaryRead(t,&cookie,1,PETSC_INT))   return -1;  /* finished reading file */
     42:   if (cookie != VEC_FILE_COOKIE) PETSC_MEX_ERROR("could not read vector cookie");
     43:   if (PetscBinaryRead(t,&M,1,PETSC_INT))        PETSC_MEX_ERROR("reading number rows");
     44:
     45:   if (dim == 1) {
     46:     plhs[0]  = mxCreateDoubleMatrix(M,1,mxREAL);
     47:   } else if (dim == 2) {
     48:     if (dims[0]*dims[1] != M) {
     49:       printf("PETSC_MEX_ERROR: m %d * n %d != M %d\n",dims[0],dims[1],M);
     50:       return -1;
     51:     }
     52:     plhs[0]  = mxCreateDoubleMatrix(dims[0],dims[1],mxREAL);
     53:   } else {
     54:     plhs[0] = mxCreateNumericArray(dim,dims,mxDOUBLE_CLASS,mxREAL);
     55:   }

     57:   /* read in matrix */
     58:   if (!compx) { /* real */
     59:     if (PetscBinaryRead(t,mxGetPr(plhs[0]),M,PETSC_DOUBLE)) PETSC_MEX_ERROR("read dense matrix");
     60:   } else { /* complex, currently not used */
     61:     for (i=0; i<M; i++) {
     62:       if (PetscBinaryRead(t,mxGetPr(plhs[0])+i,1,PETSC_DOUBLE)) PETSC_MEX_ERROR("read dense matrix");
     63:       if (PetscBinaryRead(t,mxGetPi(plhs[0])+i,1,PETSC_DOUBLE)) PETSC_MEX_ERROR("read dense matrix");
     64:     }
     65:   }
     66:   return 0;
     67: }

     69: #undef PETSC_MEX_ERROR
     70: #define PETSC_MEX_ERROR(a) {fprintf(stdout,"ReadInVecs %s \n",a); return;}
     71: /*-----------------------------------------------------------------*/

     75: void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
     76: {
     77:   static int fd = -1,dims[4],dim = 1,dof;
     78:   char       filename[1024],buffer[1024];
     79:   int        err,d2,d3,d4;
     80:   FILE       *file;

     82:   /* check output parameters */
     83:   if (nlhs != 1) PETSC_MEX_ERROR("Receive requires one output argument.");
     84:   if (fd == -1) {
     85:     if (!mxIsChar(prhs[0])) PETSC_MEX_ERROR("First arg must be string.");
     86:
     87:     /* open the file */
     88:     mxGetString(prhs[0],filename,256);
     89:     fd = open(filename,O_RDONLY,0);

     91:     strcat(filename,".info");
     92:     file = fopen(filename,"r");
     93:     if (file) {
     94:       fgets(buffer,1024,file);
     95:       if (!strncmp(buffer,"-daload_info",12)) {
     96:         sscanf(buffer,"-daload_info %d,%d,%d,%d,%d,%d,%d,%d\n",&dim,&dims[0],&dims[1],&dims[2],&dof,&d2,&d3,&d4);
     97:         if (dof > 1) {
     98:           dim++;
     99:           dims[3] = dims[2];
    100:           dims[2] = dims[1];
    101:           dims[1] = dims[0];
    102:           dims[0] = dof;
    103:         }
    104:       }
    105:       fclose(file);
    106:     }
    107:   }

    109:   /* read in the next vector */
    110:   err = ReadInVecs(plhs,fd,dim,dims);

    112:   if (err) {  /* file is finished so close and allow a restart */
    113:     close(fd);
    114:     fd = -1;
    115:   }
    116:   return;
    117: }

    120: