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

      SAVEMATC.C

      This is a set of routines for building arrays in
      matfiles, a column or row at a time.
      The routine can handle matrices in double precision,
      floating point, long integer, and signed and unsigned
      short integer.  Only real number types are permitted,
      since complex arrays are stored as separate real and
      imaginary parts; this could be done a column at a time
      only if the number of columns were known beforehand,
      whereas this set of routines is mainly designed for
      the case in which one does not know how many columns
      will be added.

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

            Eric Firing
            Thu  03-09-1989

     9/21/89 - JR - added function tonan_f, which is identical
                    to tonan_d, except the source is a floating
                    point array;
                  - for both tonan_f and tonan_d, changed the
                    comparison of source against ADJ_BADFLOAT
                    from <= to < (shouldn't really matter)
     7/20/94 - JR - start_mat_array modified to print error message
		    and exit if row-oriented flag is passed as argument;
		    Matlab version 4 no longer supports this
**************************************************************/

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

#if PROTOTYPE_ALLOWED
int convert_type(char type, int machine, int *nbytes, long *P);
#else
int convert_type();
#endif

#define ADJ_BADFLOAT 0.999999e37

#if PROTOTYPE_ALLOWED
MAT_ARRAY_FP_TYPE *start_mat_array(FILE *fp,
                                   char type, char *pname, char rc_flag,
                                   int n_rc)
#else
MAT_ARRAY_FP_TYPE *start_mat_array(fp, type, pname, rc_flag, n_rc)
FILE *fp;       /* File pointer */
char type;      /* d, f, l, i, u, t:
                   double, float, long, short, ushort, text */
char *pname;    /* pointer to matrix name */
char rc_flag;   /* r if rows will be added, c otherwise */
int n_rc;       /* number of columns of rows will be added, or the converse */
#endif
{
   Fmatrix x;
   MAT_ARRAY_FP_TYPE *mat_ptr;

   mat_ptr = (MAT_ARRAY_FP_TYPE *)malloc(sizeof(MAT_ARRAY_FP_TYPE));
   if (mat_ptr == NULL)
   {
      printf("\nSAVEMATC error, no memory for MAT_ARRAY_FP_TYPE\n");
      exit(-1);
   }
   mat_ptr->fp = fp;
   mat_ptr->Fmatrix_origin = ftell(fp);
   mat_ptr->count = 0;
   mat_ptr->n_rc = n_rc;

   if (convert_type(type, MACHINE_TYPE, &(mat_ptr->nbytes), &(x.type))) exit(-1);

   switch (rc_flag)
   {
      case 'r':
	 printf("\nSAVEMATC error, Matlab version 4 does not allow row-wise array storage\n");
	 exit(-1);
	 /*
         x.type += ROW_WISE;
         x.ncols = n_rc;
         x.mrows = 0;
         mat_ptr->count_offset = (char *)&(x.mrows) - (char *)&(x);*/
         break;
      case 'c':
         x.mrows = n_rc;
         x.ncols = 0;
         mat_ptr->count_offset = (char *)&(x.ncols) - (char *)&(x);
         break;
      default:
         printf("\nSAVEMATC error, invalid rc_flag: %c\n ", rc_flag);
         exit(-1);
   }

   x.type += MACHINE_TYPE;
   x.imagf = 0;
   x.namlen = strlen(pname) + 1;

   fwrite(&x, sizeof(Fmatrix), 1, fp);
   fwrite(pname, sizeof(char), (int)x.namlen, fp);

   return(mat_ptr);
}                       /* start_mat_array() */

#if PROTOTYPE_ALLOWED
MAT_ARRAY_FP_TYPE *continue_mat_array(FILE *fp,
                                      char type, char rc_flag,
                                      int n_rc)
#else
MAT_ARRAY_FP_TYPE *continue_mat_array(fp, type, rc_flag, n_rc)
FILE *fp;       /* File pointer */
char type;      /* d, f, l, i, u, t:
                   double, float, long, short, ushort, text */
char rc_flag;   /* r if rows will be added, c otherwise */
int n_rc;       /* number of columns in rows to be added, or the converse */
#endif
{
   MAT_ARRAY_FP_TYPE *mat_ptr;
   Fmatrix x;
   long p;

   if (fseek(fp, 0L, SEEK_SET))
   {
      fprintf(stderr, "\n ERROR: Seek to matfile beginning failed\n");
      return(NULL);
   }
   if (fread( (char *) &x, sizeof(Fmatrix), 1, fp) != 1 )
   {
      fprintf(stderr, "\n ERROR: Can't read Fmatrix in continue_mat_array\n");
      return(NULL);
   }
   if ( (mat_ptr = (MAT_ARRAY_FP_TYPE *) malloc(sizeof(MAT_ARRAY_FP_TYPE))) == NULL )
   {
      fprintf(stderr, "\n ERROR: Insufficent memory in continue_mat_array\n");
      return(NULL);
   }
   mat_ptr->fp = fp;
   mat_ptr->Fmatrix_origin = 0L;

   /*------------------------- check for compatibility ----------------------*/
   if ( type_M(x) != MACHINE_TYPE )
   {
      fprintf(stderr, "\n ERROR: Matfile written by different machine\n");
      return(NULL);
   }
   if ( type_O(x) == COLUMN_WISE )
   {
      if (rc_flag != 'c' || x.mrows != n_rc)
      {
         fprintf(stderr, "\n ERROR: Orientation/columnsize mismatch\n");
         return(NULL);
      }
      mat_ptr->count_offset = (char *) &(x.ncols) - (char *) &x;
      mat_ptr->count        = x.ncols;
   }
   else /* ROW_WISE */
   {
      if (rc_flag != 'r' || x.ncols != n_rc)
      {
         fprintf(stderr, "\n ERROR: Orientation/rowsize mismatch\n");
         return(NULL);
      }
      mat_ptr->count_offset = (char *) &(x.mrows) - (char *) &x;
      mat_ptr->count        = x.mrows;
   }
   if (convert_type(type, MACHINE_TYPE, &(mat_ptr->nbytes), &p)) return(NULL);
   if (type_P(x) != p)
   {
      fprintf(stderr, "\n ERROR: Variable type mismatch\n");
      return(NULL);
   }
   if (type_T(x))
   {
      fprintf(stderr, "\n ERROR: Matrix is of type text\n");
      return(NULL);
   }
   if (x.imagf)
   {
      fprintf(stderr, "\n ERROR: Matrix has imaginary component\n");
      return(NULL);
   }
   mat_ptr->n_rc = n_rc;
   if (fseek(fp, 0L, SEEK_END))
   {
      fprintf(stderr, "\n ERROR: Seeking to matfile end\n");
      return(NULL);
   }
   return(mat_ptr);
}

#if PROTOTYPE_ALLOWED
void add_to_mat_array(MAT_ARRAY_FP_TYPE *mat_ptr, char *data_ptr)
#else
void add_to_mat_array(mat_ptr, data_ptr)
MAT_ARRAY_FP_TYPE *mat_ptr;
char *data_ptr;
#endif
{
   fwrite(data_ptr, mat_ptr->nbytes, mat_ptr->n_rc, mat_ptr->fp);
   mat_ptr->count++;
}

#if PROTOTYPE_ALLOWED
void end_mat_array(MAT_ARRAY_FP_TYPE *mat_ptr)
#else
void end_mat_array(mat_ptr)
MAT_ARRAY_FP_TYPE *mat_ptr;
#endif
{
   fseek(mat_ptr->fp, mat_ptr->Fmatrix_origin + mat_ptr->count_offset, 0);
   fwrite(&(mat_ptr->count), sizeof(mat_ptr->count), 1, mat_ptr->fp);
   fseek(mat_ptr->fp, 0L, 2);
   free(mat_ptr);
}

#if PROTOTYPE_ALLOWED
void tonan_f(double *dest, float *source, int n)
#else
void tonan_f(dest, source, n)
double *dest;
float *source;
int n;
#endif
{
   int i;
   union { double d; long l[2]; } NaN;

   NaN.l[0] = NaN_d_0;
   NaN.l[1] = NaN_d_1;
   for (i=0; i<n; i++, dest++, source++)
   {
      *dest =  (*source < ADJ_BADFLOAT) ? ((double) *source) : NaN.d;
   }
}

#if PROTOTYPE_ALLOWED
void tonan_d(double *dest, double *source, int n)
#else
void tonan_d(dest, source, n)
double *dest, *source;
int n;
#endif
{
   int i;
   union { double d; long l[2]; } NaN;

   NaN.l[0] = NaN_d_0;
   NaN.l[1] = NaN_d_1;
   for (i=0; i<n; i++, dest++, source++)
   {
      *dest =  (*source < ADJ_BADFLOAT) ? *source : NaN.d;
   }
}

#if PROTOTYPE_ALLOWED
void tonan_mask_d(double *dest, double *source, int *mask, int minmask, int n)
#else
void tonan_mask_d(dest, source, mask, minmask, n)
double *dest, *source;
int *mask, minmask, n;
#endif
{
   int i;
   union { double d; long l[2]; } NaN;

   NaN.l[0] = NaN_d_0;
   NaN.l[1] = NaN_d_1;
   for (i=0; i<n; i++, dest++, source++, mask++)
   {
      *dest =  (*mask >= minmask) ? *source : NaN.d;
   }
}

/* Test routine. */
#ifdef EXE
double d_array[12][10];
#if PROTOTYPE_ALLOWED
void main(void)
#else
void main()
#endif
{
   int ncols;
   int mrows;
   int i;
   FILE *fp;
   MAT_ARRAY_FP_TYPE *mat_ptr;

   fp = fopen("testsm.mat", "wb");
   ncols = 10;
   mrows = 12;
   for (i=0; i<ncols*mrows; i++)
   {
      d_array[0][i] = i+1;
   }
   mat_ptr = start_mat_array(fp, 'd', "testrows", 'r', 10);
   for (i=0; i<mrows; i++)
      add_to_mat_array(mat_ptr, (char *)&(d_array[i][0]));
   end_mat_array(mat_ptr);
   mat_ptr = start_mat_array(fp, 'd', "testcols", 'c', 10);
   for (i=0; i<mrows; i++)
      add_to_mat_array(mat_ptr, (char *)&(d_array[i][0]));
   end_mat_array(mat_ptr);
   fclose(fp);
}

#endif

#if PROTOTYPE_ALLOWED
int convert_type(char type, int machine, int *nbytes, long *P)
#else
int convert_type(type, machine, nbytes, P)
char type;
int machine, *nbytes;
long *P;      /* to match with Matlab Fmatrix type (MOPT) */
#endif
{
   switch (type)
   {
      case 'd':  *nbytes = 8; *P = 00; break;
      case 'f':  *nbytes = 4; *P = 10; break;
      case 'l':  *nbytes = 4; *P = 20; break;
      case 'i':  if (HOST_BYTES_PER_WORD == 2) {*nbytes = 2; *P = 30;}
                 else {*nbytes = 4; *P = 20;}
                 break;
      case 's':  *nbytes = 2; *P = 30; break;
      case 'u':  *nbytes = 2; *P = 40; break;
      default :  fprintf(stderr, "\n ERROR: Invalid type %c\n ", type);
                 return(-1);
   }
   return(0);
}