NetCDF4: added support for writting complex numbers (status: experimental).

ChangeLog

+2020-04-20  Uwe Schulzweida
+
+	* NetCDF4: added support for writting complex numbers (status: experimental)
+
 2020-04-08  Uwe Schulzweida
 
 	* Fix: --without-threads failed (use HAVE_LIBPTHREAD in async_worker.c)

src/cdf_int.c

@@ -429,8 +429,8 @@ void cdf_put_vara_double(int ncid, int varid, const size_t start[],
 }
 
 
-void  cdf_put_vara_float(int ncid, int varid, const size_t start[],
-                         const size_t count[], const float *fp)
+void cdf_put_vara_float(int ncid, int varid, const size_t start[],
+                        const size_t count[], const float *fp)
 {
   int status = nc_put_vara_float(ncid, varid, start, count, fp);
 
@@ -454,7 +454,18 @@ void  cdf_put_vara_float(int ncid, int varid, const size_t start[],
 }
 
 
-void  cdf_get_vara_int(int ncid, int varid, const size_t start[], const size_t count[], int *dp)
+void cdf_put_vara(int ncid, int varid, const size_t start[], const size_t count[], const void *cp)
+{
+  int status = nc_put_vara(ncid, varid, start, count, cp);
+
+  if ( CDF_Debug || status != NC_NOERR )
+    Message("ncid = %d varid = %d", ncid, varid);
+
+  if ( status != NC_NOERR ) Error("%s", nc_strerror(status));
+}
+
+
+void cdf_get_vara_int(int ncid, int varid, const size_t start[], const size_t count[], int *dp)
 {
   int status = nc_get_vara_int(ncid, varid, start, count, dp);
 
@@ -465,7 +476,7 @@ void  cdf_get_vara_int(int ncid, int varid, const size_t start[], const size_t c
 }
 
 
-void  cdf_get_vara_double(int ncid, int varid, const size_t start[], const size_t count[], double *dp)
+void cdf_get_vara_double(int ncid, int varid, const size_t start[], const size_t count[], double *dp)
 {
   int status = nc_get_vara_double(ncid, varid, start, count, dp);
 
@@ -476,7 +487,7 @@ void  cdf_get_vara_double(int ncid, int varid, const size_t start[], const size_
 }
 
 
-void  cdf_get_vara_float(int ncid, int varid, const size_t start[], const size_t count[], float *fp)
+void cdf_get_vara_float(int ncid, int varid, const size_t start[], const size_t count[], float *fp)
 {
   int status = nc_get_vara_float(ncid, varid, start, count, fp);
 
@@ -487,7 +498,7 @@ void  cdf_get_vara_float(int ncid, int varid, const size_t start[], const size_t
 }
 
 
-void  cdf_get_vara_text(int ncid, int varid, const size_t start[], const size_t count[], char *tp)
+void cdf_get_vara_text(int ncid, int varid, const size_t start[], const size_t count[], char *tp)
 {
   int status = nc_get_vara_text(ncid, varid, start, count, tp);
 
@@ -498,7 +509,7 @@ void  cdf_get_vara_text(int ncid, int varid, const size_t start[], const size_t
 }
 
 
-void  cdf_get_vara_uchar(int ncid, int varid, const size_t start[], const size_t count[], unsigned char *tp)
+void cdf_get_vara_uchar(int ncid, int varid, const size_t start[], const size_t count[], unsigned char *tp)
 {
   int status = nc_get_vara_uchar(ncid, varid, start, count, tp);
 

src/cdf_int.h

@@ -5,32 +5,32 @@
 
 #include <netcdf.h>
 
-void cdf_create (const char *path, int cmode, int *idp);
-int  cdf_open   (const char *path, int omode, int *idp);
-void cdf_close  (int ncid);
+void cdf_create(const char *path, int cmode, int *idp);
+int  cdf_open(const char *path, int omode, int *idp);
+void cdf_close(int ncid);
 
 void cdf_redef(int ncid);
 void cdf_enddef(int ncid);
 void cdf__enddef(const int ncid, const size_t hdr_pad);
 void cdf_sync(int ncid);
 
-void cdf_inq (int ncid, int *ndimsp, int *nvarsp, int *ngattsp, int *unlimdimidp);
+void cdf_inq(int ncid, int *ndimsp, int *nvarsp, int *ngattsp, int *unlimdimidp);
 
-void cdf_def_dim (int ncid, const char *name, size_t len, int *idp);
-void cdf_inq_dimid (int ncid, const char *name, int *dimidp);
-void cdf_inq_dim (int ncid, int dimid, char *name, size_t * lengthp);
-void cdf_inq_dimname (int ncid, int dimid, char *name);
-void cdf_inq_dimlen (int ncid, int dimid, size_t * lengthp);
-void cdf_def_var (int ncid, const char *name, nc_type xtype, int ndims, const int dimids[], int *varidp);
+void cdf_def_dim(int ncid, const char *name, size_t len, int *idp);
+void cdf_inq_dimid(int ncid, const char *name, int *dimidp);
+void cdf_inq_dim(int ncid, int dimid, char *name, size_t * lengthp);
+void cdf_inq_dimname(int ncid, int dimid, char *name);
+void cdf_inq_dimlen(int ncid, int dimid, size_t * lengthp);
+void cdf_def_var(int ncid, const char *name, nc_type xtype, int ndims, const int dimids[], int *varidp);
 void cdf_def_var_serial(int ncid, const char *name, nc_type xtype, int ndims, const int dimids[], int *varidp);
 void cdf_inq_varid(int ncid, const char *name, int *varidp);
 void cdf_inq_nvars(int ncid, int *nvarsp);
 void cdf_inq_var(int ncid, int varid, char *name, nc_type *xtypep, int *ndimsp, int dimids[], int *nattsp);
-void cdf_inq_varname (int ncid, int varid, char *name);
-void cdf_inq_vartype (int ncid, int varid, nc_type *xtypep);
-void cdf_inq_varndims (int ncid, int varid, int *ndimsp);
-void cdf_inq_vardimid (int ncid, int varid, int dimids[]);
-void cdf_inq_varnatts (int ncid, int varid, int *nattsp);
+void cdf_inq_varname(int ncid, int varid, char *name);
+void cdf_inq_vartype(int ncid, int varid, nc_type *xtypep);
+void cdf_inq_varndims(int ncid, int varid, int *ndimsp);
+void cdf_inq_vardimid(int ncid, int varid, int dimids[]);
+void cdf_inq_varnatts(int ncid, int varid, int *nattsp);
 
 void cdf_copy_att (int ncid_in, int varid_in, const char *name, int ncid_out, int varid_out);
 void cdf_put_var_text   (int ncid, int varid, const char *tp);
@@ -64,8 +64,10 @@ void cdf_put_vara_double(int ncid, int varid, const size_t start[], const size_t
 
 void cdf_get_vara_float(int ncid, int varid, const size_t start[], const size_t count[], float *fp);
 void cdf_put_vara_float(int ncid, int varid, const size_t start[], const size_t count[], const float *fp);
-void  cdf_get_vara_int(int ncid, int varid, const size_t start[],
-                       const size_t count[], int *dp);
+
+void cdf_get_vara_int(int ncid, int varid, const size_t start[], const size_t count[], int *dp);
+
+void cdf_put_vara(int ncid, int varid, const size_t start[], const size_t count[], const void *cp);
 
 void cdf_put_att_text(int ncid, int varid, const char *name, size_t len, const char *tp);
 void cdf_put_att_int(int ncid, int varid, const char *name, nc_type xtype, size_t len, const int *ip);

src/cdf_write.c

@@ -38,6 +38,50 @@ void cdfDefVarDeflate(int ncid, int ncvarID, int deflate_level)
 #endif
 }
 
+#ifdef HAVE_NETCDF4
+static
+nc_type cdfTypeComplexFloat(stream_t *streamptr)
+{
+  if (streamptr->nc_complex_float_id == CDI_UNDEFID)
+    {
+      typedef struct complex_float { float r, i; } complex_float;
+      const int fileID  = streamptr->fileID;
+      int nc_complex_id;
+      int status;
+      status = nc_def_compound(fileID, sizeof(complex_float), "complex_float", &nc_complex_id);
+      if (status != NC_NOERR) Error("%s", nc_strerror(status));
+      status = nc_insert_compound(fileID, nc_complex_id, "r", NC_COMPOUND_OFFSET(complex_float, r), NC_FLOAT);
+      if (status != NC_NOERR) Error("%s", nc_strerror(status));
+      status = nc_insert_compound(fileID, nc_complex_id, "i", NC_COMPOUND_OFFSET(complex_float, i), NC_FLOAT);
+      if (status != NC_NOERR) Error("%s", nc_strerror(status));
+      streamptr->nc_complex_float_id = nc_complex_id;
+    }
+
+  return (nc_type) streamptr->nc_complex_float_id;
+}
+
+static
+nc_type cdfTypeComplexDouble(stream_t *streamptr)
+{
+  if (streamptr->nc_complex_double_id == CDI_UNDEFID)
+    {
+      typedef struct complex_double { double r, i; } complex_double;
+      const int fileID  = streamptr->fileID;
+      int nc_complex_id;
+      int status;
+      status = nc_def_compound(fileID, sizeof(complex_double), "complex_double", &nc_complex_id);
+      if (status != NC_NOERR) Error("%s", nc_strerror(status));
+      status = nc_insert_compound(fileID, nc_complex_id, "r", NC_COMPOUND_OFFSET(complex_double, r), NC_DOUBLE);
+      if (status != NC_NOERR) Error("%s", nc_strerror(status));
+      status = nc_insert_compound(fileID, nc_complex_id, "i", NC_COMPOUND_OFFSET(complex_double, i), NC_DOUBLE);
+      if (status != NC_NOERR) Error("%s", nc_strerror(status));
+      streamptr->nc_complex_double_id = nc_complex_id;
+    }
+
+  return (nc_type) streamptr->nc_complex_double_id;
+}
+#endif
+
 nc_type cdfDefDatatype(int datatype, stream_t *streamptr)
 {
   nc_type xtype = NC_FLOAT;
@@ -51,8 +95,8 @@ nc_type cdfDefDatatype(int datatype, stream_t *streamptr)
       else if ( datatype == CDI_DATATYPE_UINT8  ) xtype = NC_UBYTE;
       else if ( datatype == CDI_DATATYPE_UINT16 ) xtype = NC_USHORT;
       else if ( datatype == CDI_DATATYPE_UINT32 ) xtype = NC_UINT;
-      else if ( datatype == CDI_DATATYPE_CPX32 || datatype == CDI_DATATYPE_CPX64 )
-        Error("CDI library does not support complex numbers with NetCDF4/HDF5!");
+      else if ( datatype == CDI_DATATYPE_CPX32  ) xtype = cdfTypeComplexFloat(streamptr);
+      else if ( datatype == CDI_DATATYPE_CPX64  ) xtype = cdfTypeComplexDouble(streamptr);
 #else
       else if ( datatype == CDI_DATATYPE_UINT8  ) xtype = NC_SHORT;
       else if ( datatype == CDI_DATATYPE_UINT16 ) xtype = NC_INT;
@@ -93,6 +137,7 @@ void cdfDefVarMissval(stream_t *streamptr, int varID, int dtype, int lcheck)
       if ( lcheck && streamptr->ncmode == 2 ) cdf_redef(fileID);
 
       nc_type xtype = cdfDefDatatype(dtype, streamptr);
+      printf("xtype %d\n", xtype);
       if ( xtype == NC_BYTE && missval > 127 && missval < 256 ) xtype = NC_INT;
 
       if ( lcheck == 0 ||
@@ -1033,10 +1078,43 @@ void cdf_write_var_data(int fileID, int vlistID, int varID, int ncvarID, int dty
         }
     }
 
-  if ( memtype == MEMTYPE_FLOAT )
-    cdf_put_vara_float(fileID, ncvarID, start, count, pdata_sp);
+  if (dtype == CDI_DATATYPE_CPX32 || dtype == CDI_DATATYPE_CPX64)
+    {
+      void *pdata = (memtype == MEMTYPE_FLOAT) ? (void*)pdata_sp : (void*)pdata_dp;
+      float *cdata_sp = NULL;
+      double *cdata_dp = NULL;
+      if (memtype == MEMTYPE_FLOAT && dtype == CDI_DATATYPE_CPX64)
+        {
+          cdata_dp = (double *) Malloc(2*nvals*sizeof(double));
+          for (size_t i = 0; i < nvals; i++)
+            {
+              cdata_dp[2*i] = (double)(pdata_sp[2*i]);
+              cdata_dp[2*i+1] = (double)(pdata_sp[2*i+1]);
+            }
+          pdata = cdata_dp;
+        }
+      else if (memtype == MEMTYPE_DOUBLE && dtype == CDI_DATATYPE_CPX32)
+        {
+          cdata_sp = (float *) Malloc(2*nvals*sizeof(float));
+          for (size_t i = 0; i < nvals; i++)
+            {
+              cdata_sp[2*i] = (float)(pdata_dp[2*i]);
+              cdata_sp[2*i+1] = (float)(pdata_dp[2*i+1]);
+            }
+          pdata = cdata_sp;
+        }
+
+      cdf_put_vara(fileID, ncvarID, start, count, pdata);
+      if (cdata_sp) Free(cdata_sp);
+      if (cdata_dp) Free(cdata_dp);
+    }
   else
-    cdf_put_vara_double(fileID, ncvarID, start, count, pdata_dp);
+    {
+      if ( memtype == MEMTYPE_FLOAT )
+        cdf_put_vara_float(fileID, ncvarID, start, count, pdata_sp);
+      else
+        cdf_put_vara_double(fileID, ncvarID, start, count, pdata_dp);
+    }
 
   if ( mdata_dp ) Free(mdata_dp);
   if ( sdata_dp ) Free(sdata_dp);

src/cdi_int.h

@@ -282,6 +282,8 @@ typedef struct {
   int         ncmode;
   int         vlistID;
 #ifdef HAVE_LIBNETCDF
+  int         nc_complex_float_id;
+  int         nc_complex_double_id;
   ncgrid_t    ncgrid[MAX_GRIDS_PS];
   int         zaxisID[MAX_ZAXES_PS];	// Warning: synchronous array to vlist_to_pointer(vlistID)->zaxisIDs
   int         nczvarID[MAX_ZAXES_PS];

src/stream.c

@@ -890,6 +890,9 @@ void streamDefaultValue ( stream_t * streamptr )
   basetimeInit(&streamptr->basetime);
 
 #ifdef HAVE_LIBNETCDF
+  streamptr->nc_complex_float_id  = CDI_UNDEFID;
+  streamptr->nc_complex_double_id = CDI_UNDEFID;
+
   for ( int i = 0; i < MAX_ZAXES_PS; i++ ) streamptr->zaxisID[i]  = CDI_UNDEFID;
   for ( int i = 0; i < MAX_ZAXES_PS; i++ ) streamptr->nczvarID[i] = CDI_UNDEFID;