/*

Copyright (C) 2003 Motorola Inc
Copyright (C) 2003 David Bateman

This program is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.

This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with this program; see the file COPYING.  If not, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

In addition to the terms of the GPL, you are permitted to link
this program with any Open Source program, as defined by the
Open Source Initiative (www.opensource.org)

*/

#if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
#pragma implementation
#endif

#include <climits>
#include <iostream>
#include <iomanip>

#include <octave/config.h>
#include <octave/oct-obj.h>
#include <octave/ov.h>
#include <octave/lo-ieee.h>
#include <octave/gripes.h>
#include <octave/unwind-prot.h>
#include <octave/cmd-edit.h>
#include <octave/utils.h>
#include <octave/variables.h>

#ifdef CLASS_HAS_LOAD_SAVE
#include <octave/data-conv.h>
#include <octave/byte-swap.h>
#include <octave/ls-oct-ascii.h>
#include <octave/ls-hdf5.h>
#endif

#ifdef NEED_OCTAVE_QUIT
#define OCTAVE_QUIT do {} while (0)
#else
#include <octave/quit.h>
#endif

#include <memory>
#ifndef OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size) \
  std::auto_ptr<T> buf ## _auto_ptr (new T [size]); \
  T *buf = buf ## _auto_ptr.get ()
#endif

#include "fixed-def.h"
#include "ov-base-fixed-mat.h"
#include "ov-base-fixed-mat.cc"
#include "int/fixed.h"
#include "fixed-var.h"
#include "ov-fixed.h"
#include "ov-fixed-mat.h"
#include "ov-fixed-complex.h"
#include "ov-fixed-cx-mat.h"
#include "fixed-conv.h"

#if ! defined (UCHAR_MAX)
#define UCHAR_MAX 255
#endif

template class octave_base_fixed_matrix<FixedComplexMatrix>;

DEFINE_OCTAVE_ALLOCATOR (octave_fixed_complex_matrix);

#ifdef TYPEID_HAS_CLASS
DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_fixed_complex_matrix, 
				     "fixed complex matrix",
				     "FixedPoint");
#else
DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_fixed_complex_matrix, 
				     "fixed complex matrix");
#endif

#ifdef HAVE_ND_ARRAYS
NDArray
octave_fixed_complex_matrix::array_value (bool force_conversion) const
{
  NDArray retval;
  int flag = force_conversion;

#if defined(HAVE_OK_TO_LOSE_IMAGINARY_PART)
  if (! flag)
    flag = Vok_to_lose_imaginary_part;
#else
  if (! flag)
    flag = (Vwarn_imag_to_real ? -1 : 1);
#endif

  if (flag > 0)
    {
      int nr = rows ();
      int nc = columns ();
      dim_vector dv(nr,nc);
      retval.resize (dv);

      for (int i=0; i<nr; i++)
	for (int j=0; j<nc; j++)
	  retval(i + j*nr) = real (matrix(i,j).fixedpoint());
    }
  else
    gripe_implicit_conversion ("fixed complex", "matrix");

  return retval;
}

ComplexNDArray
octave_fixed_complex_matrix::complex_array_value (bool) const
{
  int nr = rows ();
  int nc = columns ();
  dim_vector dv(nr,nc);
  ComplexNDArray retval (dv);

  for (int i=0; i<nr; i++)
    for (int j=0; j<nc; j++)
      retval(i + j*nr) = matrix(i,j).fixedpoint();

  return retval;
}
#endif

#if defined (HAVE_OCTAVE_CONCAT) || defined (HAVE_OLD_OCTAVE_CONCAT)
octave_value
octave_fixed_complex_matrix::resize (const dim_vector& dv) const
{
  if (dv.length() > 2)
    {
      error ("Can not resize fixed point to NDArray");
      return octave_value ();
    }
  FixedComplexMatrix retval (matrix); 
  retval.resize (dv(0), dv(1)); 
  return new octave_fixed_complex_matrix (retval);
}
#endif

FixedComplexMatrix
octave_fixed_complex_matrix::do_index_intern (const octave_value_list& idx,
				     int resize_ok)
{
  FixedComplexMatrix retval;

  int len = idx.length ();

  switch (len)
    {
    case 2:
      {
	idx_vector i = idx (0).index_vector ();
	idx_vector j = idx (1).index_vector ();
	retval = FixedComplexMatrix( matrix.index (i, j, resize_ok));
      }
      break;

    case 1:
      {
	idx_vector i = idx (0).index_vector ();
	retval = FixedComplexMatrix( matrix.index (i, resize_ok));
      }
      break;

    default:
      {
	std::string n = type_name ();

	error ("invalid number of indices (%d) for %s value",
	       len, n.c_str ());
      }
      break;
    }

  return retval;
}

octave_value
octave_fixed_complex_matrix::do_index_op (const octave_value_list& idx,
				     int resize_ok)
{
  octave_value retval;

  FixedComplexMatrix new_matrix = do_index_intern (idx, resize_ok);
  
  if (!error_state) {
    retval = new octave_fixed_complex_matrix ( new_matrix);
    retval.maybe_mutate();
  }

  return retval;
}

octave_value
octave_fixed_complex_matrix::subsasgn (const std::string& type,
				  const std::list<octave_value_list>& idx,
				  const octave_value& rhs)
{
  octave_value retval;

  switch (type[0])
    {
    case '(':
      {
	if (type.length () == 1)
	  retval = numeric_assign (type, idx, rhs);
	else if (type.length () == 2) 
	  {
	    std::list<octave_value_list>::const_iterator p = idx.begin ();
	    octave_value_list key_idx = *++p;

	    assert (key_idx.length () == 1);

	    std::string key = key_idx(0).string_value ();

	    if (key == __FIXED_SIGN_STR)
	      error("can not directly change the sign in a fixed structure");
	    else if (key == __FIXED_VALUE_STR) 
	      error("can not directly change the value of a fixed structure");
	    else if (key == __FIXED_DECSIZE_STR) {
	      if (rhs.is_matrix_type()) {
		FixedComplexMatrix old_matrix = 
		  do_index_intern(idx.front(), 0);
		octave_value new_matrix = new octave_fixed_complex_matrix(
		    old_matrix.chdecsize(rhs.complex_matrix_value()));
		retval = numeric_assign (type, idx, new_matrix);
	      } else {
		FixedComplexMatrix old_matrix = 
		  do_index_intern(idx.front(), 0);
		octave_value new_matrix = new octave_fixed_complex_matrix(
		    old_matrix.chdecsize(rhs.complex_value()));
		retval = numeric_assign (type, idx, new_matrix);
	      }
	    } else if (key == __FIXED_INTSIZE_STR) {
	      if (rhs.is_matrix_type()) {
		FixedComplexMatrix old_matrix = 
		  do_index_intern(idx.front(), 0);
		octave_value new_matrix = new octave_fixed_complex_matrix(
		    old_matrix.chintsize(rhs.complex_matrix_value()));
		retval = numeric_assign (type, idx, new_matrix);
	      } else {
		FixedComplexMatrix old_matrix = 
		  do_index_intern(idx.front(), 0);
		octave_value new_matrix = new octave_fixed_complex_matrix(
		    old_matrix.chintsize(rhs.complex_value()));
		retval = numeric_assign (type, idx, new_matrix);
	      }
	    } else
	      error ("fixed point structure has no member `%s'", 
		     key.c_str ());    
	  }
	else
	  {
	    std::string nm = type_name ();
	    error ("in indexed assignment of %s, illegal assignment", 
		   nm.c_str ());
	  }
      }
      break;

    case '.':
      {
	octave_value_list key_idx = idx.front ();

	assert (key_idx.length () == 1);

	std::string key = key_idx(0).string_value ();

	if (key == __FIXED_SIGN_STR)
	  error("can not directly change the sign in a fixed structure");
	else if (key == __FIXED_VALUE_STR) 
	  error("can not directly change the value of a fixed structure");
	else if (key == __FIXED_DECSIZE_STR) {
	  if (rhs.is_matrix_type())
	    retval = new octave_fixed_complex_matrix (matrix.chdecsize(
				rhs.complex_matrix_value()));
	  else
	    retval = new octave_fixed_complex_matrix (matrix.chdecsize(
				rhs.complex_value()));
	} else if (key == __FIXED_INTSIZE_STR) {
	  if (rhs.is_matrix_type())
	    retval = new octave_fixed_complex_matrix (matrix.chintsize(
				rhs.complex_matrix_value()));
	  else
	    retval = new octave_fixed_complex_matrix (matrix.chintsize(
				rhs.complex_value()));
	} else
	  error ("fixed point structure has no member `%s'", key.c_str ());    
      }
      break;

    case '{':
      {
	std::string nm = type_name ();
	error ("%s cannot be indexed with %c", nm.c_str (), type[0]);
      }
      break;

    default:
      panic_impossible ();
    }

  return retval;
}

octave_value *
octave_fixed_complex_matrix::try_narrowing_conversion (void)
{
  octave_value *retval = 0;

  int nr = matrix.rows ();
  int nc = matrix.cols ();

  // Why doesn't this work!!! Got to get it to work !!!
  if (nr == 1 && nc == 1)
    {
      FixedPointComplex c = matrix (0, 0);

      if (imag (c) == 0.0)
	retval = new octave_fixed (real (c));
      else
	retval = new octave_fixed_complex (c);
    }
  else if (nr == 0 || nc == 0)
    retval = new octave_fixed_matrix (FixedMatrix (nr, nc));
  else if (matrix.all_elements_are_real ())
    retval = new octave_fixed_matrix (real (matrix));

  return retval;
}

double
octave_fixed_complex_matrix::double_value (bool force_conversion) const
{
  double retval = lo_ieee_nan_value ();

  int flag = force_conversion;

#if defined(HAVE_OK_TO_LOSE_IMAGINARY_PART)
  if (! flag)
    flag = Vok_to_lose_imaginary_part;
#else
  if (! flag)
    flag = (Vwarn_imag_to_real ? -1 : 1);
#endif

  if (flag < 0)
    gripe_implicit_conversion ("fixed complex matrix", "real scalar");

  if (flag)
    {

#if defined(HAVE_DO_FORTRAN_INDEXING)
      if ((rows () == 1 && columns () == 1)
	  || (Vdo_fortran_indexing && rows () > 0 && columns () > 0))
	retval = std::real (matrix (0, 0) .fixedpoint());
#else
      if (rows () > 0 && columns () > 0)
	{
	  if (Vwarn_fortran_indexing)
	    gripe_implicit_conversion ("real matrix", "real scalar");
	  
	  retval = std::real (matrix (0, 0) .fixedpoint());
	}
#endif
      else
	gripe_invalid_conversion ("fixed complex matrix", "real scalar");
    }
  else
    gripe_invalid_conversion ("fixed complex matrix", "real scalar");

  return retval;
}

FixedPoint
octave_fixed_complex_matrix::fixed_value (bool force_conversion) const
{
  FixedPoint retval;

  int flag = force_conversion;

#if defined(HAVE_OK_TO_LOSE_IMAGINARY_PART)
  if (! flag)
    flag = Vok_to_lose_imaginary_part;
#else
  if (! flag)
    flag = (Vwarn_imag_to_real ? -1 : 1);
#endif

  if (flag < 0)
    gripe_implicit_conversion ("fixed complex matrix", "fixed scalar");

  if (flag)
    {
#if defined(HAVE_DO_FORTRAN_INDEXING)
      if ((rows () == 1 && columns () == 1)
	  || (Vdo_fortran_indexing && rows () > 0 && columns () > 0))
	retval = real( matrix (0, 0));
#else
      if (rows () > 0 && columns () > 0)
	{
	  if (Vwarn_fortran_indexing)
	    gripe_implicit_conversion ("real matrix", "real scalar");
	  
	  retval = real( matrix (0, 0));
	}
#endif
      else
	gripe_invalid_conversion ("fixed complex matrix", "fixed scalar");
    }
  else
    gripe_invalid_conversion ("fixed complex matrix", "fixed scalar");

  return retval;
}

Matrix
octave_fixed_complex_matrix::matrix_value (bool force_conversion) const
{
  Matrix retval;

  int flag = force_conversion;

#if defined(HAVE_OK_TO_LOSE_IMAGINARY_PART)
  if (! flag)
    flag = Vok_to_lose_imaginary_part;
#else
  if (! flag)
    flag = (Vwarn_imag_to_real ? -1 : 1);
#endif

  if (flag < 0)
    gripe_implicit_conversion ("fixed complex matrix", "real matrix");

  if (flag)
    retval = ::real (matrix.fixedpoint());
  else
    gripe_invalid_conversion ("fixed complex matrix", "real matrix");

  return retval;
}

FixedMatrix
octave_fixed_complex_matrix::fixed_matrix_value (bool force_conversion) const
{
  FixedMatrix retval;

  int flag = force_conversion;

#if defined(HAVE_OK_TO_LOSE_IMAGINARY_PART)
  if (! flag)
    flag = Vok_to_lose_imaginary_part;
#else
  if (! flag)
    flag = (Vwarn_imag_to_real ? -1 : 1);
#endif

  if (flag < 0)
    gripe_implicit_conversion ("fixed complex matrix", "fixed matrix");

  if (flag)
    retval = real (matrix);
  else
    gripe_invalid_conversion ("fixed complex matrix", "fixed matrix");

  return retval;
}

Complex
octave_fixed_complex_matrix::complex_value (bool) const
{
  double tmp = lo_ieee_nan_value ();

  Complex retval (tmp, tmp);

#if defined(HAVE_DO_FORTRAN_INDEXING)
      if ((rows () == 1 && columns () == 1)
	  || (Vdo_fortran_indexing && rows () > 0 && columns () > 0))
	retval = matrix (0, 0) .fixedpoint();
#else
      if (rows () > 0 && columns () > 0)
	{
	  if (Vwarn_fortran_indexing)
	    gripe_implicit_conversion ("real matrix", "real scalar");
	  
	  retval = matrix (0, 0) .fixedpoint();
	}
#endif
  else
    gripe_invalid_conversion ("fixed matrix", "complex scalar");

  return retval;
}

void
octave_fixed_complex_matrix::print_raw (std::ostream& os,
				   bool pr_as_read_syntax) const
{
  int prec = check_preference("output_precision");
  double min_num = std::max(abs(real(matrix)).row_min().min().fixedpoint(),
			    abs(imag(matrix)).row_min().min().fixedpoint());
  int new_prec = (int)std::max(real(matrix).getdecsize().row_max().max(),
			       imag(matrix).getdecsize().row_max().max()) +
    (min_num >= 1. ? (int)log10(min_num) + 1 : 0);

  bind_builtin_variable ("output_precision", new_prec); 

  octave_print_internal (os, complex_matrix_value(), false, 
			 current_print_indent_level ());

  bind_builtin_variable ("output_precision", prec);
}

#ifdef CLASS_HAS_LOAD_SAVE
bool 
octave_fixed_complex_matrix::save_ascii (std::ostream& os, 
					 bool& infnan_warned, 
					 bool strip_nan_and_inf)
{
  dim_vector d = dims ();
  os << "# ndims: " << d.length () << "\n";

  for (int i=0; i < d.length (); i++)
    os << " " << d (i);

  FixedMatrix re (real (matrix));
  FixedMatrix im (imag (matrix));
  os << "\n" << re.getintsize () << im.getintsize ()
     << re.getdecsize () << im.getdecsize () 
     << re.fixedpoint() << im.fixedpoint ();

  return true;
}

bool 
octave_fixed_complex_matrix::load_ascii (std::istream& is)
{
  int mdims;
  bool success = true;

  if (extract_keyword (is, "ndims", mdims))
    {
      dim_vector dv;
      dv.resize (mdims);
      
      for (int i = 0; i < mdims; i++)
	is >> dv(i);

      if (dv.length() != 2)
	{
	  error ("load: N-D fixed arrays not supported");
	  success = false;
	}
      else
	{
	  Matrix rintsize (dv(0), dv(1)), rdecsize (dv(0), dv(1)), 
	    rnumber (dv(0), dv(1));
	  Matrix iintsize (dv(0), dv(1)), idecsize (dv(0), dv(1)), 
	    inumber (dv(0), dv(1));

	  is >> rintsize >> iintsize >> rdecsize >> idecsize 
	     >> rnumber >> inumber;

	  if (!is) 
	    {
	      error ("load: failed to load matrix constant");
	      success = false;
	    }
	  matrix = 
	    FixedComplexMatrix (FixedMatrix (rintsize, rdecsize, rnumber),
				FixedMatrix (iintsize, idecsize, inumber));

	}
    }
  else 
    {
      error ("load: failed to extract dimension of fixed point variable");
      success = false;
    }

  return success;;
}

bool 
octave_fixed_complex_matrix::save_binary (std::ostream& os, 
					  bool& save_as_floats)
{
  dim_vector d = dims ();

  // Only treat 2-D array for now
  if (d.length() != 2)
    return false;

  // Use negative value for ndims to be consistent with other types
  FOUR_BYTE_INT tmp = - d.length();
  os.write (X_CAST (char *, &tmp), 4);
  for (int i=0; i < d.length (); i++)
    {
      tmp = d(i);
      os.write (X_CAST (char *, &tmp), 4);
    }

  char size = (char) sizeof (unsigned int);
  os.write (X_CAST (char *, &size), 1);

  // intsize and decsize are integers in the range [0:32], so store as char
  FixedMatrix re (real (matrix)), im (imag (matrix));
  LS_DO_WRITE (char, re.getintsize ().fortran_vec (), 1, d.numel (), os);
  LS_DO_WRITE (char, im.getintsize ().fortran_vec (), 1, d.numel (), os);
  LS_DO_WRITE (char, re.getdecsize ().fortran_vec (), 1, d.numel (), os);
  LS_DO_WRITE (char, im.getdecsize ().fortran_vec (), 1, d.numel (), os);
  LS_DO_WRITE (unsigned int, re.getnumber ().fortran_vec (), 
	       sizeof (unsigned int), d.numel (), os);
  LS_DO_WRITE (unsigned int, im.getnumber ().fortran_vec (), 
	       sizeof (unsigned int), d.numel (), os);

  return true;
}

#ifdef HAVE_SWAP_BYTES
static inline void swap_4_bytes (volatile void *ptr)
{
  swap_bytes <4> (ptr);
}
#endif

bool 
octave_fixed_complex_matrix::load_binary (std::istream& is, bool swap,
				 oct_mach_info::float_format fmt)
{
  FOUR_BYTE_INT mdims;
  if (! is.read (X_CAST (char *, &mdims), 4))
    return false;
  if (swap)
    swap_4_bytes (X_CAST (char *, &mdims));

  if (mdims != -2)
    return false;

  mdims = - mdims;
  FOUR_BYTE_INT di;
  dim_vector dv;
  dv.resize (mdims);

  for (int i = 0; i < mdims; i++)
    {
      if (! is.read (X_CAST (char *, &di), 4))
	return false;
      if (swap)
	swap_4_bytes (X_CAST (char *, &di));
      dv(i) = di;
    }

  char size;
  Matrix rintsize (dv(0), dv(1)), rdecsize (dv(0), dv(1)), 
    rnumber (dv(0), dv(1));
  Matrix iintsize (dv(0), dv(1)), idecsize (dv(0), dv(1)), 
    inumber (dv(0), dv(1));

  if (! is.read (X_CAST (char *, &size), 1))
    return false;

  LS_DO_READ_1(rintsize.fortran_vec (), dv.numel (), is);
  LS_DO_READ_1(iintsize.fortran_vec (), dv.numel (), is);
  LS_DO_READ_1(rdecsize.fortran_vec (), dv.numel (), is);
  LS_DO_READ_1(idecsize.fortran_vec (), dv.numel (), is);

  if (size == 4)
    {
      LS_DO_READ(unsigned int, swap, rnumber.fortran_vec (), 4, 
		 dv.numel (), is);
      LS_DO_READ(unsigned int, swap, inumber.fortran_vec (), 4, 
		 dv.numel (), is);
    }
  else if (size == 8)
    {
      LS_DO_READ(unsigned int, swap, rnumber.fortran_vec (), 8, 
		 dv.numel (), is);
      LS_DO_READ(unsigned int, swap, inumber.fortran_vec (), 8, 
		 dv.numel (), is);
    }
  else
    return false;

  if (error_state || ! is)
    return false;
  
  // This is ugly, is there a better way?
  matrix.resize (dv(0), dv(1));
  for (int i = 0; i < dv(0); i++)
    for (int j = 0; j < dv(1); j++)
      matrix (i, j) = 
	FixedPointComplex (FixedPoint ((unsigned int)rintsize (i, j), 
				       (unsigned int)rdecsize (i, j), 
				       (unsigned int)rnumber (i, j)),
			   FixedPoint ((unsigned int)iintsize (i, j), 
				       (unsigned int)idecsize (i, j), 
				       (unsigned int)inumber (i, j)));
  
  return true;
}

#if defined (HAVE_HDF5)
bool
octave_fixed_complex_matrix::save_hdf5 (hid_t loc_id, const char *name, 
					bool save_as_floats)
{
  hid_t group_hid = -1;
  group_hid = H5Gcreate (loc_id, name, 0);
  if (group_hid < 0 ) return false;

  dim_vector d = dims ();
  hsize_t hdims[d.length () > 2 ? d.length () : 3];
  hid_t space_hid = -1, data_hid = -1, type_hid = -1;
  int rank = ( (d (0) == 1) && (d.length () == 2) ? 1 : d.length ());
  bool retval = true;

  // Octave uses column-major, while HDF5 uses row-major ordering
  for (int i = 0, j = d.length() - 1; i < d.length (); i++, j--)
    hdims[i] = d (j);

  space_hid = H5Screate_simple (rank, hdims, (hsize_t*) 0);
  if (space_hid < 0) 
    {
      H5Gclose (group_hid);
      return false;
    }

  type_hid = hdf5_make_fixed_complex_type (H5T_NATIVE_UCHAR, 1);
  if (type_hid < 0) 
    {
      H5Sclose (space_hid);
      H5Gclose (group_hid);
      return false;
    }

  data_hid = H5Dcreate (group_hid, "int", type_hid, space_hid, H5P_DEFAULT);
  if (data_hid < 0) 
    {
      H5Sclose (space_hid);
      H5Tclose (type_hid);
      H5Gclose (group_hid);
      return false;
    }
  
  Complex * m = matrix.getintsize ().fortran_vec ();
  OCTAVE_LOCAL_BUFFER(unsigned char, tmp, 2 * d.numel ());
  for (int i = 0; i < d.numel (); i++)
    {
      tmp[i<<1] = (unsigned char) real (m[i]);
      tmp[(i<<1)+1] = (unsigned char) imag (m[i]);
    }
  retval = H5Dwrite (data_hid, type_hid, H5S_ALL, H5S_ALL, H5P_DEFAULT,
		     (void*) tmp) >= 0;
  H5Dclose (data_hid);
  if (!retval)
    {
      H5Sclose (space_hid);
      H5Tclose (type_hid);
      H5Gclose (group_hid);
      return false;
    }    

  data_hid = H5Dcreate (group_hid, "dec", type_hid, space_hid, H5P_DEFAULT);
  if (data_hid < 0) 
    {
      H5Sclose (space_hid);
      H5Tclose (type_hid);
      H5Gclose (group_hid);
      return false;
    }
  
  m = matrix.getdecsize (). fortran_vec ();
  for (int i = 0; i < d.numel (); i++)
    {
      tmp[i<<1] = (unsigned char) real (m[i]);
      tmp[(i<<1)+1] = (unsigned char) imag (m[i]);
    }
  retval = H5Dwrite (data_hid, type_hid, H5S_ALL, H5S_ALL, H5P_DEFAULT,
		     (void*) tmp) >= 0;
  H5Dclose (data_hid);
  H5Tclose (type_hid);

  if (!retval)
    {
      H5Sclose (space_hid);
      H5Gclose (group_hid);
      return false;
    }    

  type_hid = hdf5_make_fixed_complex_type (H5T_NATIVE_UINT, 
					   sizeof (unsigned int));
  if (type_hid < 0) 
    {
      H5Sclose (space_hid);
      H5Gclose (group_hid);
      return false;
    }

  data_hid = H5Dcreate (group_hid, "num", type_hid, space_hid, H5P_DEFAULT);
  if (data_hid < 0) 
    {
      H5Sclose (space_hid);
      H5Tclose (type_hid);
      H5Gclose (group_hid);
      return false;
    }
  
  m = matrix.getnumber ().fortran_vec ();
  OCTAVE_LOCAL_BUFFER(unsigned int, num, 2*d.numel ());
  for (int i = 0; i < d.numel (); i++)
    {
      num[i<<1] = (unsigned int) real (m[i]);
      num[(i<<1)+1] = (unsigned int) imag (m[i]);
    }
  retval = H5Dwrite (data_hid, type_hid, H5S_ALL, H5S_ALL, H5P_DEFAULT,
		     (void*) num) >= 0;
  H5Dclose (data_hid);
  H5Sclose (space_hid);
  H5Tclose (type_hid);
  H5Gclose (group_hid);
  return retval;
}

bool
octave_fixed_complex_matrix::load_hdf5 (hid_t loc_id, const char *name,
					bool have_h5giterate_bug)
{
  herr_t retval = -1;
  hid_t group_hid, data_hid, space_id, type_hid;
  hsize_t rank, rank_old;

  group_hid = H5Gopen (loc_id, name);
  if (group_hid < 0 ) return false;

  hid_t complex_type = hdf5_make_fixed_complex_type (H5T_NATIVE_UINT, 
						     sizeof(unsigned int));

  data_hid = H5Dopen (group_hid, "int");
  type_hid = H5Dget_type (data_hid);

  if (! hdf5_types_compatible (type_hid, complex_type))
    {
      H5Tclose(complex_type);
      H5Tclose(type_hid);
      H5Dclose (data_hid);
      H5Gclose (group_hid);
      return false;
    }

  space_id = H5Dget_space (data_hid);
  rank = H5Sget_simple_extent_ndims (space_id);
  rank_old = rank;

  if (rank < 1 || rank > 2)
    {
      // No N-D array yet
      H5Tclose(complex_type);
      H5Tclose(type_hid);
      H5Sclose (space_id);
      H5Dclose (data_hid);
      H5Gclose (group_hid);
      return false;
    }

  OCTAVE_LOCAL_BUFFER (hsize_t, hdims, rank);
  OCTAVE_LOCAL_BUFFER (hsize_t, maxdims, rank);

  H5Sget_simple_extent_dims (space_id, hdims, maxdims);

  dim_vector dv;
  dim_vector dv_old;

  // Octave uses column-major, while HDF5 uses row-major ordering
  if (rank == 1)
    {
      dv.resize (2);
      dv(0) = 1;
      dv(1) = hdims[0];
    }
  else
    {
      dv.resize (rank);
      for (int i = 0, j = rank - 1; i < (int)rank; i++, j--)
	dv(j) = hdims[i];
    }
  dv_old = dv;

  OCTAVE_LOCAL_BUFFER (unsigned int, intsize, 2 * dv.numel ());
  if (H5Dread (data_hid, complex_type, H5S_ALL, H5S_ALL, H5P_DEFAULT, 
	       (void *) intsize) < 0) 
    {
      H5Tclose(complex_type);
      H5Tclose(type_hid);
      H5Sclose (space_id);
      H5Dclose (data_hid);
      H5Gclose (group_hid);
      return false;
    }
  H5Tclose(type_hid);
  H5Dclose (data_hid);

  data_hid = H5Dopen (group_hid, "dec");
  type_hid = H5Dget_type (data_hid);

  if (! hdf5_types_compatible (type_hid, complex_type))
    {
      H5Tclose(complex_type);
      H5Tclose(type_hid);
      H5Dclose (data_hid);
      H5Gclose (group_hid);
      return false;
    }

  space_id = H5Dget_space (data_hid);
  rank = H5Sget_simple_extent_ndims (space_id);

  if (rank != rank_old)
    {
      H5Tclose(complex_type);
      H5Tclose(type_hid);
      H5Sclose (space_id);
      H5Dclose (data_hid);
      H5Gclose (group_hid);
      return false;
    }

  H5Sget_simple_extent_dims (space_id, hdims, maxdims);

  // Octave uses column-major, while HDF5 uses row-major ordering
  if (rank == 1)
    {
      dv.resize (2);
      dv(0) = 1;
      dv(1) = hdims[0];
    }
  else
    {
      dv.resize (rank);
      for (int i = 0, j = rank - 1; i < (int)rank; i++, j--)
	dv(j) = hdims[i];
    }

  if (dv_old != dv)
    {
      H5Tclose(complex_type);
      H5Tclose(type_hid);
      H5Sclose (space_id);
      H5Dclose (data_hid);
      H5Gclose (group_hid);
      return false;
    }

  OCTAVE_LOCAL_BUFFER (unsigned int, decsize, 2 * dv.numel ());
  if (H5Dread (data_hid, complex_type, H5S_ALL, H5S_ALL, H5P_DEFAULT, 
	       (void *) decsize) < 0) 
    {
      H5Tclose(complex_type);
      H5Tclose(type_hid);
      H5Sclose (space_id);
      H5Dclose (data_hid);
      H5Gclose (group_hid);
      return false;
    }
  H5Tclose(type_hid);
  H5Dclose (data_hid);

  data_hid = H5Dopen (group_hid, "num");
  type_hid = H5Dget_type (data_hid);

  if (! hdf5_types_compatible (type_hid, complex_type))
    {
      H5Tclose(complex_type);
      H5Tclose(type_hid);
      H5Dclose (data_hid);
      H5Gclose (group_hid);
      return false;
    }

  space_id = H5Dget_space (data_hid);
  rank = H5Sget_simple_extent_ndims (space_id);

  if (rank != rank_old)
    {
      H5Tclose(complex_type);
      H5Tclose(type_hid);
      H5Sclose (space_id);
      H5Dclose (data_hid);
      H5Gclose (group_hid);
      return false;
    }

  H5Sget_simple_extent_dims (space_id, hdims, maxdims);

  // Octave uses column-major, while HDF5 uses row-major ordering
  if (rank == 1)
    {
      dv.resize (2);
      dv(0) = 1;
      dv(1) = hdims[0];
    }
  else
    {
      dv.resize (rank);
      for (int i = 0, j = rank - 1; i < (int)rank; i++, j--)
	dv(j) = hdims[i];
    }

  if (dv_old != dv)
    {
      H5Tclose(complex_type);
      H5Tclose(type_hid);
      H5Sclose (space_id);
      H5Dclose (data_hid);
      H5Gclose (group_hid);
      return false;
    }

  OCTAVE_LOCAL_BUFFER (unsigned int, number, 2 * dv.numel ());
  retval = H5Dread (data_hid, complex_type, H5S_ALL, H5S_ALL, 
		    H5P_DEFAULT, (void *) number); 
  H5Tclose(complex_type);
  H5Tclose(type_hid);
  H5Dclose (data_hid);
  H5Sclose (space_id);
  H5Gclose (group_hid);
  if (retval < 0)
    return false;

  // This is ugly, is there a better way?
  matrix.resize (dv(0), dv(1));
  unsigned int * ivec = intsize;
  unsigned int * dvec = decsize;
  unsigned int * nvec = number;
  for (int j = 0; j < dv(1); j++)
    for (int i = 0; i < dv(0); i++)
      matrix (i, j) = 
	FixedPointComplex (FixedPoint (*ivec++, *dvec++, *nvec++),
			   FixedPoint (*ivec++, *dvec++, *nvec++));


  return true;
}
#endif
#endif

/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/