class NMatrix

# File lib/nmatrix/shortcuts.rb, line 100
def [](*params)
  options = params.last.is_a?(Hash) ? params.pop : {}

  # First find the dimensions of the array.
  i = 0
  shape = []
  row = params
  while row.is_a?(Array)
    shape[i] = row.length
    row = row[0]
    i += 1
  end

  # A row vector should be stored as 1xN, not N
  #shape.unshift(1) if shape.size == 1

  # Then flatten the array.
  NMatrix.new(shape, params.flatten, options)
end

# File lib/nmatrix/shortcuts.rb, line 313
def block_diagonal(*params)
  options = params.last.is_a?(Hash) ? params.pop : {}

  params.each_index do |i|
    params[i] = params[i].to_nm if params[i].is_a?(Array) # Convert Array to NMatrix
    params[i] = NMatrix.new([1,1], [params[i]]) if params[i].is_a?(Numeric) # Convert number to NMatrix
  end

  block_sizes = [] #holds the size of each matrix block
  params.each do |b|
    unless b.is_a?(NMatrix)
      raise(ArgumentError, "Only NMatrix or appropriate Array objects or single numbers allowed")
    end
    raise(ArgumentError, "Only 2D matrices or 2D arrays allowed") unless b.shape.size == 2
    raise(ArgumentError, "Only square-shaped blocks allowed") unless b.shape[0] == b.shape[1]
    block_sizes << b.shape[0]
  end

  block_diag_mat = NMatrix.zeros(block_sizes.sum, options)
  (0...params.length).each do |n|
    # First determine the size and position of the n'th block in the block-diagonal matrix
    block_size = block_sizes[n]
    block_pos = block_sizes[0...n].sum
    # populate the n'th block in the block-diagonal matrix
    (0...block_size).each do |i|
      (0...block_size).each do |j|
        block_diag_mat[block_pos+i,block_pos+j] = params[n][i,j]
      end
    end
  end

  return block_diag_mat
end

# File lib/nmatrix/shortcuts.rb, line 265
def diagonal(entries, opts={})
  m = NMatrix.zeros(entries.size,
                    {:dtype => guess_dtype(entries[0]), :capacity => entries.size + 1}.merge(opts)
                   )
  entries.each_with_index do |n, i|
    m[i,i] = n
  end
  m
end

# File lib/nmatrix/shortcuts.rb, line 229
def eye(shape, opts={})
  # Fill the diagonal with 1's.
  m = NMatrix.zeros(shape, {:dtype => :float64}.merge(opts))
  (0...m.shape[0]).each do |i|
    m[i, i] = 1
  end

  m
end

static VALUE nm_guess_dtype(VALUE self, VALUE v) {
  return ID2SYM(rb_intern(DTYPE_NAMES[nm_dtype_guess(v)]));
}

# File lib/nmatrix/nmatrix.rb, line 88
def load_matlab_file(file_path)
  NMatrix::IO::Mat5Reader.new(File.open(file_path, 'rb')).to_ruby
end

# File lib/nmatrix/nmatrix.rb, line 99
def load_pcd_file(file_path)
  NMatrix::IO::PointCloudReader::MetaReader.new(file_path)
end

# File lib/nmatrix/nmatrix.rb, line 136
def meshgrid(vectors, options = {})
  raise(ArgumentError, 'Expected at least 2 arrays.') if vectors.size < 2
  options[:indexing] ||= :xy
  raise(ArgumentError, 'Indexing must be :xy of :ij') unless [:ij, :xy].include? options[:indexing]
  mats = vectors.map { |arr| arr.respond_to?(:flatten) ? arr.flatten : arr.to_flat_array }
  mats[0], mats[1] = mats[1], mats[0] if options[:indexing] == :xy
  new_dim = mats.size
  lengths = mats.map(&:size)
  result = mats.map.with_index do |matrix, axis|
    if options[:sparse]
      new_shape = Array.new(new_dim, 1)
      new_shape[axis] = lengths[axis]
      new_elements = matrix
    else
      before_axis = lengths[0...axis].reduce(:*)
      after_axis = lengths[(axis+1)..-1].reduce(:*)
      new_shape = lengths
      new_elements = after_axis ? matrix.map{ |el| [el] * after_axis }.flatten : matrix
      new_elements *= before_axis if before_axis
    end
    NMatrix.new(new_shape, new_elements)
  end
  result[0], result[1] = result[1], result[0] if options[:indexing] == :xy
  result
end

static VALUE nm_min_dtype(VALUE self, VALUE v) {
  return ID2SYM(rb_intern(DTYPE_NAMES[nm_dtype_min(v)]));
}

static VALUE nm_init(int argc, VALUE* argv, VALUE nm) {
  NM_CONSERVATIVE(nm_register_value(&nm));
  NM_CONSERVATIVE(nm_register_values(argv, argc));

  if (argc <= 3) { // Call the new constructor unless all four arguments are given (or the 7-arg version is given)
    NM_CONSERVATIVE(nm_unregister_values(argv, argc));
    NM_CONSERVATIVE(nm_unregister_value(&nm));
        return nm_init_new_version(argc, argv, nm);
  }

  /* First, determine stype (dense by default) */
  nm::stype_t stype;
  size_t  offset = 0;

  if (!SYMBOL_P(argv[0]) && TYPE(argv[0]) != T_STRING) {
    stype = nm::DENSE_STORE;

  } else {
    // 0: String or Symbol
    stype  = interpret_stype(argv[0]);
    offset = 1;
  }

  // If there are 7 arguments and Yale, refer to a different init function with fewer sanity checks.
  if (argc == 7) {
    if (stype == nm::YALE_STORE) {
      NM_CONSERVATIVE(nm_unregister_values(argv, argc));
      NM_CONSERVATIVE(nm_unregister_value(&nm));
      return nm_init_yale_from_old_yale(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], nm);

    } else {
      NM_CONSERVATIVE(nm_unregister_values(argv, argc));
      NM_CONSERVATIVE(nm_unregister_value(&nm));
      rb_raise(rb_eArgError, "Expected 2-4 arguments (or 7 for internal Yale creation)");
    }
  }

  // 1: Array or Fixnum
  size_t dim;
  size_t* shape = interpret_shape(argv[offset], &dim);

  // 2-3: dtype
  nm::dtype_t dtype = interpret_dtype(argc-1-offset, argv+offset+1, stype);

  size_t init_cap = 0, init_val_len = 0;
  void* init_val  = NULL;
  if (!SYMBOL_P(argv[1+offset]) || TYPE(argv[1+offset]) == T_ARRAY) {
        // Initial value provided (could also be initial capacity, if yale).

    if (stype == nm::YALE_STORE && NM_RUBYVAL_IS_NUMERIC(argv[1+offset])) {
      init_cap = FIX2UINT(argv[1+offset]);

    } else {
        // 4: initial value / dtype
      init_val = interpret_initial_value(argv[1+offset], dtype);

      if (TYPE(argv[1+offset]) == T_ARRAY)      init_val_len = RARRAY_LEN(argv[1+offset]);
      else                                  init_val_len = 1;
    }

  } else {
        // DType is RUBYOBJ.

    if (stype == nm::DENSE_STORE) {
        /*
         * No need to initialize dense with any kind of default value unless it's
         * an RUBYOBJ matrix.
         */
      if (dtype == nm::RUBYOBJ) {
        // Pretend [nil] was passed for RUBYOBJ.
        init_val = NM_ALLOC(VALUE);
        *(VALUE*)init_val = Qnil;

        init_val_len = 1;

      } else {
        init_val = NULL;
      }
    } else if (stype == nm::LIST_STORE) {
      init_val = NM_ALLOC_N(char, DTYPE_SIZES[dtype]);
      std::memset(init_val, 0, DTYPE_SIZES[dtype]);
    }
  }

  if (dtype == nm::RUBYOBJ) {
    nm_register_values(reinterpret_cast<VALUE*>(init_val), init_val_len);
  }

  // TODO: Update to allow an array as the initial value.
  NMATRIX* nmatrix;
  UnwrapNMatrix(nm, nmatrix);

  nmatrix->stype = stype;

  switch (stype) {
    case nm::DENSE_STORE:
      nmatrix->storage = (STORAGE*)nm_dense_storage_create(dtype, shape, dim, init_val, init_val_len);
      break;

    case nm::LIST_STORE:
      nmatrix->storage = (STORAGE*)nm_list_storage_create(dtype, shape, dim, init_val);
      break;

    case nm::YALE_STORE:
      nmatrix->storage = (STORAGE*)nm_yale_storage_create(dtype, shape, dim, init_cap);
      nm_yale_storage_init((YALE_STORAGE*)(nmatrix->storage), NULL);
      break;
  }

  if (dtype == nm::RUBYOBJ) {
    nm_unregister_values(reinterpret_cast<VALUE*>(init_val), init_val_len);
  }

  NM_CONSERVATIVE(nm_unregister_values(argv, argc));
  NM_CONSERVATIVE(nm_unregister_value(&nm));

  return nm;
}

# File lib/nmatrix/shortcuts.rb, line 171
def ones(shape, opts={})
  NMatrix.new(shape, 1, {:dtype => :float64, :default => 1}.merge(opts))
end

# File lib/nmatrix/shortcuts.rb, line 184
def ones_like(nm)
  NMatrix.ones(nm.shape, dtype: nm.dtype, stype: nm.stype, capacity: nm.capacity, default: 1)
end

# File lib/nmatrix/shortcuts.rb, line 370
def random(shape, opts={})
  scale = opts.delete(:scale) || 1.0

  rng = Random.new

  random_values = []


  # Construct the values of the final matrix based on the dimension.
  if opts[:dtype] == :complex64 || opts[:dtype] == :complex128
    NMatrix.size(shape).times { |i| random_values << Complex(rng.rand(scale), rng.rand(scale)) }
  else
    NMatrix.size(shape).times { |i| random_values << rng.rand(scale) }
  end

  NMatrix.new(shape, random_values, {:dtype => :float64, :stype => :dense}.merge(opts))
end

static VALUE nm_read(int argc, VALUE* argv, VALUE self) {
  using std::ifstream;

  NM_CONSERVATIVE(nm_register_values(argv, argc));
  NM_CONSERVATIVE(nm_register_value(&self));

  VALUE file, force_;

  // Read the arguments
  rb_scan_args(argc, argv, "11", &file, &force_);
  bool force   = (force_ != Qnil && force_ != Qfalse);


  if (!RB_FILE_EXISTS(file)) { // FIXME: Errno::ENOENT
    NM_CONSERVATIVE(nm_unregister_values(argv, argc));
    NM_CONSERVATIVE(nm_unregister_value(&self));
    rb_raise(rb_get_errno_exc("ENOENT"), "%s", RSTRING_PTR(file));
  }

  // Open a file stream
  ifstream f(RSTRING_PTR(file), std::ios::in | std::ios::binary);

  uint16_t major, minor, release;
  get_version_info(major, minor, release); // compare to NMatrix version

  uint16_t fmajor, fminor, frelease, null16;

  // READ FIRST 64-BIT BLOCK
  f.read(reinterpret_cast<char*>(&fmajor),   sizeof(uint16_t));
  f.read(reinterpret_cast<char*>(&fminor),   sizeof(uint16_t));
  f.read(reinterpret_cast<char*>(&frelease), sizeof(uint16_t));
  f.read(reinterpret_cast<char*>(&null16),   sizeof(uint16_t));

  int ver  = major * 10000 + minor * 100 + release,
      fver = fmajor * 10000 + fminor * 100 + release;
  if (fver > ver && force == false) {
    NM_CONSERVATIVE(nm_unregister_values(argv, argc));
    NM_CONSERVATIVE(nm_unregister_value(&self));
    rb_raise(rb_eIOError, "File was created in newer version of NMatrix than current (%u.%u.%u)", fmajor, fminor, frelease);
  }
  if (null16 != 0) rb_warn("nm_read: Expected zero padding was not zero (0)\n");

  uint8_t dt, st, it, sm;
  uint16_t dim;

  // READ SECOND 64-BIT BLOCK
  f.read(reinterpret_cast<char*>(&dt), sizeof(uint8_t));
  f.read(reinterpret_cast<char*>(&st), sizeof(uint8_t));
  f.read(reinterpret_cast<char*>(&it), sizeof(uint8_t)); // FIXME: should tell how few bytes indices are stored as
  f.read(reinterpret_cast<char*>(&sm), sizeof(uint8_t));
  f.read(reinterpret_cast<char*>(&null16), sizeof(uint16_t));
  f.read(reinterpret_cast<char*>(&dim), sizeof(uint16_t));

  if (null16 != 0) rb_warn("nm_read: Expected zero padding was not zero (1)");
  nm::stype_t stype = static_cast<nm::stype_t>(st);
  nm::dtype_t dtype = static_cast<nm::dtype_t>(dt);
  nm::symm_t  symm  = static_cast<nm::symm_t>(sm);
  //nm::itype_t itype = static_cast<nm::itype_t>(it);

  // READ NEXT FEW 64-BIT BLOCKS
  size_t* shape = NM_ALLOC_N(size_t, dim);
  read_padded_shape(f, dim, shape);

  STORAGE* s;
  if (stype == nm::DENSE_STORE) {
    s = nm_dense_storage_create(dtype, shape, dim, NULL, 0);
    nm_register_storage(stype, s);

    read_padded_dense_elements(f, reinterpret_cast<DENSE_STORAGE*>(s), symm, dtype);

  } else if (stype == nm::YALE_STORE) {
    uint32_t ndnz, length;

    // READ YALE-SPECIFIC 64-BIT BLOCK
    f.read(reinterpret_cast<char*>(&ndnz),     sizeof(uint32_t));
    f.read(reinterpret_cast<char*>(&length),   sizeof(uint32_t));

    s = nm_yale_storage_create(dtype, shape, dim, length); // set length as init capacity

    nm_register_storage(stype, s);

    read_padded_yale_elements(f, reinterpret_cast<YALE_STORAGE*>(s), length, symm, dtype);
  } else {
    NM_CONSERVATIVE(nm_unregister_values(argv, argc));
    NM_CONSERVATIVE(nm_unregister_value(&self));
    rb_raise(nm_eStorageTypeError, "please convert to yale or dense before saving");
  }

  NMATRIX* nm = nm_create(stype, s);

  // Return the appropriate matrix object (Ruby VALUE)
  // FIXME: This should probably return CLASS_OF(self) instead of cNMatrix, but I don't know how that works for
  // FIXME: class methods.
  nm_register_nmatrix(nm);
  VALUE to_return = Data_Wrap_Struct(cNMatrix, nm_mark, nm_delete, nm);

  nm_unregister_nmatrix(nm);
  NM_CONSERVATIVE(nm_unregister_values(argv, argc));
  NM_CONSERVATIVE(nm_unregister_value(&self));
  nm_unregister_storage(stype, s);

  switch(stype) {
  case nm::DENSE_STORE:
  case nm::YALE_STORE:
    return to_return;
  default: // this case never occurs (due to earlier rb_raise)
    return Qnil;
  }

}

# File lib/nmatrix/lapack_ext_common.rb, line 30
def NMatrix.register_lapack_extension(name)
  if (defined? @@lapack_extension)
    raise "Attempting to load #{name} when #{@@lapack_extension} is already loaded. You can only load one LAPACK extension."
  end

  @@lapack_extension = name
end

# File lib/nmatrix/shortcuts.rb, line 418
def seq(shape, options={})

  # Construct the values of the final matrix based on the dimension.
  values = (0 ... NMatrix.size(shape)).to_a

  # It'll produce :int32, except if a dtype is provided.
  NMatrix.new(shape, values, {:stype => :dense}.merge(options))
end

# File lib/nmatrix/nmatrix.rb, line 104
def size(shape)
  shape = [shape,shape] unless shape.is_a?(Array)
  (0...shape.size).inject(1) { |x,i| x * shape[i] }
end

# File lib/nmatrix/homogeneous.rb, line 127
def translation *args
  xyz = args.shift if args.first.is_a?(NMatrix) || args.first.is_a?(Array)
  default_dtype = xyz.respond_to?(:dtype) ? xyz.dtype : NMatrix.guess_dtype(xyz)
  opts = {dtype: default_dtype}
  opts = opts.merge(args.pop) if args.size > 0 && args.last.is_a?(Hash)
  xyz ||= args

  n = if args.size > 0
    NMatrix.eye(4, opts)
  else
    NMatrix.eye(4, opts)
  end
  n[0..2,3] = xyz
  n
end

static VALUE nm_upcast(VALUE self, VALUE t1, VALUE t2) {
  nm::dtype_t d1    = nm_dtype_from_rbsymbol(t1),
              d2    = nm_dtype_from_rbsymbol(t2);

  return ID2SYM(rb_intern( DTYPE_NAMES[ Upcast[d1][d2] ] ));
}

# File lib/nmatrix/homogeneous.rb, line 66
def x_rotation angle_in_radians, opts={}
  c = Math.cos(angle_in_radians)
  s = Math.sin(angle_in_radians)
  NMatrix.new(4, [1.0, 0.0, 0.0, 0.0,
                  0.0, c,   -s,  0.0,
                  0.0, s,    c,  0.0,
                  0.0, 0.0, 0.0, 1.0], {dtype: :float64}.merge(opts))
end

# File lib/nmatrix/homogeneous.rb, line 75
def y_rotation angle_in_radians, opts={}
  c = Math.cos(angle_in_radians)
  s = Math.sin(angle_in_radians)
  NMatrix.new(4, [ c,  0.0,  s,  0.0,
                  0.0, 1.0, 0.0, 0.0,
                  -s,  0.0,  c,  0.0,
                  0.0, 0.0, 0.0, 1.0], {dtype: :float64}.merge(opts))
end

# File lib/nmatrix/homogeneous.rb, line 84
def z_rotation angle_in_radians, opts={}
  c = Math.cos(angle_in_radians)
  s = Math.sin(angle_in_radians)
  NMatrix.new(4, [ c,  -s,  0.0, 0.0,
                   s,   c,  0.0, 0.0,
                  0.0, 0.0, 1.0, 0.0,
                  0.0, 0.0, 0.0, 1.0], {dtype: :float64}.merge(opts))
end

# File lib/nmatrix/shortcuts.rb, line 146
def zeros(shape, opts = {})
  NMatrix.new(shape, 0, {:dtype => :float64}.merge(opts))
end

# File lib/nmatrix/shortcuts.rb, line 197
def zeros_like(nm)
  NMatrix.zeros(nm.shape, dtype: nm.dtype, stype: nm.stype, capacity: nm.capacity, default: 0)
end

static VALUE nm_eqeq(VALUE left, VALUE right) {
  NM_CONSERVATIVE(nm_register_value(&left));
  NM_CONSERVATIVE(nm_register_value(&right));

  NMATRIX *l, *r;

  CheckNMatrixType(left);
  CheckNMatrixType(right);

  UnwrapNMatrix(left, l);
  UnwrapNMatrix(right, r);

  bool result = false;

  // Check that the shapes match before going any further.
  if (l->storage->dim != r->storage->dim) {
    NM_CONSERVATIVE(nm_unregister_value(&left));
    NM_CONSERVATIVE(nm_unregister_value(&right));
    rb_raise(nm_eShapeError, "cannot compare matrices with different dimension");
  }

  size_t dim = l->storage->dim;
  for (size_t i=0; i<dim; i++) {
    if (l->storage->shape[i] != r->storage->shape[i]) {
      NM_CONSERVATIVE(nm_unregister_value(&left));
      NM_CONSERVATIVE(nm_unregister_value(&right));
      rb_raise(nm_eShapeError, "cannot compare matrices with different shapes");
    }
  }

  if (l->stype != r->stype) { // DIFFERENT STYPES

    if (l->stype == nm::DENSE_STORE)
      result = rb_funcall(left, rb_intern("dense_eql_sparse?"), 1, right);
    else if (r->stype == nm::DENSE_STORE)
      result = rb_funcall(right, rb_intern("dense_eql_sparse?"), 1, left);
    else
      result = rb_funcall(left, rb_intern("sparse_eql_sparse?"), 1, right);

  } else {

    switch(l->stype) {       // SAME STYPES
    case nm::DENSE_STORE:
      result = nm_dense_storage_eqeq(l->storage, r->storage);
      break;
    case nm::LIST_STORE:
      result = nm_list_storage_eqeq(l->storage, r->storage);
      break;
    case nm::YALE_STORE:
      result = nm_yale_storage_eqeq(l->storage, r->storage);
      break;
    }
  }

  NM_CONSERVATIVE(nm_unregister_value(&left));
  NM_CONSERVATIVE(nm_unregister_value(&right));

  return result ? Qtrue : Qfalse;
}

static VALUE nm_mref(int argc, VALUE* argv, VALUE self) {
  static void* (*ttable[nm::NUM_STYPES])(const STORAGE*, SLICE*) = {
    nm_dense_storage_ref,
    nm_list_storage_ref,
    nm_yale_storage_ref
  };
  nm::stype_t stype = NM_STYPE(self);
  return nm_xslice(argc, argv, ttable[stype], nm_delete_ref, self);
}

static VALUE nm_mset(int argc, VALUE* argv, VALUE self) {

  size_t dim = NM_DIM(self); // last arg is the value

  VALUE to_return = Qnil;

  if ((size_t)(argc) > NM_DIM(self)+1) {
    rb_raise(rb_eArgError, "wrong number of arguments (%d for %lu)", argc, effective_dim(NM_STORAGE(self))+1);
  } else {
    NM_CONSERVATIVE(nm_register_value(&self));
    NM_CONSERVATIVE(nm_register_values(argv, argc));

    SLICE* slice = get_slice(dim, argc-1, argv, NM_STORAGE(self)->shape);

    static void (*ttable[nm::NUM_STYPES])(VALUE, SLICE*, VALUE) = {
      nm_dense_storage_set,
      nm_list_storage_set,
      nm_yale_storage_set
    };

    ttable[NM_STYPE(self)](self, slice, argv[argc-1]);

    free_slice(slice);

    to_return = argv[argc-1];

    NM_CONSERVATIVE(nm_unregister_value(&self));
    NM_CONSERVATIVE(nm_unregister_values(argv, argc));
  }

  return to_return;
}

# File lib/nmatrix/math.rb, line 784
def abs
  if stype == :dense
    self.__dense_map__ { |v| v.abs }
  elsif stype == :list
    # FIXME: Need __list_map_stored__, but this will do for now.
    self.__list_map_merged_stored__(nil, nil) { |v,dummy| v.abs }
  else
    self.__yale_map_stored__ { |v| v.abs }
  end.cast(self.stype, abs_dtype)
end

# File lib/nmatrix/math.rb, line 768
def abs_dtype
  if self.dtype == :complex64
    :float32
  elsif self.dtype == :complex128
    :float64
  else
    self.dtype
  end
end

# File lib/nmatrix/homogeneous.rb, line 228
def angle_vector
  raise(ShapeError, "Expected length-4 vector or matrix (quaternion)") if self.shape[0] != 4
  raise("Expected unit quaternion") if self[0] > 1

  xyz = NMatrix.new([3], dtype: self.dtype)

  angle = 2 * Math.acos(self[0])
  s = Math.sqrt(1.0 - self[0]*self[0])

  xyz[0..2] = self[1..3]
  xyz /= s if s >= 0.001 # avoid divide by zero
  return [angle, xyz]
end

# File lib/nmatrix/math.rb, line 805
def asum incx=1, n=nil
  if self.shape == [1]
    return self[0].abs unless self.complex_dtype?
    return self[0].real.abs + self[0].imag.abs
  end
  return method_missing(:asum, incx, n) unless vector?
  NMatrix::BLAS::asum(self, incx, self.size / incx)
end

# File lib/nmatrix/nmatrix.rb, line 913
def binned_sorted_indices
  return method_missing(:sorted_indices) unless vector?
  ary = self.to_flat_array
  ary2 = []
  last_bin = ary.each_index.sort_by { |i| [ary[i]] }.inject([]) do |result, element|
    if result.empty? || ary[result[-1]] == ary[element]
      result << element
    else
      ary2 << result
      [element]
    end
  end
  ary2 << last_bin unless last_bin.empty?
  ary2
end

static VALUE nm_capacity(VALUE self) {
  NM_CONSERVATIVE(nm_register_value(&self));
  VALUE cap;

  switch(NM_STYPE(self)) {
  case nm::YALE_STORE:
    cap = UINT2NUM(reinterpret_cast<YALE_STORAGE*>(NM_STORAGE_YALE(self)->src)->capacity);
    break;

  case nm::DENSE_STORE:
    cap = UINT2NUM(nm_storage_count_max_elements( NM_STORAGE_DENSE(self) ));
    break;

  case nm::LIST_STORE:
    cap = UINT2NUM(nm_list_storage_count_elements( NM_STORAGE_LIST(self) ));
    break;

  default:
    NM_CONSERVATIVE(nm_unregister_value(&self));
    rb_raise(nm_eStorageTypeError, "unrecognized stype in nm_capacity()");
  }

  NM_CONSERVATIVE(nm_unregister_value(&self));
  return cap;
}

# File lib/nmatrix/nmatrix.rb, line 226
def cast(*params)
  if (params.size > 0 && params[0].is_a?(Hash))
    opts = {
        :stype => self.stype,
        :dtype => self.dtype,
        :default => self.stype == :dense ? 0 : self.default_value
    }.merge(params[0])

    self.cast_full(opts[:stype], opts[:dtype], opts[:default])
  else
    params << self.stype if params.size == 0
    params << self.dtype if params.size == 1
    #HACK: the default value can cause an exception if dtype is not complex
    #and default_value is. (The ruby C code apparently won't convert these.)
    #Perhaps this should be fixed in the C code (in rubyval_to_cval).
    default_value = maybe_get_noncomplex_default_value(params[1])
    params << (self.stype == :dense ? 0 : default_value) if params.size == 2
    self.cast_full(*params)
  end

end

VALUE nm_cast(VALUE self, VALUE new_stype_symbol, VALUE new_dtype_symbol, VALUE init) {
  NM_CONSERVATIVE(nm_register_value(&self));
  NM_CONSERVATIVE(nm_register_value(&init));

  nm::dtype_t new_dtype = nm_dtype_from_rbsymbol(new_dtype_symbol);
  nm::stype_t new_stype = nm_stype_from_rbsymbol(new_stype_symbol);

  CheckNMatrixType(self);
  NMATRIX *rhs;

  UnwrapNMatrix( self, rhs );

  void* init_ptr = NM_ALLOCA_N(char, DTYPE_SIZES[new_dtype]);
  rubyval_to_cval(init, new_dtype, init_ptr);

  NMATRIX* m = nm_cast_with_ctype_args(rhs, new_stype, new_dtype, init_ptr);
  nm_register_nmatrix(m);

  VALUE to_return = Data_Wrap_Struct(CLASS_OF(self), nm_mark, nm_delete, m);

  nm_unregister_nmatrix(m);
  NM_CONSERVATIVE(nm_unregister_value(&self));
  NM_CONSERVATIVE(nm_unregister_value(&init));
  return to_return;

}

# File lib/nmatrix/nmatrix.rb, line 989
def clone_structure(capacity = nil)
  opts = {stype: self.stype, default: self.default_value, dtype: self.dtype}
  opts = {capacity: capacity}.merge(opts) if self.yale?
  NMatrix.new(self.shape, opts)
end

# File lib/nmatrix/nmatrix.rb, line 267
def cols
  shape[1]
end

# File lib/nmatrix/nmatrix.rb, line 516
def column(column_number, get_by = :copy)
  rank(1, column_number, get_by)
end

static VALUE nm_complex_conjugate(VALUE self) {
  VALUE copy;
  return nm_complex_conjugate_bang(nm_init_copy(copy,self));
}

static VALUE nm_complex_conjugate_bang(VALUE self) {

  NMATRIX* m;
  void* elem;
  size_t size, p;

  UnwrapNMatrix(self, m);

  if (m->stype == nm::DENSE_STORE) {

    size = nm_storage_count_max_elements(NM_STORAGE(self));
    elem = NM_STORAGE_DENSE(self)->elements;

  } else if (m->stype == nm::YALE_STORE) {

    size = nm_yale_storage_get_size(NM_STORAGE_YALE(self));
    elem = NM_STORAGE_YALE(self)->a;

  } else {
    rb_raise(rb_eNotImpError, "please cast to yale or dense (complex) first");
  }

  // Walk through and negate the imaginary component
  if (NM_DTYPE(self) == nm::COMPLEX64) {

    for (p = 0; p < size; ++p) {
      reinterpret_cast<nm::Complex64*>(elem)[p].i = -reinterpret_cast<nm::Complex64*>(elem)[p].i;
    }

  } else if (NM_DTYPE(self) == nm::COMPLEX128) {

    for (p = 0; p < size; ++p) {
      reinterpret_cast<nm::Complex128*>(elem)[p].i = -reinterpret_cast<nm::Complex128*>(elem)[p].i;
    }

  }
  return self;
}

# File lib/nmatrix/nmatrix.rb, line 364
def complex_dtype?
  [:complex64, :complex128].include?(self.dtype)
end

# File lib/nmatrix/nmatrix.rb, line 662
def concat(*matrices)
  rank = nil
  rank = matrices.pop unless matrices.last.is_a?(NMatrix)

  # Find the first matching dimension and concatenate along that (unless rank is specified)
  if rank.nil?
    rank = self.dim-1
    self.shape.reverse_each.with_index do |s,i|
      matrices.each do |m|
        if m.shape[i] != s
          rank -= 1
          break
        end
      end
    end
  elsif rank.is_a?(Symbol) # Convert to numeric
    rank = {:row => 0, :column => 1, :col => 1, :lay => 2, :layer => 2}[rank]
  end

  # Need to figure out the new shape.
  new_shape = self.shape.dup
  new_shape[rank] = matrices.inject(self.shape[rank]) { |total,m| total + m.shape[rank] }

  # Now figure out the options for constructing the concatenated matrix.
  opts = {stype: self.stype, default: self.default_value, dtype: self.dtype}
  if self.yale?
    # We can generally predict the new capacity for Yale. Subtract out the number of rows
    # for each matrix being concatenated, and then add in the number of rows for the new
    # shape. That takes care of the diagonal. The rest of the capacity is represented by
    # the non-diagonal non-default values.
    new_cap = matrices.inject(self.capacity - self.shape[0]) do |total,m|
      total + m.capacity - m.shape[0]
    end - self.shape[0] + new_shape[0]
    opts = {capacity: new_cap}.merge(opts)
  end

  # Do the actual construction.
  n = NMatrix.new(new_shape, opts)

  # Figure out where to start and stop the concatenation. We'll use NMatrices instead of
  # Arrays because then we can do elementwise addition.
  ranges = self.shape.map.with_index { |s,i| 0...self.shape[i] }

  matrices.unshift(self)
  matrices.each do |m|
    n[*ranges] = m

    # move over by the requisite amount
    ranges[rank]  = (ranges[rank].first + m.shape[rank])...(ranges[rank].last + m.shape[rank])
  end

  n
end

# File lib/nmatrix/math.rb, line 624
def conjugate_transpose
  self.transpose.complex_conjugate!
end

# File lib/nmatrix/math.rb, line 523
def corr
  raise NotImplementedError, "Does not work for complex dtypes" if complex_dtype?
  standard_deviation = std
  cov / (standard_deviation.transpose.dot(standard_deviation))
end

# File lib/nmatrix/math.rb, line 510
def cov(opts={})
  raise TypeError, "Only works for non-integer dtypes" if integer_dtype?
   opts = {
    for_sample_data: true
  }.merge(opts)
  
  denominator      = opts[:for_sample_data] ? rows - 1 : rows
  ones             = NMatrix.ones [rows,1] 
  deviation_scores = self - ones.dot(ones.transpose).dot(self) / rows
  deviation_scores.transpose.dot(deviation_scores) / denominator
end

static VALUE nm_data_pointer(VALUE self) {
  //if (NM_DTYPE(self) == nm::LIST_STORE)
  //  rb_warn("pointer requested for list storage, which may be meaningless");

  // This is actually pretty easy, since all of the storage types have their elements positioned in the same place
  // relative to one another. So yes, believe it or not, this should work just as well for Yale or list storage as for
  // dense.
  return INT2FIX(NM_STORAGE_DENSE(self)->elements);
}

# File lib/nmatrix/nmatrix.rb, line 727
def dconcat(*matrices)
  concat(*matrices, :layer)
end

static VALUE nm_default_value(VALUE self) {
  switch(NM_STYPE(self)) {
  case nm::YALE_STORE:
    return nm_yale_default_value(self);
  case nm::LIST_STORE:
    return nm_list_default_value(self);
  case nm::DENSE_STORE:
  default:
    return Qnil;
  }
}

# File lib/nmatrix/shortcuts.rb, line 41
def dense?; return stype == :dense; end

# File lib/nmatrix/math.rb, line 454
def det
  raise(ShapeError, "determinant can be calculated only for square matrices") unless self.dim == 2 && self.shape[0] == self.shape[1]

  # Cast to a dtype for which getrf is implemented
  new_dtype = self.integer_dtype? ? :float64 : self.dtype
  copy = self.cast(:dense, new_dtype)

  # Need to know the number of permutations. We'll add up the diagonals of
  # the factorized matrix.
  pivot = copy.getrf!

  num_perm = 0 #number of permutations
  pivot.each_with_index do |swap, i|
    #pivot indexes rows starting from 1, instead of 0, so need to subtract 1 here
    num_perm += 1 if swap-1 != i
  end
  prod = num_perm % 2 == 1 ? -1 : 1 # odd permutations => negative
  [shape[0],shape[1]].min.times do |i|
    prod *= copy[i,i]
  end

  # Convert back to an integer if necessary
  new_dtype != self.dtype ? prod.round : prod #prevent rounding errors
end

static VALUE nm_det_exact(VALUE self) {

  if (NM_STYPE(self) != nm::DENSE_STORE) {
    rb_raise(rb_eNotImpError, "can only calculate exact determinant for dense matrices");
    return Qnil;
  }
  if (NM_DIM(self) != 2 || NM_SHAPE0(self) != NM_SHAPE1(self)) {
    rb_raise(nm_eShapeError, "matrices must be square to have a determinant defined");
    return Qnil;
  }

  NM_CONSERVATIVE(nm_register_value(&self));

  // Calculate the determinant and then assign it to the return value
  void* result = NM_ALLOCA_N(char, DTYPE_SIZES[NM_DTYPE(self)]);
  nm::dtype_t dtype = NM_DTYPE(self);
  nm_math_det_exact(NM_SHAPE0(self), NM_STORAGE_DENSE(self)->elements, NM_SHAPE0(self), NM_DTYPE(self), result);

  if (dtype == nm::RUBYOBJ) {
    nm_register_values(reinterpret_cast<VALUE*>(result), 1);
  }
  VALUE to_return = rubyobj_from_cval(result, NM_DTYPE(self)).rval;
  if (dtype == nm::RUBYOBJ) {
    nm_unregister_values(reinterpret_cast<VALUE*>(result), 1);
  }
  NM_CONSERVATIVE(nm_unregister_value(&self));

  return to_return;
}

# File lib/nmatrix/nmatrix.rb, line 281
def diagonal main_diagonal=true
  diag_size = [cols, rows].min
  diag = NMatrix.new [diag_size], dtype: dtype

  if main_diagonal
    0.upto(diag_size-1) do |i|
      diag[i] = self[i,i]
    end
  else
    row = 0
    (diag_size-1).downto(0) do |col|
      diag[row] = self[row,col]
      row += 1
    end
  end

  diag
end

static VALUE nm_dim(VALUE self) {
  return INT2FIX(NM_STORAGE(self)->dim);
}

static VALUE nm_multiply(VALUE left_v, VALUE right_v) {
  NM_CONSERVATIVE(nm_register_value(&left_v));
  NM_CONSERVATIVE(nm_register_value(&right_v));

  NMATRIX *left, *right;

  UnwrapNMatrix( left_v, left );

  if (NM_RUBYVAL_IS_NUMERIC(right_v)) {
    NM_CONSERVATIVE(nm_unregister_value(&left_v));
    NM_CONSERVATIVE(nm_unregister_value(&right_v));
    return matrix_multiply_scalar(left, right_v);
  }

  else if (TYPE(right_v) == T_ARRAY) {
    NM_CONSERVATIVE(nm_unregister_value(&left_v));
    NM_CONSERVATIVE(nm_unregister_value(&right_v));
    rb_raise(rb_eNotImpError, "please convert array to nx1 or 1xn NMatrix first");
  }

  else { // both are matrices (probably)
    CheckNMatrixType(right_v);
    UnwrapNMatrix( right_v, right );

    // work like vector dot product for 1dim
    if (left->storage->dim == 1 && right->storage->dim == 1) {
      if (left->storage->shape[0] != right->storage->shape[0]) {
        NM_CONSERVATIVE(nm_unregister_value(&left_v));
        NM_CONSERVATIVE(nm_unregister_value(&right_v));
        rb_raise(rb_eArgError, "The left- and right-hand sides of the operation must have the same dimensionality.");
      } else {
        VALUE result = elementwise_op(nm::EW_MUL, left_v, right_v);
        VALUE to_return = rb_funcall(result, rb_intern("sum"),0);
        NM_CONSERVATIVE(nm_unregister_value(&left_v));
        NM_CONSERVATIVE(nm_unregister_value(&right_v));
        return to_return;
      }
    }

    if (left->storage->shape[1] != right->storage->shape[0]) {
      NM_CONSERVATIVE(nm_unregister_value(&left_v));
      NM_CONSERVATIVE(nm_unregister_value(&right_v));
      rb_raise(rb_eArgError, "incompatible dimensions");
    }

    if (left->stype != right->stype) {
      NM_CONSERVATIVE(nm_unregister_value(&left_v));
      NM_CONSERVATIVE(nm_unregister_value(&right_v));
      rb_raise(rb_eNotImpError, "matrices must have same stype");
    }

    NM_CONSERVATIVE(nm_unregister_value(&left_v));
    NM_CONSERVATIVE(nm_unregister_value(&right_v));
    return matrix_multiply(left, right);

  }

  NM_CONSERVATIVE(nm_unregister_value(&left_v));
  NM_CONSERVATIVE(nm_unregister_value(&right_v));

  return Qnil;
}

static VALUE nm_dtype(VALUE self) {
  ID dtype = rb_intern(DTYPE_NAMES[NM_DTYPE(self)]);
  return ID2SYM(dtype);
}

# File lib/nmatrix/enumerate.rb, line 40
def each &bl
  if self.stype == :dense
    self.__dense_each__(&bl)
  elsif block_given?
    self.each_with_indices(&bl)
  else # Handle case where no block is given
    Enumerator.new do |yielder|
      self.each_with_indices do |params|
        yielder.yield params
      end
    end
  end
end

# File lib/nmatrix/enumerate.rb, line 144
def each_column(get_by=:reference)
  return enum_for(:each_column, get_by) unless block_given?
  (0...self.shape[1]).each do |j|
    yield self.column(j, get_by)
  end
  self
end

# File lib/nmatrix/enumerate.rb, line 160
def each_layer(get_by=:reference)
  return enum_for(:each_layer, get_by) unless block_given?
  (0...self.shape[2]).each do |k|
    yield self.layer(k, get_by)
  end
  self
end

static VALUE nm_each_ordered_stored_with_indices(VALUE nmatrix) {
  NM_CONSERVATIVE(nm_register_value(&nmatrix));
  VALUE to_return = Qnil;

  switch(NM_STYPE(nmatrix)) {
  case nm::YALE_STORE:
    to_return = nm_yale_each_ordered_stored_with_indices(nmatrix);
    break;
  case nm::DENSE_STORE:
    to_return = nm_dense_each_with_indices(nmatrix);
    break;
  case nm::LIST_STORE:
    to_return = nm_list_each_with_indices(nmatrix, true);
    break;
  default:
    NM_CONSERVATIVE(nm_unregister_value(&nmatrix));
    rb_raise(nm_eDataTypeError, "Not a proper storage type");
  }

  NM_CONSERVATIVE(nm_unregister_value(&nmatrix));
  return to_return;
}

# File lib/nmatrix/enumerate.rb, line 117
def each_rank(dimen=0, get_by=:reference)
  return enum_for(:each_rank, dimen, get_by) unless block_given?
  (0...self.shape[dimen]).each do |idx|
    yield self.rank(dimen, idx, get_by)
  end
  self
end

# File lib/nmatrix/enumerate.rb, line 131
def each_row(get_by=:reference)
  return enum_for(:each_row, get_by) unless block_given?
  (0...self.shape[0]).each do |i|
    yield self.row(i, get_by)
  end
  self
end

# File lib/nmatrix/enumerate.rb, line 175
def each_stored_with_index(&block)
  raise(NotImplementedError, "only works for dim 2 vectors") unless self.dim <= 2
  return enum_for(:each_stored_with_index) unless block_given?

  self.each_stored_with_indices do |v, i, j|
    if shape[0] == 1
      yield(v,j)
    elsif shape[1] == 1
      yield(v,i)
    else
      method_missing(:each_stored_with_index, &block)
    end
  end
  self
end

static VALUE nm_each_stored_with_indices(VALUE nmatrix) {
  NM_CONSERVATIVE(nm_register_value(&nmatrix));
  VALUE to_return = Qnil;

  switch(NM_STYPE(nmatrix)) {
  case nm::YALE_STORE:
    to_return = nm_yale_each_stored_with_indices(nmatrix);
    break;
  case nm::DENSE_STORE:
    to_return = nm_dense_each_with_indices(nmatrix);
    break;
  case nm::LIST_STORE:
    to_return = nm_list_each_with_indices(nmatrix, true);
    break;
  default:
    NM_CONSERVATIVE(nm_unregister_value(&nmatrix));
    rb_raise(nm_eDataTypeError, "Not a proper storage type");
  }

  NM_CONSERVATIVE(nm_unregister_value(&nmatrix));
  return to_return;
}

static VALUE nm_each_with_indices(VALUE nmatrix) {
  NM_CONSERVATIVE(nm_register_value(&nmatrix));
  VALUE to_return = Qnil;

  switch(NM_STYPE(nmatrix)) {
  case nm::YALE_STORE:
    to_return = nm_yale_each_with_indices(nmatrix);
    break;
  case nm::DENSE_STORE:
    to_return = nm_dense_each_with_indices(nmatrix);
    break;
  case nm::LIST_STORE:
    to_return = nm_list_each_with_indices(nmatrix, false);
    break;
  default:
    NM_CONSERVATIVE(nm_unregister_value(&nmatrix));
    rb_raise(nm_eDataTypeError, "Not a proper storage type");
  }

  NM_CONSERVATIVE(nm_unregister_value(&nmatrix));
  return to_return;
}

static VALUE nm_effective_dim(VALUE self) {
  return INT2FIX(effective_dim(NM_STORAGE(self)));
}

# File lib/nmatrix/math.rb, line 204
def factorize_cholesky
  raise "Matrix must be symmetric/Hermitian for Cholesky factorization" unless self.hermitian?
  l = self.clone.potrf_lower!.tril!
  u = l.conjugate_transpose
  [u,l]
end

# File lib/nmatrix/math.rb, line 224
def factorize_lu with_permutation_matrix=nil
  raise(NotImplementedError, "only implemented for dense storage") unless self.stype == :dense
  raise(NotImplementedError, "matrix is not 2-dimensional") unless self.dimensions == 2

  t     = self.clone
  pivot = t.getrf!
  return t unless with_permutation_matrix

  [t, FactorizeLUMethods.permutation_matrix_from(pivot)]
end

# File lib/nmatrix/nmatrix.rb, line 354
def float_dtype?
  [:float32, :float64].include?(dtype)
end

# File lib/nmatrix/math.rb, line 350
def gesdd(workspace_size=nil)
  self.clone.gesdd!(workspace_size)
end

# File lib/nmatrix/math.rb, line 335
def gesdd!(workspace_size=nil)
  NMatrix::LAPACK::gesdd(self, workspace_size)
end

# File lib/nmatrix/math.rb, line 318
def gesvd(workspace_size=1)
  self.clone.gesvd!(workspace_size)
end

# File lib/nmatrix/math.rb, line 304
def gesvd!(workspace_size=1)
  NMatrix::LAPACK::gesvd(self, workspace_size)
end

# File lib/nmatrix/lapacke.rb, line 169
def getrf!
  raise(StorageTypeError, "LAPACK functions only work on dense matrices") unless self.dense?

  ipiv = NMatrix::LAPACK::lapacke_getrf(:row, self.shape[0], self.shape[1], self, self.shape[1])

  return ipiv
end

# File lib/nmatrix/nmatrix.rb, line 717
def hconcat(*matrices)
  concat(*matrices, :column)
end

static VALUE nm_hermitian(VALUE self) {
  return is_symmetric(self, true);
}

# File lib/nmatrix/math.rb, line 241
def hessenberg
  clone.hessenberg!
end

# File lib/nmatrix/math.rb, line 246
def hessenberg!
  raise ShapeError, "Trying to reduce non 2D matrix to hessenberg form" if 
    shape.size != 2
  raise ShapeError, "Trying to reduce non-square matrix to hessenberg form" if 
    shape[0] != shape[1]
  raise StorageTypeError, "Matrix must be dense" if stype != :dense
  raise TypeError, "Works with float matrices only" unless 
    [:float64,:float32].include?(dtype)

  __hessenberg__(self)
  self
end

# File lib/nmatrix/nmatrix.rb, line 966
def index(value)
  index = nil

  self.each_with_indices do |yields|
    if yields.first == value
      yields.shift
      index = yields
      break
    end
  end 

  index
end

static VALUE nm_init_copy(VALUE copy, VALUE original) {
  NM_CONSERVATIVE(nm_register_value(&copy));
  NM_CONSERVATIVE(nm_register_value(&original));

  NMATRIX *lhs, *rhs;

  CheckNMatrixType(original);

  if (copy == original) {
    NM_CONSERVATIVE(nm_unregister_value(&copy));
    NM_CONSERVATIVE(nm_unregister_value(&original));
    return copy;
  }

  UnwrapNMatrix( original, rhs );
  UnwrapNMatrix( copy,     lhs );

  lhs->stype = rhs->stype;

  // Copy the storage
  CAST_TABLE(ttable);
  lhs->storage = ttable[lhs->stype][rhs->stype](rhs->storage, rhs->storage->dtype, NULL);

  NM_CONSERVATIVE(nm_unregister_value(&copy));
  NM_CONSERVATIVE(nm_unregister_value(&original));

  return copy;
}

# File lib/nmatrix/nmatrix.rb, line 958
def inject(sym)
  return super(sym) unless self.yale?
  return self.map_stored.inject(sym)
end

# File lib/nmatrix/enumerate.rb, line 217
def inject_rank(dimen=0, initial=nil, dtype=nil)

  raise(RangeError, "requested dimension (#{dimen}) does not exist (shape: #{shape})") if dimen > self.dim

  return enum_for(:inject_rank, dimen, initial, dtype) unless block_given?

  new_shape = shape
  new_shape[dimen] = 1

  first_as_acc = false

  if initial then
    acc = NMatrix.new(new_shape, initial, :dtype => dtype || self.dtype, stype: self.stype)
  else
    each_rank(dimen) do |sub_mat|
      acc = (sub_mat.is_a?(NMatrix) and !dtype.nil? and dtype != self.dtype) ? sub_mat.cast(self.stype, dtype) : sub_mat
      break
    end
    first_as_acc = true
  end

  each_rank(dimen) do |sub_mat|
    if first_as_acc
      first_as_acc = false
      next
    end
    acc = yield(acc, sub_mat)
  end

  acc
end

# File lib/nmatrix/nmatrix.rb, line 345
def integer_dtype?
  [:byte, :int8, :int16, :int32, :int64].include?(self.dtype)
end

# File lib/nmatrix/math.rb, line 102
def invert
  #write this in terms of invert! so plugins will only have to overwrite
  #invert! and not invert
  if self.integer_dtype?
    cloned = self.cast(dtype: :float64)
    cloned.invert!
  else
    cloned = self.clone
    cloned.invert!
  end
end

# File lib/nmatrix/atlas.rb, line 186
def invert!
  raise(StorageTypeError, "invert only works on dense matrices currently") unless self.dense?
  raise(ShapeError, "Cannot invert non-square matrix") unless shape[0] == shape[1]
  raise(DataTypeError, "Cannot invert an integer matrix in-place") if self.integer_dtype?

  # Even though we are using the ATLAS plugin, we still might be missing
  # CLAPACK (and thus clapack_getri) if we are on OS X.
  if NMatrix.has_clapack?
    # Get the pivot array; factor the matrix
    # We can't used getrf! here since it doesn't have the clapack behavior,
    # so it doesn't play nicely with clapack_getri
    n = self.shape[0]
    pivot = NMatrix::LAPACK::clapack_getrf(:row, n, n, self, n)
    # Now calculate the inverse using the pivot array
    NMatrix::LAPACK::clapack_getri(:row, n, self, n, pivot)
    self
  else
    __inverse__(self,true)
  end
end

static VALUE nm_is_ref(VALUE self) {
  if (NM_SRC(self) == NM_STORAGE(self)) return Qfalse;
  return Qtrue;
}

# File lib/nmatrix/math.rb, line 587
def kron_prod(mat)
  unless self.dimensions==2 and mat.dimensions==2
    raise ShapeError, "Implemented for 2D NMatrix objects only."
  end

  # compute the shape [n,m] of the product matrix
  n, m = self.shape[0]*mat.shape[0], self.shape[1]*mat.shape[1]
  # compute the entries of the product matrix
  kron_prod_array = []
  if self.yale?
    # +:yale+ requires to get the row by copy in order to apply +#transpose+ to it
    self.each_row(getby=:copy) do |selfr|
      mat.each_row do |matr|
        kron_prod_array += (selfr.transpose.dot matr).to_flat_a
      end
    end
  else
    self.each_row do |selfr|
      mat.each_row do |matr|
        kron_prod_array += (selfr.transpose.dot matr).to_flat_a
      end
    end
  end

  NMatrix.new([n,m], kron_prod_array) 
end

# File lib/nmatrix/math.rb, line 425
def laswp(ary, opts={})
  self.clone.laswp!(ary, opts)
end

# File lib/nmatrix/math.rb, line 378
def laswp!(ary, opts={})
  raise(StorageTypeError, "ATLAS functions only work on dense matrices") unless self.dense?
  opts = { convention: :intuitive }.merge(opts)
  
  if opts[:convention] == :intuitive
    if ary.length != ary.uniq.length
      raise(ArgumentError, "No duplicated entries in the order array are allowed under convention :intuitive")
    end
    n = self.shape[1]
    p = []
    order = (0...n).to_a
    0.upto(n-2) do |i|
      p[i] = order.index(ary[i])
      order[i], order[p[i]] = order[p[i]], order[i]
    end
    p[n-1] = n-1
  else
    p = ary
  end

  NMatrix::LAPACK::laswp(self, p)
end

# File lib/nmatrix/nmatrix.rb, line 853
def layer(layer_number, get_by = :copy)
  rank(2, layer_number, get_by)
end

# File lib/nmatrix/shortcuts.rb, line 53
def list?;  return stype == :list; end

static VALUE nm_unary_log(int argc, VALUE* argv, VALUE self) {
  NM_CONSERVATIVE(nm_register_values(argv, argc));
  const double default_log_base = exp(1.0);
  NMATRIX* left;
  UnwrapNMatrix(self, left);
  std::string sym;

  switch(left->stype) {
  case nm::DENSE_STORE:
    sym = "__dense_unary_log__";
    break;
  case nm::YALE_STORE:
    sym = "__yale_unary_log__";
    break;
  case nm::LIST_STORE:
    sym = "__list_unary_log__";
    break;
  }
  NM_CONSERVATIVE(nm_unregister_values(argv, argc));
  if (argc > 0) { //supplied a base
    return rb_funcall(self, rb_intern(sym.c_str()), 1, argv[0]);
  }
  return rb_funcall(self, rb_intern(sym.c_str()), 1, nm::RubyObject(default_log_base).rval);
}

# File lib/nmatrix/nmatrix.rb, line 800
def lower_triangle(k = 0)
  raise(NotImplementedError, "only implemented for 2D matrices") if self.shape.size > 2

  t = self.clone_structure
  (0...self.shape[0]).each do |i|
    if i + k >= shape[0]
      t[i, :*] = self[i, :*]
    else
      t[i, (i+k+1)...self.shape[1]] = 0
      t[i, 0..(i+k)] = self[i, 0..(i+k)]
    end
  end
  t
end

# File lib/nmatrix/nmatrix.rb, line 829
def lower_triangle!(k = 0)
  raise(NotImplementedError, "only implemented for 2D matrices") if self.shape.size > 2

  (0...self.shape[0]).each do |i|
    if i + k < shape[0]
      self[i, (i+k+1)...self.shape[1]] = 0
    end
  end
  self
end

# File lib/nmatrix/enumerate.rb, line 72
def map(&bl)
  return enum_for(:map) unless block_given?
  cp = self.cast(dtype: :object)
  cp.map!(&bl)
  cp
end

# File lib/nmatrix/enumerate.rb, line 87
def map!
  return enum_for(:map!) unless block_given?
  iterated = false
  self.each_stored_with_indices do |e, *i|
    iterated = true
    self[*i] = (yield e)
  end
  #HACK: if there's a single element in a non-dense matrix, it won't iterate and
  #won't change the default value; this ensures that it does get changed.
  unless iterated then
    self.each_with_indices do |e, *i|
      self[*i] = (yield e)
    end
  end
end

static VALUE nm_map_stored(VALUE nmatrix) {
  NM_CONSERVATIVE(nm_register_value(&nmatrix));
  VALUE to_return = Qnil;

  switch(NM_STYPE(nmatrix)) {
  case nm::YALE_STORE:
    to_return = nm_yale_map_stored(nmatrix);
    break;
  case nm::DENSE_STORE:
    to_return = nm_dense_map(nmatrix);
    break;
  case nm::LIST_STORE:
    to_return = nm_list_map_stored(nmatrix, Qnil);
    break;
  default:
    NM_CONSERVATIVE(nm_unregister_value(&nmatrix));
    rb_raise(nm_eDataTypeError, "Not a proper storage type");
  }

  NM_CONSERVATIVE(nm_unregister_value(&nmatrix));
  return to_return;
}

# File lib/nmatrix/math.rb, line 712
def max(dimen=0)
  inject_rank(dimen) do |max, sub_mat|
    if max.is_a? NMatrix then
      max * (max >= sub_mat).cast(self.stype, self.dtype) + ((max)*0.0 + (max < sub_mat).cast(self.stype, self.dtype)) * sub_mat
    else
      max >= sub_mat ? max : sub_mat
    end
  end
end

# File lib/nmatrix/math.rb, line 659
def mean(dimen=0)
  reduce_dtype = nil
  if integer_dtype? then
    reduce_dtype = :float64
  end
  inject_rank(dimen, 0.0, reduce_dtype) do |mean, sub_mat|
    mean + sub_mat
  end / shape[dimen]
end

# File lib/nmatrix/math.rb, line 693
def min(dimen=0)
  inject_rank(dimen) do |min, sub_mat|
    if min.is_a? NMatrix then
      min * (min <= sub_mat).cast(self.stype, self.dtype) + ((min)*0.0 + (min > sub_mat).cast(self.stype, self.dtype)) * sub_mat
    else
      min <= sub_mat ? min : sub_mat
    end
  end
end

# File lib/nmatrix/math.rb, line 824
def nrm2 incx=1, n=nil
  return method_missing(:nrm2, incx, n) unless vector?
  NMatrix::BLAS::nrm2(self, incx, self.size / incx)
end

# File lib/nmatrix/nmatrix.rb, line 429
def nvector?
  self.effective_dim < self.dim
end

# File lib/nmatrix/nmatrix.rb, line 374
def object_dtype?
  dtype == :object
end

static VALUE nm_offset(VALUE self) {
  NM_CONSERVATIVE(nm_register_value(&self));
  STORAGE* s   = NM_STORAGE(self);

  // Copy elements into a VALUE array and then use those to create a Ruby array with rb_ary_new4.
  VALUE* offset = NM_ALLOCA_N(VALUE, s->dim);
  nm_register_values(offset, s->dim);
  for (size_t index = 0; index < s->dim; ++index)
    offset[index] = INT2FIX(s->offset[index]);

  nm_unregister_values(offset, s->dim);
  NM_CONSERVATIVE(nm_unregister_value(&self));
  return rb_ary_new4(s->dim, offset);
}

# File lib/nmatrix/atlas.rb, line 207
def potrf!(which)
  raise(StorageTypeError, "ATLAS functions only work on dense matrices") unless self.dense?
  raise(ShapeError, "Cholesky decomposition only valid for square matrices") unless self.dim == 2 && self.shape[0] == self.shape[1]

  NMatrix::LAPACK::clapack_potrf(:row, which, self.shape[0], self, self.shape[1])
end

# File lib/nmatrix/math.rb, line 188
def potrf_lower!
  potrf! :lower
end

# File lib/nmatrix/math.rb, line 184
def potrf_upper!
  potrf! :upper
end

# File lib/nmatrix/math.rb, line 542
def pow n
  raise ShapeError, "Only works with 2D square matrices." if 
    shape[0] != shape[1] or shape.size != 2
  raise TypeError, "Only works with integer powers" unless n.is_a?(Integer)

  sequence = (integer_dtype? ? self.cast(dtype: :int64) : self).clone
  product  = NMatrix.eye shape[0], dtype: sequence.dtype, stype: sequence.stype 

  if n == 0
    return NMatrix.eye(shape, dtype: dtype, stype: stype)
  elsif n == 1
    return sequence
  elsif n < 0
    n = n.abs
    sequence.invert!
    product = NMatrix.eye shape[0], dtype: sequence.dtype, stype: sequence.stype
  end

  # Decompose n to reduce the number of multiplications.
  while n > 0
    product = product.dot(sequence) if n % 2 == 1
    n = n / 2
    sequence = sequence.dot(sequence)
  end

  product
end

# File lib/nmatrix/homogeneous.rb, line 159
def quaternion
  raise(ShapeError, "Expected square matrix") if self.shape[0] != self.shape[1]
  raise(ShapeError, "Expected 3x3 rotation (or 4x4 homogeneous) matrix") if self.shape[0] > 4 || self.shape[0] < 3

  q = NMatrix.new([4], dtype: self.dtype == :float32 ? :float32: :float64)
  rotation_trace = self[0,0] + self[1,1] + self[2,2]
  if rotation_trace >= 0
    self_w = self.shape[0] == 4 ? self[3,3] : 1.0
    root_of_homogeneous_trace = Math.sqrt(rotation_trace + self_w)
    q[0] = root_of_homogeneous_trace * 0.5
    s = 0.5 / root_of_homogeneous_trace
    q[1] = (self[2,1] - self[1,2]) * s
    q[2] = (self[0,2] - self[2,0]) * s
    q[3] = (self[1,0] - self[0,1]) * s
  else
    h = 0
    h = 1 if self[1,1] > self[0,0]
    h = 2 if self[2,2] > self[h,h]

    case_macro = Proc.new do |i,j,k,ii,jj,kk|
      qq = NMatrix.new([4], dtype: :float64)
      self_w = self.shape[0] == 4 ? self[3,3] : 1.0
      s = Math.sqrt( (self[ii,ii] - (self[jj,jj] + self[kk,kk])) + self_w)
      qq[i] = s*0.5
      s = 0.5 / s
      qq[j] = (self[ii,jj] + self[jj,ii]) * s
      qq[k] = (self[kk,ii] + self[ii,kk]) * s
      qq[0] = (self[kk,jj] - self[jj,kk]) * s
      qq
    end

    case h
    when 0
      q = case_macro.call(1,2,3, 0,1,2)
    when 1
      q = case_macro.call(2,3,1, 1,2,0)
    when 2
      q = case_macro.call(3,1,2, 2,0,1)
    end

    self_w = self.shape[0] == 4 ? self[3,3] : 1.0
    if self_w != 1
      s = 1.0 / Math.sqrt(self_w)
      q[0] *= s
      q[1] *= s
      q[2] *= s
      q[3] *= s
    end
  end

  q
end

# File lib/nmatrix/nmatrix.rb, line 481
def rank(shape_idx, rank_idx, meth = :copy)

  if shape_idx > (self.dim-1)
    raise(RangeError, "#rank call was out of bounds")
  end

  params = Array.new(self.dim)
  params.each.with_index do |v,d|
    params[d] = d == shape_idx ? rank_idx : 0...self.shape[d]
  end

  meth == :reference ? self[*params] : self.slice(*params)
end

# File lib/nmatrix/nmatrix.rb, line 1008
def repeat(count, axis)
  raise(ArgumentError, 'Matrix should be repeated at least 2 times.') if count < 2
  new_shape = shape
  new_shape[axis] *= count
  new_matrix = NMatrix.new(new_shape)
  slice = new_shape.map { |axis_size| 0...axis_size }
  start = 0
  count.times do
    slice[axis] = start...(start += shape[axis])
    new_matrix[*slice] = self
  end
  new_matrix
end

# File lib/nmatrix/nmatrix.rb, line 550
def reshape new_shape,*shapes
  if new_shape.is_a?Fixnum
    newer_shape =  [new_shape]+shapes
  else  # new_shape is an Array
    newer_shape = new_shape
  end
  t = reshape_clone_structure(newer_shape)
  left_params  = [:*]*newer_shape.size
  right_params = [:*]*self.shape.size
  t[*left_params] = self[*right_params]
  t
end

# File lib/nmatrix/nmatrix.rb, line 575
def reshape! new_shape,*shapes
  if self.is_ref?
    raise(ArgumentError, "This operation cannot be performed on reference slices")
  else
    if new_shape.is_a?Fixnum
      shape =  [new_shape]+shapes
    else  # new_shape is an Array
      shape = new_shape
    end
    self.reshape_bang(shape)
  end
end

static VALUE nm_unary_round(int argc, VALUE* argv, VALUE self) {
  NM_CONSERVATIVE(nm_register_values(argv, argc));
  const int default_precision = 0;
  NMATRIX* left;
  UnwrapNMatrix(self, left);
  std::string sym;

  switch(left->stype) {
  case nm::DENSE_STORE:
    sym = "__dense_unary_round__";
    break;
  case nm::YALE_STORE:
    sym = "__yale_unary_round__";
    break;
  case nm::LIST_STORE:
    sym = "__list_unary_round__";
    break;
  }
  NM_CONSERVATIVE(nm_unregister_values(argv, argc));
  if (argc > 0) { //supplied precision
    return rb_funcall(self, rb_intern(sym.c_str()), 1, argv[0]);
  }
  return rb_funcall(self, rb_intern(sym.c_str()), 1, nm::RubyObject(default_precision).rval);
}

# File lib/nmatrix/nmatrix.rb, line 533
def row(row_number, get_by = :copy)
  rank(0, row_number, get_by)
end

# File lib/nmatrix/nmatrix.rb, line 256
def rows
  shape[0]
end

static VALUE nm_shape(VALUE self) {
  NM_CONSERVATIVE(nm_register_value(&self));
  STORAGE* s   = NM_STORAGE(self);

  // Copy elements into a VALUE array and then use those to create a Ruby array with rb_ary_new4.
  VALUE* shape = NM_ALLOCA_N(VALUE, s->dim);
  nm_register_values(shape, s->dim);
  for (size_t index = 0; index < s->dim; ++index)
    shape[index] = INT2FIX(s->shape[index]);

  nm_unregister_values(shape, s->dim);
  NM_CONSERVATIVE(nm_unregister_value(&self));
  return rb_ary_new4(s->dim, shape);
}

# File lib/nmatrix/nmatrix.rb, line 886
def shuffle(*args)
  method_missing(:shuffle!, *args) if self.effective_dim > 1
  t = self.clone
  t.shuffle!(*args)
end

# File lib/nmatrix/nmatrix.rb, line 868
def shuffle!(*args)
  method_missing(:shuffle!, *args) if self.effective_dim > 1
  ary = self.to_flat_a
  ary.shuffle!(*args)
  ary.each.with_index { |v,idx| self[idx] = v }
  self
end

# File lib/nmatrix/nmatrix.rb, line 413
def size
  NMatrix.size(self.shape)
end

static VALUE nm_mget(int argc, VALUE* argv, VALUE self) {
  static void* (*ttable[nm::NUM_STYPES])(const STORAGE*, SLICE*) = {
    nm_dense_storage_get,
    nm_list_storage_get,
    nm_yale_storage_get
  };
  nm::stype_t stype = NM_STYPE(self);
  return nm_xslice(argc, argv, ttable[stype], nm_delete, self);
}

# File lib/nmatrix/lapacke.rb, line 198
def solve b
  raise(ShapeError, "Must be called on square matrix") unless self.dim == 2 && self.shape[0] == self.shape[1]
  raise(ShapeError, "number of rows of b must equal number of cols of self") if 
    self.shape[1] != b.shape[0]
  raise ArgumentError, "only works with dense matrices" if self.stype != :dense
  raise ArgumentError, "only works for non-integer, non-object dtypes" if 
    integer_dtype? or object_dtype? or b.integer_dtype? or b.object_dtype?

  x     = b.clone
  clone = self.clone
  n = self.shape[0]
  ipiv = NMatrix::LAPACK.lapacke_getrf(:row, n, n, clone, n)
  NMatrix::LAPACK.lapacke_getrs(:row, :no_transpose, n, b.shape[1], clone, n, ipiv, x, b.shape[1])
  x
end

# File lib/nmatrix/nmatrix.rb, line 899
def sorted_indices
  return method_missing(:sorted_indices) unless vector?
  ary = self.to_flat_array
  ary.each_index.sort_by { |i| ary[i] }  # from: http://stackoverflow.com/a/17841159/170300
end

# File lib/nmatrix/math.rb, line 757
def std(dimen=0)
  variance(dimen).sqrt
end

static VALUE nm_stype(VALUE self) {
  NM_CONSERVATIVE(nm_register_value(&self));
  VALUE stype = ID2SYM(rb_intern(STYPE_NAMES[NM_STYPE(self)]));
  NM_CONSERVATIVE(nm_unregister_value(&self));
  return stype;
}

# File lib/nmatrix/math.rb, line 677
def sum(dimen=0)
  inject_rank(dimen, 0.0) do |sum, sub_mat|
    sum + sub_mat
  end
end

static VALUE nm_supershape(VALUE self) {

  STORAGE* s   = NM_STORAGE(self);
  if (s->src == s) {
    return nm_shape(self); // easy case (not a slice)
  }
  else s = s->src;

  NM_CONSERVATIVE(nm_register_value(&self));

  VALUE* shape = NM_ALLOCA_N(VALUE, s->dim);
  nm_register_values(shape, s->dim);
  for (size_t index = 0; index < s->dim; ++index)
    shape[index] = INT2FIX(s->shape[index]);

  nm_unregister_values(shape, s->dim);
  NM_CONSERVATIVE(nm_unregister_value(&self));
  return rb_ary_new4(s->dim, shape);
}

static VALUE nm_symmetric(VALUE self) {
  return is_symmetric(self, false);
}

# File lib/nmatrix/nmatrix.rb, line 451
def to_a(dimen=nil)
  if self.dim == 2

    return self.to_flat_a if self.shape[0] == 1

    ary = []
    begin
      self.each_row do |row|
        ary << row.to_flat_a
      end
    #rescue NotImplementedError # Oops. Try copying instead
    #  self.each_row(:copy) do |row|
    #    ary << row.to_a.flatten
    #  end
    end
    ary
  else
    to_a_rec(0)
  end
end

# File lib/nmatrix/nmatrix.rb, line 388
def to_f
  raise IndexError, 'to_f only valid for matrices with a single element' unless shape.all? { |e| e == 1 }
  self[*Array.new(shape.size, 0)]
end

# File lib/nmatrix/nmatrix.rb, line 400
def to_flat_array
  ary = Array.new(self.size)
  self.each.with_index { |v,i| ary[i] = v }
  ary
end

# File lib/nmatrix/nmatrix.rb, line 306
def to_hash
  if stype == :yale
    h = {}
    each_stored_with_indices do |val,i,j|
      next if val == 0 # Don't bother storing the diagonal zero values -- only non-zeros.
      if h.has_key?(i)
        h[i][j] = val
      else
        h[i] = {j => val}
      end
    end
    h
  else # dense and list should use a C internal function.
    # FIXME: Write a C internal to_h function.
    m = stype == :dense ? self.cast(:list, self.dtype) : self
    m.__list_to_hash__
  end
end

# File lib/nmatrix/math.rb, line 640
def trace
  raise(ShapeError, "Expected square matrix") unless self.shape[0] == self.shape[1] && self.dim == 2

  (0...self.shape[0]).inject(0) do |total,i|
    total + self[i,i]
  end
end

# File lib/nmatrix/nmatrix.rb, line 602
def transpose(permute = nil)
  if self.dim == 1
    return self.clone
  elsif self.dim == 2
    new_shape = [self.shape[1], self.shape[0]]
  elsif permute.nil?
    raise(ArgumentError, "need permutation array of size #{self.dim}")
  elsif permute.sort.uniq != (0...self.dim).to_a
    raise(ArgumentError, "invalid permutation array")
  else
    # Figure out the new shape based on the permutation given as an argument.
    new_shape = permute.map { |p| self.shape[p] }
  end

  if self.dim > 2 # FIXME: For dense, several of these are basically equivalent to reshape.

    # Make the new data structure.
    t = self.reshape_clone_structure(new_shape)

    self.each_stored_with_indices do |v,*indices|
      p_indices = permute.map { |p| indices[p] }
      t[*p_indices] = v
    end
    t
  elsif self.list? # TODO: Need a C list transposition algorithm.
    # Make the new data structure.
    t = self.reshape_clone_structure(new_shape)

    self.each_column.with_index do |col,j|
      t[j,:*] = col.to_flat_array
    end
    t
  else
    # Call C versions of Yale and List transpose, which do their own copies
    self.clone_transpose
  end
end

# File lib/nmatrix/nmatrix.rb, line 745
def upper_triangle(k = 0)
  raise(NotImplementedError, "only implemented for 2D matrices") if self.shape.size > 2

  t = self.clone_structure
  (0...self.shape[0]).each do |i|
    if i - k < 0
      t[i, :*] = self[i, :*]
    else
      t[i, 0...(i-k)]             = 0
      t[i, (i-k)...self.shape[1]] = self[i, (i-k)...self.shape[1]]
    end
  end
  t
end

# File lib/nmatrix/nmatrix.rb, line 774
def upper_triangle!(k = 0)
  raise(NotImplementedError, "only implemented for 2D matrices") if self.shape.size > 2

  (0...self.shape[0]).each do |i|
    if i - k >= 0
      self[i, 0...(i-k)] = 0
    end
  end
  self
end

# File lib/nmatrix/math.rb, line 734
def variance(dimen=0)
  reduce_dtype = nil
  if integer_dtype? then
    reduce_dtype = :float64
  end
  m = mean(dimen)
  inject_rank(dimen, 0.0, reduce_dtype) do |var, sub_mat|
    var + (m - sub_mat)*(m - sub_mat)/(shape[dimen]-1)
  end
end

# File lib/nmatrix/nmatrix.rb, line 722
def vconcat(*matrices)
  concat(*matrices, :row)
end

# File lib/nmatrix/nmatrix.rb, line 439
def vector?
  self.effective_dim == 1
end

static VALUE nm_write(int argc, VALUE* argv, VALUE self) {
  using std::ofstream;

  if (argc < 1 || argc > 2) {
    rb_raise(rb_eArgError, "Expected one or two arguments");
  }

  NM_CONSERVATIVE(nm_register_values(argv, argc));
  NM_CONSERVATIVE(nm_register_value(&self));

  VALUE file = argv[0],
        symm = argc == 1 ? Qnil : argv[1];

  NMATRIX* nmatrix;
  UnwrapNMatrix( self, nmatrix );

  nm::symm_t symm_ = interpret_symm(symm);

  if (nmatrix->storage->dtype == nm::RUBYOBJ) {
    NM_CONSERVATIVE(nm_unregister_values(argv, argc));
    NM_CONSERVATIVE(nm_unregister_value(&self));
    rb_raise(rb_eNotImpError, "Ruby Object writing is not implemented yet");
  }

  // Get the dtype, stype, itype, and symm and ensure they're the correct number of bytes.
  uint8_t st = static_cast<uint8_t>(nmatrix->stype),
          dt = static_cast<uint8_t>(nmatrix->storage->dtype),
          sm = static_cast<uint8_t>(symm_);
  uint16_t dim = nmatrix->storage->dim;

  //FIXME: Cast the matrix to the smallest possible index type. Write that in the place of IType.

  // Check arguments before starting to write.
  if (nmatrix->stype == nm::LIST_STORE) {
    NM_CONSERVATIVE(nm_unregister_values(argv, argc));
    NM_CONSERVATIVE(nm_unregister_value(&self));
    rb_raise(nm_eStorageTypeError, "cannot save list matrix; cast to yale or dense first");
  }
  if (symm_ != nm::NONSYMM) {
    NM_CONSERVATIVE(nm_unregister_values(argv, argc));
    NM_CONSERVATIVE(nm_unregister_value(&self));

    if (dim != 2) rb_raise(rb_eArgError, "symmetry/triangularity not defined for a non-2D matrix");
    if (nmatrix->storage->shape[0] != nmatrix->storage->shape[1])
      rb_raise(rb_eArgError, "symmetry/triangularity not defined for a non-square matrix");
    if (symm_ == nm::HERM &&
          dt != static_cast<uint8_t>(nm::COMPLEX64) && dt != static_cast<uint8_t>(nm::COMPLEX128) && dt != static_cast<uint8_t>(nm::RUBYOBJ))
      rb_raise(rb_eArgError, "cannot save a non-complex matrix as hermitian");
  }

  ofstream f(RSTRING_PTR(file), std::ios::out | std::ios::binary);

  // Get the NMatrix version information.
  uint16_t major, minor, release, null16 = 0;
  get_version_info(major, minor, release);

  // WRITE FIRST 64-BIT BLOCK
  f.write(reinterpret_cast<const char*>(&major),   sizeof(uint16_t));
  f.write(reinterpret_cast<const char*>(&minor),   sizeof(uint16_t));
  f.write(reinterpret_cast<const char*>(&release), sizeof(uint16_t));
  f.write(reinterpret_cast<const char*>(&null16),  sizeof(uint16_t));

  uint8_t ZERO = 0;
  // WRITE SECOND 64-BIT BLOCK
  f.write(reinterpret_cast<const char*>(&dt), sizeof(uint8_t));
  f.write(reinterpret_cast<const char*>(&st), sizeof(uint8_t));
  f.write(reinterpret_cast<const char*>(&ZERO),sizeof(uint8_t));
  f.write(reinterpret_cast<const char*>(&sm), sizeof(uint8_t));
  f.write(reinterpret_cast<const char*>(&null16), sizeof(uint16_t));
  f.write(reinterpret_cast<const char*>(&dim), sizeof(uint16_t));

  // Write shape (in 64-bit blocks)
  write_padded_shape(f, nmatrix->storage->dim, nmatrix->storage->shape);

  if (nmatrix->stype == nm::DENSE_STORE) {
    write_padded_dense_elements(f, reinterpret_cast<DENSE_STORAGE*>(nmatrix->storage), symm_, nmatrix->storage->dtype);
  } else if (nmatrix->stype == nm::YALE_STORE) {
    YALE_STORAGE* s = reinterpret_cast<YALE_STORAGE*>(nmatrix->storage);
    uint32_t ndnz   = s->ndnz,
             length = nm_yale_storage_get_size(s);
    f.write(reinterpret_cast<const char*>(&ndnz),   sizeof(uint32_t));
    f.write(reinterpret_cast<const char*>(&length), sizeof(uint32_t));

    write_padded_yale_elements(f, s, length, symm_, s->dtype);
  }

  f.close();

  NM_CONSERVATIVE(nm_unregister_values(argv, argc));
  NM_CONSERVATIVE(nm_unregister_value(&self));

  return Qtrue;
}

# File lib/nmatrix/shortcuts.rb, line 47
def yale?;  return stype == :yale; end

# File lib/nmatrix/math.rb, line 921
def __dense_unary_round__(precision)
  if self.complex_dtype?
    self.__dense_map__ { |l| Complex(l.real.round(precision), l.imag.round(precision)) }
                                  .cast(stype, dtype)
  else
    self.__dense_map__ { |l| l.round(precision) }.cast(stype, dtype)
  end
end

# File lib/nmatrix/math.rb, line 903
def __list_unary_round__(precision)
  if self.complex_dtype?
    self.__list_map_stored__(nil) { |l| Complex(l.real.round(precision), l.imag.round(precision)) }
                                  .cast(stype, dtype)
  else
    self.__list_map_stored__(nil) { |l| l.round(precision) }.cast(stype, dtype)
  end
end

# File lib/nmatrix/math.rb, line 912
def __yale_unary_round__(precision)
  if self.complex_dtype?
    self.__yale_map_stored__ { |l| Complex(l.real.round(precision), l.imag.round(precision)) }
                                  .cast(stype, dtype)
  else
    self.__yale_map_stored__ { |l| l.round(precision) }.cast(stype, dtype)
  end
end

# File lib/nmatrix/math.rb, line 931
def dtype_for_floor_or_ceil
  if self.integer_dtype? or [:complex64, :complex128, :object].include?(self.dtype)
    return_dtype = dtype
  elsif [:float32, :float64].include?(self.dtype)
    return_dtype = :int64
  end

  return_dtype
end