/***********************************************************

      SAVEMAT.C

      A routine for saving a matrix as a MATfile, the
      native format of Matlab.
      The routine can handle matrices in double precision,
      floating point, long integer, and signed and unsigned
      short integer.
      It will also convert a C nul-terminated string into
      a Matlab string, which is stored as a double precision
      array.

      This is a major modification of a function written by
      J. N. Little, 86-11-03, for MATLAB.

            Eric Firing
            88-02-12

      Mon  03-13-1989  Modified to put #defines into the header.

      Tue  09-19-89 JR Modified to store arrays in either column-wise
                       or row-wise fashion, depending on the
                       setting of the tens digit of the imagf flag:
                       choose from REAL_COLS, REAL_ROWS, IMAG_COLS,
                       IMAG_ROWS.
           93/02/22 JR Modified to actually transpose any row-wise
		       arrays to suit Matlab 4.0, which requires
		       all arrays to be column-wise, regardless of
		       the orientation bit in the file header.

**************************************************************/

#include "common.h"
#include "matfile.h"


/**** Convert a nul-terminated string to a Matlab string.
      Allocate heap space for the Matlab string.
****/
double *matstring(char *s, int *n)
{
   int i;
   double *mstring;

   *n = strlen(s);
   mstring = (double *) calloc(*n, sizeof(double));
   if (mstring != NULL)  {
      for (i=0; i<*n; i++)
         mstring[i] = (double) s[i];
   }
   return(mstring);  /* will be NULL if there is an error */
}                       /* matstring() */

void savemat(FILE *fp,
             char type, char *pname,
             int mrows, int ncols, int imagf,
             char *preal, char *pimag)
{                        /* savemat() */
   Fmatrix x;
   int mn;
   int nbytes, free_preal = 0, free_pimag = 0;
   char *mstring;  /* temporary pointer */

   x.type = 0;
   switch (type) {
      case 't':
         /* Note that preal is changed to point to a new
            block of memory allocated on the heap.
         */
         mstring = (char *) matstring(preal, &ncols);
         if (mstring == NULL)  {
            printf("\nSAVEMAT error, not enough memory for MAT string.\n");
            exit(-1);
         }
         preal = mstring;
	 free_preal = 1;
         mrows = 1;
         nbytes = 8;
         x.type += 1;
         break;
      case 'd':
         nbytes = 8;
         x.type += 0;
         break;
      case 'f':
         nbytes = 4;
         x.type += 10;
         break;
      case 'l':
         nbytes = 4;
         x.type += 20;
         break;
      case 'i':
         nbytes = 4;
         x.type += 20;
         break;
      case 's':
         nbytes = 2;
         x.type += 30;
         break;
      case 'u':
         nbytes = 2;
         x.type += 40;
         break;
      default: { printf("\nSAVEMAT error, invalid type: %c\n ", type);
         exit(-1);
      }
   }
   x.mrows = mrows;
   x.ncols = ncols;
   x.imagf = imagf % 10;            /* LS Bit is imaginary flag */
   x.namlen = strlen(pname) + 1;
   mn = x.mrows * x.ncols * nbytes;

   x.type += COLUMN_WISE;  /* Matlab 4.0 supports only column-wise! */
   if ( (imagf == REAL_ROWS || imagf == IMAG_ROWS) && (x.mrows > 1 || x.ncols > 1) )
   {
      if ( (preal = transpose(preal, mrows, ncols, nbytes)) == NULL ) exit(-1);
      free_preal = 1;
      if (x.imagf)
      {
         if ( (pimag = transpose(pimag, mrows, ncols, nbytes)) == NULL ) exit(-1);
	 free_pimag = 1;
      }
   }

   x.type += MACHINE_TYPE;

   fwrite(&x, sizeof(Fmatrix), 1, fp);
   fwrite(pname, sizeof(char), (int)x.namlen, fp);
   fwrite(preal, 1, mn, fp);
   if (free_preal) free(preal);
   if (x.imagf)
      fwrite(pimag, 1, mn, fp);
   if (free_pimag) free(pimag);
}                       /* savemat() */

#define Mrows Ccols
#define Mcols Crows

char *transpose(char pdata[], int Crows, int Ccols, int nbytes)
{
   int i, j, k, Cbytes_per_row, Mbytes_per_row;
   char *buff;

   if ( (buff = malloc(Crows * Ccols * nbytes)) == NULL )
   {
      printf("\n savemat transpose error: insufficient memory\n\n");
      return(NULL);
   }

   Cbytes_per_row = nbytes * Ccols;
   Mbytes_per_row = nbytes * Mcols;

   for (i = 0; i < Mrows; i++)
      for (j = 0; j < Mcols; j++)
	 for (k = 0; k < nbytes; k++)
	    buff[i * Mbytes_per_row + j * nbytes + k] =
           pdata[j * Cbytes_per_row + i * nbytes + k];

   return(buff);
}