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

      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"

int convert_type(char type, int machine, int *nbytes, int *P);

MAT_ARRAY_FP_TYPE *start_mat_array(FILE *fp,
                                   char type, char *pname, char rc_flag,
                                   int n_rc)
{
   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() */

MAT_ARRAY_FP_TYPE *continue_mat_array(FILE *fp,
                                      char type, char rc_flag,
                                      int n_rc)
{
   MAT_ARRAY_FP_TYPE *mat_ptr;
   Fmatrix x;
   int 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);
}

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

void end_mat_array(MAT_ARRAY_FP_TYPE *mat_ptr)
{
   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);
}

void tonan_f(double *dest, float *source, int n)
{
   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;
   }
}

void tonan_d(double *dest, double *source, int n)
{
   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;
   }
}

void tonan_mask_d(double *dest, double *source, int *mask, int minmask, int n)
{
   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];
int main(void)
{
   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);
   return 0;
}

#endif

int convert_type(char type, int machine, int *nbytes, int *P)
{
   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':  *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);
}