/*
  tensor_source.c
  Ruby/GSL: Ruby extension library for GSL (GNU Scientific Library)
    (C) Copyright 2004 by Yoshiki Tsunesada

  Ruby/GSL is free software: you can redistribute it and/or modify it
  under the terms of the GNU General Public License.
  This library is distributed in the hope that it will be useful, but
  WITHOUT ANY WARRANTY.
*/

/*
  The gsl_tensor package is developed by J. Burguet, and
  distributed separately as an add-on package.
 */

#ifdef HAVE_GSL_TENSOR_GSL_TENSOR_H

#include "rb_gsl_config.h"
#include "rb_gsl_tensor.h"
#include "rb_gsl_common.h"

#ifdef BASE_DOUBLE
VALUE cgsl_tensor, cgsl_tensor_int;
VALUE cgsl_tensor_view, cgsl_tensor_int_view;
#define NUMCONV(x) NUM2DBL(x)
#define C_TO_VALUE(x) rb_float_new(x)
#define CHECK_TEN(x) CHECK_TENSOR(x)
#define TEN_P(x) TENSOR_P(x)
#define VEC_P(x) VECTOR_P(x)
#define MAT_P(x) MATRIX_P(x)
#else defined(BASE_INT)
#define NUMCONV(x) FIX2INT(x)
#define C_TO_VALUE(x) INT2FIX(x)
#define CHECK_TEN(x) CHECK_TENSOR_INT(x)
#define TEN_P(x) TENSOR_INT_P(x)
#define VEC_P(x) VECTOR_INT_P(x)
#define MAT_P(x) MATRIX_INT_P(x)
#endif

GSL_TYPE(rbgsl_tensor)* FUNCTION(rbgsl_tensor,alloc)(const unsigned int rank, 
						     const size_t dimension)
{
  GSL_TYPE(rbgsl_tensor) *t;
  t = ALLOC(GSL_TYPE(rbgsl_tensor));
  t->tensor = FUNCTION(gsl_tensor,alloc)(rank, dimension);
  if (rank == 0) 
    t->indices = gsl_permutation_alloc(1);
  else
    t->indices = gsl_permutation_alloc(rank);
  return t;
}

GSL_TYPE(rbgsl_tensor)* FUNCTION(rbgsl_tensor,calloc)(const unsigned int rank, 
						      const size_t dimension)
{
  GSL_TYPE(rbgsl_tensor) *t;
  t = ALLOC(GSL_TYPE(rbgsl_tensor));
  t->tensor = FUNCTION(gsl_tensor,calloc)(rank, dimension);
  if (rank == 0) 
    t->indices = gsl_permutation_alloc(1);
  else
    t->indices = gsl_permutation_alloc(rank);
  return t;
}

GSL_TYPE(rbgsl_tensor)* FUNCTION(rbgsl_tensor,copy)(const GSL_TYPE(rbgsl_tensor) *t)
{
  GSL_TYPE(rbgsl_tensor) *tnew;
  tnew = ALLOC(GSL_TYPE(rbgsl_tensor));
  if (t->tensor->rank == 0) 
    tnew->indices = gsl_permutation_alloc(1);
  else
    tnew->indices = gsl_permutation_alloc(t->tensor->rank);
  tnew->tensor = FUNCTION(gsl_tensor,copy)(t->tensor);
  return tnew;
}

void FUNCTION(rbgsl_tensor,free)(GSL_TYPE(rbgsl_tensor) *t)
{
  gsl_permutation_free(t->indices);
  FUNCTION(gsl_tensor,free)(t->tensor);
  free((GSL_TYPE(rbgsl_tensor) *) t);
}

void FUNCTION(rbgsl_tensor,free2)(GSL_TYPE(rbgsl_tensor) *t)
{
  gsl_permutation_free(t->indices);
  free((GSL_TYPE(gsl_tensor)*) t->tensor);
  free((GSL_TYPE(rbgsl_tensor) *) t);
}

/* singleton methods */
static VALUE FUNCTION(rb_gsl_tensor,new)(int argc, VALUE *argv, VALUE klass)
{
  unsigned int rank;
  size_t dim;
  GSL_TYPE(rbgsl_tensor) *t = NULL;
  switch (argc) {
  case 2:
    rank = FIX2UINT(argv[0]);
    dim = (size_t) FIX2UINT(argv[1]);
    t = FUNCTION(rbgsl_tensor,alloc)(rank, dim);
    break;
  default:
    rb_raise(rb_eArgError, "wrong number of arguments (%d for 2, rank and dimension)",
	     argc);
    break;
  }
  return Data_Wrap_Struct(GSL_TYPE(cgsl_tensor), 0, FUNCTION(rbgsl_tensor,free), t);
}

static VALUE FUNCTION(rb_gsl_tensor,calloc)(VALUE klass, VALUE r, VALUE s)
{
  unsigned int rank;
  size_t dim;
  GSL_TYPE(rbgsl_tensor) *t = NULL;
  rank = FIX2UINT(r);
  dim = (size_t) FIX2UINT(s);
  t = FUNCTION(rbgsl_tensor,calloc)(rank, dim);
  return Data_Wrap_Struct(GSL_TYPE(cgsl_tensor), 0, FUNCTION(rbgsl_tensor,free), t);
}

static VALUE FUNCTION(rb_gsl_tensor,copy_singleton)(VALUE klass, VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t = NULL, *tnew = NULL;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  tnew = FUNCTION(rbgsl_tensor,copy(t));
  return Data_Wrap_Struct(GSL_TYPE(cgsl_tensor), 0, FUNCTION(rbgsl_tensor,free), tnew);
}

static VALUE FUNCTION(rb_gsl_tensor,memcpy_singleton)(VALUE klass, VALUE a, VALUE b)
{
  GSL_TYPE(rbgsl_tensor) *dst, *src;
  CHECK_TEN(b);
  Data_Get_Struct(a, GSL_TYPE(rbgsl_tensor), dst);
  Data_Get_Struct(b, GSL_TYPE(rbgsl_tensor), src);
  return INT2FIX(FUNCTION(gsl_tensor,memcpy)(dst->tensor, src->tensor));
}

static VALUE FUNCTION(rb_gsl_tensor,swap_singleton)(VALUE klass, VALUE a, VALUE b)
{
  GSL_TYPE(rbgsl_tensor) *t1, *t2;
  CHECK_TEN(b);
  Data_Get_Struct(a, GSL_TYPE(rbgsl_tensor), t1);
  Data_Get_Struct(b, GSL_TYPE(rbgsl_tensor), t2);
  return INT2FIX(FUNCTION(gsl_tensor,swap)(t1->tensor, t2->tensor));
}

/*****/
static VALUE FUNCTION(rb_gsl_tensor,copy)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t = NULL, *tnew = NULL;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  tnew = FUNCTION(rbgsl_tensor,copy)(t);
  return Data_Wrap_Struct(GSL_TYPE(cgsl_tensor), 0, FUNCTION(rbgsl_tensor,free), tnew);
}

static VALUE FUNCTION(rb_gsl_tensor,set_zero)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t = NULL;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  FUNCTION(gsl_tensor,set_zero)(t->tensor);
  return obj;
}

static VALUE FUNCTION(rb_gsl_tensor,set_all)(VALUE obj, VALUE xx)
{
  GSL_TYPE(rbgsl_tensor) *t = NULL;
  BASE x;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  x = NUMCONV(xx);
  FUNCTION(gsl_tensor,set_all)(t->tensor, x);
  return obj;
}

static void rb_gsl_tensor_get_indices_array(gsl_tensor_indices *v, VALUE ary);
static void rbgsl_tensor_get_indices(int argc, VALUE *argv, gsl_tensor_indices *indices, 
				     size_t *n);
#ifdef BASE_DOUBLE
static void rb_gsl_tensor_get_indices_array(gsl_tensor_indices *v, VALUE ary)
{
  size_t i, nn;
  nn = (size_t) GSL_MIN_INT((int) v->size, (int) RARRAY(ary)->len);
  for (i = 0; i < nn; i++)
    v->data[i] = FIX2UINT(rb_ary_entry(ary, i));
}

static void rbgsl_tensor_get_indices(int argc, VALUE *argv, 
				     gsl_tensor_indices *indices, size_t *n)
{
  size_t i;
  for (i = 0; i < indices->size; i++) indices->data[i] = 0;
  switch (argc) {
  case 1:
    switch (TYPE(argv[0])) {
    case T_ARRAY:
      *n = (size_t) GSL_MIN_INT((int) indices->size, (int) RARRAY(argv[0])->len);
      rb_gsl_tensor_get_indices_array(indices, argv[0]);
      break;
    case T_FIXNUM:
      *n = 1;
      indices->data[0] = FIX2INT(argv[0]);
      break;
    default:
      rb_raise(rb_eTypeError, "wrong argument type %s (Array expected)",
	       rb_class2name(CLASS_OF(argv[0])));
      break;
    }
    break;
  default:
    *n = (size_t) GSL_MIN_INT(argc, (int) indices->size);
    for (i = 0; i < *n; i++) {
      CHECK_FIXNUM(argv[i]);
      indices->data[i] = FIX2INT(argv[i]);
    }
    break;
  }
}
#endif

size_t FUNCTION(gsl_tensor,position)(const size_t * indices,
				     const GSL_TYPE(gsl_tensor) * t);
static VALUE FUNCTION(rb_gsl_tensor,position)(int argc, VALUE *argv, VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t = NULL;
  size_t n, position;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);  
  rbgsl_tensor_get_indices(argc, argv, t->indices, &n);  
  position = (size_t) FUNCTION(gsl_tensor,position)(t->indices->data,t->tensor);
  return INT2FIX(position);
}
static VALUE FUNCTION(rb_gsl_tensor,subtensor)(int argc, VALUE *argv, VALUE obj);
static VALUE FUNCTION(rb_gsl_tensor,get)(int argc, VALUE *argv, VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t = NULL;
  BASE x;
  size_t n;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);  
  rbgsl_tensor_get_indices(argc, argv, t->indices, &n);  
  if (n < t->tensor->rank) {
    return FUNCTION(rb_gsl_tensor,subtensor)(argc, argv, obj);
  } else {
    x = FUNCTION(gsl_tensor,get)(t->tensor, t->indices->data);
    return C_TO_VALUE(x);
  }
  return Qnil;
}

static VALUE FUNCTION(rb_gsl_tensor,set)(int argc, VALUE *argv, VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t = NULL;
  size_t n;
  BASE x;
  if (argc < 2) 
    rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 2)", argc);
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);  
  rbgsl_tensor_get_indices(argc-1, argv, t->indices, &n);
  x = NUMCONV(argv[argc-1]);
  FUNCTION(gsl_tensor,set)(t->tensor, t->indices->data, x);
  return obj;
}

static VALUE FUNCTION(rb_gsl_tensor,fread)(VALUE obj, VALUE io)
{
  GSL_TYPE(rbgsl_tensor) *t = NULL;
  FILE *f = NULL;
  int status, flag = 0;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  f = rb_gsl_open_readfile(io, &flag);
  status = FUNCTION(gsl_tensor,fread)(f, t->tensor);
  if (flag == 1) fclose(f);
  return INT2FIX(status);
}

static VALUE FUNCTION(rb_gsl_tensor,fwrite)(VALUE obj, VALUE io)
{
  GSL_TYPE(rbgsl_tensor) *t = NULL;
  FILE *f = NULL;
  int status, flag = 0;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  f = rb_gsl_open_writefile(io, &flag);
  status = FUNCTION(gsl_tensor,fwrite)(f, t->tensor);
  if (flag == 1) fclose(f);
  return INT2FIX(status);
}

static VALUE FUNCTION(rb_gsl_tensor,fprintf)(int argc, VALUE *argv, VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *h = NULL;
  FILE *fp = NULL;
  int status, flag = 0;
  if (argc != 1 && argc != 2) 
    rb_raise(rb_eArgError, "wrong number of arguments (%d for 1 or 2)", argc);
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), h);
  fp = rb_gsl_open_writefile(argv[0], &flag);
  switch (argc) {
  case 2:
    if (TYPE(argv[1]) == T_STRING)
      status = FUNCTION(gsl_tensor,fprintf)(fp, h->tensor, STR2CSTR(argv[1]));
    else
      rb_raise(rb_eTypeError, "argv 2 String expected");
    break;
  default:
    status = FUNCTION(gsl_tensor,fprintf)(fp, h->tensor, OUT_FORMAT);
    break;
  }
  if (flag == 1) fclose(fp);
  return INT2FIX(status);
}

static VALUE FUNCTION(rb_gsl_tensor,printf)(int argc, VALUE *argv, VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *h = NULL;
  int status;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), h);
  if (argc == 1) {
    if (TYPE(argv[0]) != T_STRING) 
      rb_raise(rb_eTypeError, "String expected");
    else
      status = FUNCTION(gsl_tensor,fprintf)(stdout, h->tensor, STR2CSTR(argv[0]));
  } else {
    status = FUNCTION(gsl_tensor,fprintf)(stdout, h->tensor, OUT_FORMAT);
  }
  return INT2FIX(status);
}

static VALUE FUNCTION(rb_gsl_tensor,fscanf)(VALUE obj, VALUE io)
{
  GSL_TYPE(rbgsl_tensor) *h = NULL;
  FILE *fp = NULL;
  int status, flag = 0;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), h);
  fp = rb_gsl_open_readfile(io, &flag);
  status = FUNCTION(gsl_tensor,fscanf)(fp, h->tensor);
  if (flag == 1) fclose(fp);
  return INT2FIX(status);
}

static VALUE FUNCTION(rb_gsl_tensor,swap_indices)(VALUE obj, VALUE ii, VALUE jj)
{
  GSL_TYPE(rbgsl_tensor) *t, *tnew;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  tnew = ALLOC(GSL_TYPE(rbgsl_tensor));
  if (t->tensor->rank == 0) 
    tnew->indices = gsl_permutation_alloc(1);
  else
    tnew->indices = gsl_permutation_alloc(t->tensor->rank);
  tnew->tensor = FUNCTION(gsl_tensor,swap_indices)(t->tensor, FIX2INT(ii), FIX2INT(jj));
  return Data_Wrap_Struct(GSL_TYPE(cgsl_tensor), 0, FUNCTION(rbgsl_tensor,free), tnew);
}

static VALUE FUNCTION(rb_gsl_tensor,max)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  return C_TO_VALUE(FUNCTION(gsl_tensor,max)(t->tensor));
}

static VALUE FUNCTION(rb_gsl_tensor,min)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  return C_TO_VALUE(FUNCTION(gsl_tensor,min)(t->tensor));
}

static VALUE FUNCTION(rb_gsl_tensor,minmax)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  BASE min, max;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  FUNCTION(gsl_tensor,minmax)(t->tensor, &min, &max);
  return rb_ary_new3(2, C_TO_VALUE(min), C_TO_VALUE(max));
}

static VALUE FUNCTION(rb_gsl_tensor,max_index)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  FUNCTION(gsl_tensor,max_index)(t->tensor, t->indices->data);
  return Data_Wrap_Struct(cgsl_index, 0, NULL, t->indices);
}

static VALUE FUNCTION(rb_gsl_tensor,min_index)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  FUNCTION(gsl_tensor,min_index)(t->tensor, t->indices->data);
  return Data_Wrap_Struct(cgsl_index, 0, NULL, t->indices);
}

static VALUE FUNCTION(rb_gsl_tensor,minmax_index)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  gsl_permutation *min, *max;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  if (t->tensor->rank == 0) {
    min = gsl_permutation_alloc(1);
    max = gsl_permutation_alloc(1);
  } else {
    min = gsl_permutation_alloc(t->tensor->rank);
    max = gsl_permutation_alloc(t->tensor->rank);
  }
  FUNCTION(gsl_tensor,minmax_index)(t->tensor, min->data, max->data);
  return rb_ary_new3(2,
		     Data_Wrap_Struct(cgsl_index, 0, gsl_permutation_free, min), 
		     Data_Wrap_Struct(cgsl_index, 0, gsl_permutation_free, max));
}

static VALUE FUNCTION(rb_gsl_tensor,isnull)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  return INT2FIX(FUNCTION(gsl_tensor,isnull)(t->tensor));
}

static VALUE FUNCTION(rb_gsl_tensor,isnull2)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  int status;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  status = FUNCTION(gsl_tensor,isnull)(t->tensor);
  if (status) return Qtrue;
  else return Qfalse;
}

static VALUE FUNCTION(rb_gsl_tensor,oper)(VALUE obj, VALUE bb,
					   int flag)
{
  GSL_TYPE(rbgsl_tensor) *a, *b, *anew;
  BASE x;
  int (*f)(GSL_TYPE(gsl_tensor)*, const GSL_TYPE(gsl_tensor)*);
  int (*f2)(GSL_TYPE(gsl_tensor)*, const double);
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), a);
  anew = FUNCTION(rbgsl_tensor,copy)(a);
  if (TEN_P(bb)) {
    Data_Get_Struct(bb, GSL_TYPE(rbgsl_tensor), b);
    switch (flag) {
    case GSL_TENSOR_ADD: f = FUNCTION(&gsl_tensor,add); break;
    case GSL_TENSOR_SUB: f = FUNCTION(&gsl_tensor,sub); break;
    case GSL_TENSOR_MUL_ELEMENTS: f = FUNCTION(&gsl_tensor,mul_elements); break;
    case GSL_TENSOR_DIV_ELEMENTS: f = FUNCTION(&gsl_tensor,div_elements); break;
    default: rb_raise(rb_eRuntimeError, "unknown operation"); break;
    } 
    (*f)(anew->tensor, b->tensor);
  } else {
    switch (flag) {
    case GSL_TENSOR_ADD:
    case GSL_TENSOR_ADD_CONSTANT: 
      x = NUMCONV(bb);
      f2 = FUNCTION(&gsl_tensor,add_constant); 
      break;
    case GSL_TENSOR_SUB:
      x = -NUMCONV(bb);
      f2 = FUNCTION(&gsl_tensor,add_constant); 
      break;
    case GSL_TENSOR_ADD_DIAGONAL: f2 = FUNCTION(&gsl_tensor,add_diagonal); break;
    case GSL_TENSOR_MUL_ELEMENTS:
    case GSL_TENSOR_SCALE:
      x = NUMCONV(bb); 
      f2 = FUNCTION(&gsl_tensor,scale); 
      break;
    case GSL_TENSOR_DIV_ELEMENTS: 
      x = 1.0/NUMCONV(bb); 
      f2 = FUNCTION(&gsl_tensor,scale); 
      break;
    default: rb_raise(rb_eRuntimeError, "unknown operation"); break;
    }
    (*f2)(anew->tensor, x);
  }
  return Data_Wrap_Struct(GSL_TYPE(cgsl_tensor), 0, FUNCTION(rbgsl_tensor,free), anew);
}


static VALUE FUNCTION(rb_gsl_tensor,oper_bang)(VALUE obj, VALUE bb, int flag)
{
  GSL_TYPE(rbgsl_tensor) *a, *b;
  BASE x;
  int (*f)(GSL_TYPE(gsl_tensor)*, const GSL_TYPE(gsl_tensor)*);
  int (*f2)(GSL_TYPE(gsl_tensor)*, const double);
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), a);
  if (TEN_P(bb)) {
    Data_Get_Struct(bb, GSL_TYPE(rbgsl_tensor), b);
    switch (flag) {
    case GSL_TENSOR_ADD: f = FUNCTION(&gsl_tensor,add); break;
    case GSL_TENSOR_SUB: f = FUNCTION(&gsl_tensor,sub); break;
    case GSL_TENSOR_MUL_ELEMENTS: f = FUNCTION(&gsl_tensor,mul_elements); break;
    case GSL_TENSOR_DIV_ELEMENTS: f = FUNCTION(&gsl_tensor,div_elements); break;
    default: rb_raise(rb_eRuntimeError, "unknown operation"); break;
    } 
    (*f)(a->tensor, b->tensor);
  } else {
    switch (flag) {
    case GSL_TENSOR_ADD:
    case GSL_TENSOR_ADD_CONSTANT: 
      x = NUMCONV(bb);
      f2 = FUNCTION(&gsl_tensor,add_constant); 
      break;
    case GSL_TENSOR_SUB:
      x = -NUMCONV(bb);
      f2 = FUNCTION(&gsl_tensor,add_constant); 
      break;
    case GSL_TENSOR_ADD_DIAGONAL: f2 = FUNCTION(&gsl_tensor,add_diagonal); break;
    case GSL_TENSOR_MUL_ELEMENTS:
    case GSL_TENSOR_SCALE:
      x = NUMCONV(bb); 
      f2 = FUNCTION(&gsl_tensor,scale); 
      break;
    case GSL_TENSOR_DIV_ELEMENTS: 
      x = 1.0/NUMCONV(bb); 
      f2 = FUNCTION(&gsl_tensor,scale); 
      break;
    default: rb_raise(rb_eRuntimeError, "unknown operation"); break;
    }
    (*f2)(a->tensor, x);
  }
  return obj;
}

static VALUE FUNCTION(rb_gsl_tensor,add)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper)(obj, bb, GSL_TENSOR_ADD);
}

static VALUE FUNCTION(rb_gsl_tensor,sub)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper)(obj, bb, GSL_TENSOR_SUB);
}

static VALUE FUNCTION(rb_gsl_tensor,mul_elements)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper)(obj, bb, GSL_TENSOR_MUL_ELEMENTS);
}

static VALUE FUNCTION(rb_gsl_tensor,div_elements)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper)(obj, bb, GSL_TENSOR_DIV_ELEMENTS);
}

static VALUE FUNCTION(rb_gsl_tensor,add_constant)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper)(obj, bb, GSL_TENSOR_ADD_CONSTANT);
}

static VALUE FUNCTION(rb_gsl_tensor,add_diagonal)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper)(obj, bb, GSL_TENSOR_ADD_DIAGONAL);
}

static VALUE FUNCTION(rb_gsl_tensor,scale)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper)(obj, bb, GSL_TENSOR_SCALE);
}

/***/
static VALUE FUNCTION(rb_gsl_tensor,add_bang)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper_bang)(obj, bb, GSL_TENSOR_ADD);
}

static VALUE FUNCTION(rb_gsl_tensor,sub_bang)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper_bang)(obj, bb, GSL_TENSOR_SUB);
}

static VALUE FUNCTION(rb_gsl_tensor,mul_elements_bang)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper_bang)(obj, bb, GSL_TENSOR_MUL_ELEMENTS);
}

static VALUE FUNCTION(rb_gsl_tensor,div_elements_bang)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper_bang)(obj, bb, GSL_TENSOR_DIV_ELEMENTS);
}

static VALUE FUNCTION(rb_gsl_tensor,add_constant_bang)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper_bang)(obj, bb, GSL_TENSOR_ADD_CONSTANT);
}

static VALUE FUNCTION(rb_gsl_tensor,add_diagonal_bang)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper_bang)(obj, bb, GSL_TENSOR_ADD_DIAGONAL);
}

static VALUE FUNCTION(rb_gsl_tensor,scale_bang)(VALUE obj, VALUE bb)
{
  return FUNCTION(rb_gsl_tensor,oper_bang)(obj, bb, GSL_TENSOR_SCALE);
}

/*****/
static VALUE FUNCTION(rb_gsl_tensor,product_singleton)(VALUE obj, VALUE aa, VALUE bb)
{
  GSL_TYPE(rbgsl_tensor) *a, *b, *c;
  CHECK_TEN(aa);
  CHECK_TEN(bb);
  Data_Get_Struct(aa, GSL_TYPE(rbgsl_tensor), a);
  switch (TYPE(bb)) {
  case T_FIXNUM:
  case T_BIGNUM:
  case T_FLOAT:
    return FUNCTION(rb_gsl_tensor,mul_elements)(aa, bb);
    break;
  default:
    Data_Get_Struct(bb, GSL_TYPE(rbgsl_tensor), b);
    c = ALLOC(GSL_TYPE(rbgsl_tensor));
    c->tensor = FUNCTION(gsl_tensor,product(a->tensor, b->tensor));
    if (c->tensor->rank == 0)
      c->indices = gsl_permutation_alloc(1);
    else
      c->indices = gsl_permutation_alloc(c->tensor->rank);
    return Data_Wrap_Struct(GSL_TYPE(cgsl_tensor), 0, FUNCTION(rbgsl_tensor,free), c);
    break;
  }
}

static VALUE FUNCTION(rb_gsl_tensor,product)(VALUE obj, VALUE bb)
{
  GSL_TYPE(rbgsl_tensor) *a, *b, *c;
  switch (TYPE(bb)) {
  case T_FIXNUM:
  case T_BIGNUM:
  case T_FLOAT:
    return FUNCTION(rb_gsl_tensor,mul_elements)(obj, bb);
    break;
  default:
    CHECK_TEN(bb);
    Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), a);
    Data_Get_Struct(bb, GSL_TYPE(rbgsl_tensor), b);
    c = ALLOC(GSL_TYPE(rbgsl_tensor));
    c->tensor = FUNCTION(gsl_tensor,product(a->tensor, b->tensor));
    if (c->tensor->rank == 0)
      c->indices = gsl_permutation_alloc(1);
    else
      c->indices = gsl_permutation_alloc(c->tensor->rank);
    return Data_Wrap_Struct(GSL_TYPE(cgsl_tensor), 0, FUNCTION(rbgsl_tensor,free), c);
    break;
  }
}

static VALUE FUNCTION(rb_gsl_tensor,contract)(VALUE obj, VALUE ii, VALUE jj)
{
  GSL_TYPE(rbgsl_tensor) *t, *tnew;
  size_t rank;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  tnew = ALLOC(GSL_TYPE(rbgsl_tensor));
  tnew->tensor = FUNCTION(gsl_tensor,contract)(t->tensor, FIX2INT(ii), FIX2INT(jj));
  if (tnew->tensor->rank == 0) rank = 1;
  else rank = tnew->tensor->rank;
  tnew->indices = gsl_permutation_alloc(rank);
  return Data_Wrap_Struct(GSL_TYPE(cgsl_tensor), 0, FUNCTION(rbgsl_tensor,free), tnew);
}

static VALUE FUNCTION(rb_gsl_tensor,size)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  return INT2FIX(t->tensor->size);
}

static VALUE FUNCTION(rb_gsl_tensor,rank)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  return INT2FIX(t->tensor->rank);
}

static VALUE FUNCTION(rb_gsl_tensor,dimension)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  return INT2FIX(t->tensor->dimension);
}

static VALUE FUNCTION(rb_gsl_tensor,data)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  QUALIFIED_VIEW(gsl_vector,view) *v = NULL;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  v = FUNCTION(rb_gsl_make_vector,view)(t->tensor->data, t->tensor->size, 1);
  return Data_Wrap_Struct(QUALIFIED_VIEW(cgsl_vector,view), 0, free, v);
}

static VALUE FUNCTION(rb_gsl_tensor,2matrix)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  GSL_TYPE(gsl_matrix) *m = NULL;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  m = (GSL_TYPE(gsl_matrix)*) FUNCTION(gsl_tensor,2matrix)(t->tensor);
  return Data_Wrap_Struct(QUALIFIED_VIEW(cgsl_matrix,view), 0, FUNCTION(gsl_matrix,free), m);
}

static VALUE FUNCTION(rb_gsl_tensor,2vector)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  GSL_TYPE(gsl_vector) *v = NULL;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  v = (GSL_TYPE(gsl_vector) *) FUNCTION(gsl_tensor,2vector)(t->tensor);
  return Data_Wrap_Struct(QUALIFIED_VIEW(cgsl_vector,view), 0, FUNCTION(gsl_vector,free), v);
}

static VALUE FUNCTION(rb_gsl_tensor,to_v)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  GSL_TYPE(gsl_vector) *v;

  v = FUNCTION(gsl_vector,alloc)(t->tensor->size);
  memcpy(v->data, t->tensor->data, sizeof(BASE)*v->size);
  return Data_Wrap_Struct(GSL_TYPE(cgsl_vector), 0, FUNCTION(gsl_vector,free), v);
}

/*
  Creates a subtensor slicing the existing tensor.
  NOTE: no new data region is malloced.
  t: Tensor
  rank: rank of the tensor created
*/
GSL_TYPE(gsl_tensor) FUNCTION(gsl_tensor,subtensor)(const GSL_TYPE(gsl_tensor) *t,
						     const unsigned int rank,
						     size_t *indices)
{
  GSL_TYPE(gsl_tensor) tnew;
  size_t position;
  tnew.rank = rank;
  tnew.dimension = t->dimension;
  tnew.size =  quick_pow(t->dimension, rank);
  position = FUNCTION(gsl_tensor,position)(indices, t);
  if (position >= t->size)
    rb_raise(rb_eRangeError, "wrong indices given");
  tnew.data = t->data + position;
  return tnew;
}

static VALUE FUNCTION(rb_gsl_tensor,subtensor)(int argc, VALUE *argv, VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t, *tnew;
  unsigned int rank;
  size_t n;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  /* n: number of indices given */
  rbgsl_tensor_get_indices(argc, argv, t->indices, &n);  
  rank = t->tensor->rank - n; 
  tnew = ALLOC(GSL_TYPE(rbgsl_tensor));
  tnew->tensor = (GSL_TYPE(gsl_tensor)*) malloc(sizeof(GSL_TYPE(gsl_tensor)));
  *(tnew->tensor) = FUNCTION(gsl_tensor,subtensor)(t->tensor, rank, t->indices->data);
  if (rank == 0) 
    tnew->indices = gsl_permutation_alloc(1);
  else
    tnew->indices = gsl_permutation_alloc(rank);
  return Data_Wrap_Struct(QUALIFIED_VIEW(cgsl_tensor,view), 0, FUNCTION(rbgsl_tensor,free2), tnew);	       
}

#ifdef BASE_DOUBLE
#define SHOW_ELM 6
#define PRINTF_FORMAT "%4.3e "
#else
#define SHOW_ELM 15
#define PRINTF_FORMAT "%d "
#endif
static VALUE FUNCTION(rb_gsl_tensor,to_s)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  QUALIFIED_VIEW(gsl_matrix,view) matrix;
  QUALIFIED_VIEW(gsl_vector,view) vector;
  GSL_TYPE(gsl_matrix) *m;
  GSL_TYPE(gsl_vector) *v;
  char buf[16];
  size_t i, j;
  VALUE str;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  str = rb_str_new2("[ ");
  switch (t->tensor->rank) {
  case 2:
    matrix.matrix.data = t->tensor->data;
    matrix.matrix.size1 = t->tensor->dimension;
    matrix.matrix.size2 = t->tensor->dimension;
    matrix.matrix.tda = t->tensor->dimension;
    matrix.matrix.block = 0;
    matrix.matrix.owner = 0;
    m = &(matrix.matrix);
    for (i = 0; i < m->size1; i++) {
      if (i != 0) {
	strcpy(buf, "  ");
	rb_str_cat(str, buf, strlen(buf));
      }
      for (j = 0; j < m->size2; j++) {
	sprintf(buf, PRINTF_FORMAT, FUNCTION(gsl_matrix,get)(m, i, j));
	rb_str_cat(str, buf, strlen(buf));
	if (j == SHOW_ELM) {
	  strcpy(buf, "... ");
	  rb_str_cat(str, buf, strlen(buf));
	  break;
	}
      }
      if (i == 6) {
	strcpy(buf, "\n  ... ]");
	rb_str_cat(str, buf, strlen(buf));
	break;
      }
      if (i == m->size1 - 1) {
	strcpy(buf, "]");
	rb_str_cat(str, buf, strlen(buf));
      } else {
	strcpy(buf, "\n");
	rb_str_cat(str, buf, strlen(buf));
      }
    }
    return str;
    break;
  default:
    vector.vector.data = t->tensor->data;
    vector.vector.stride = 1;
    vector.vector.size = t->tensor->size;
    vector.vector.owner = 0;
    vector.vector.block = 0;
    v = &(vector.vector);
    sprintf(buf,  PRINTF_FORMAT, FUNCTION(gsl_vector,get)(v, 0));
    rb_str_cat(str, buf, strlen(buf));
    for (i = 1; i < v->size; i++) {
      sprintf(buf,  PRINTF_FORMAT, FUNCTION(gsl_vector,get)(v, i));
      rb_str_cat(str, buf, strlen(buf));
      if (i == SHOW_ELM && i != v->size-1) {
        strcpy(buf, "... ");
        rb_str_cat(str, buf, strlen(buf));
        break;
      }
    }
    sprintf(buf, "]");
    rb_str_cat(str, buf, strlen(buf));
    return str;
    break;
  }
}
#undef SHOW_ELM
#undef PRINTF_FORMAT

static VALUE FUNCTION(rb_gsl_tensor,inspect)(VALUE obj)
{
  VALUE str;
  char buf[64];
  sprintf(buf, "%s\n", rb_class2name(CLASS_OF(obj)));
  str = rb_str_new2(buf);
  return rb_str_concat(str, FUNCTION(rb_gsl_tensor,to_s)(obj));
}

VALUE FUNCTION(rb_gsl_tensor,equal)(int argc, VALUE *argv, VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *a, *b;
  GSL_TYPE(gsl_matrix) *m;
  GSL_TYPE(gsl_vector) *v;
  VALUE other;
  double eps = 1e-10;
  size_t i;
  switch (argc) {
  case 2:
    other = argv[0];
    eps = NUM2DBL(argv[1]);
    break;
  case 1:
    other = argv[0];
    break;
  default:
    rb_raise(rb_eArgError, "wrong number of arguments (%d for 1 or 2)", argc);
  }
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), a);
  if (TEN_P(other)) {
    Data_Get_Struct(other, GSL_TYPE(rbgsl_tensor), b);
    if (a->tensor->rank != b->tensor->rank) return Qfalse;
    if (a->tensor->dimension != b->tensor->dimension) return Qfalse;
    if (a->tensor->size != b->tensor->size) return Qfalse;
    for (i = 0; i < a->tensor->size; i++)
      if (fabs(a->tensor->data[i]-b->tensor->data[i]) > eps) 
	return Qfalse;
    return Qtrue;
  } else if (MAT_P(other)) {
    if (a->tensor->rank != 2) return Qfalse;
    Data_Get_Struct(other, GSL_TYPE(gsl_matrix), m);
    if (a->tensor->dimension != m->size1 || a->tensor->dimension != m->size2)
      return Qfalse;
    for (i = 0; i < a->tensor->size; i++)
      if (fabs(a->tensor->data[i]-m->data[i]) > eps) 
	return Qfalse;
    return Qtrue;
  } else if (VEC_P(other)) {
    Data_Get_Struct(other, GSL_TYPE(gsl_vector), v);
    if (a->tensor->size != v->size) return Qfalse;
    for (i = 0; i < a->tensor->size; i++)
      if (fabs(a->tensor->data[i]-v->data[i]) > eps) 
	return Qfalse;
    return Qtrue;
  } else {
    rb_raise(rb_eTypeError, "wrong argument type %s (Tensor, Matrix or Vector expected)", 
	     rb_class2name(CLASS_OF(other)));
  }
}

static VALUE FUNCTION(rb_gsl_tensor,uplus)(VALUE obj)
{
  return obj;
}

static VALUE FUNCTION(rb_gsl_tensor,uminus)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t, *tnew;
  size_t i;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  tnew = FUNCTION(rbgsl_tensor,copy)(t);
  for (i = 0; i < tnew->tensor->size; i++)
    tnew->tensor->data[i] *= -1;
  return Data_Wrap_Struct(GSL_TYPE(cgsl_tensor), 0, FUNCTION(rbgsl_tensor,free), tnew);
}

static VALUE FUNCTION(rb_gsl_tensor,coerce)(VALUE obj, VALUE other)
{
  GSL_TYPE(rbgsl_tensor) *t, *tnew;
  VALUE tt;
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  switch (TYPE(other)) {
  case T_FLOAT:
  case T_FIXNUM:
  case T_BIGNUM:
    tnew = FUNCTION(rbgsl_tensor,alloc)(t->tensor->rank, t->tensor->dimension);
    FUNCTION(gsl_tensor,set_all)(tnew->tensor, NUMCONV(other));
    tt = Data_Wrap_Struct(GSL_TYPE(cgsl_tensor), 0, FUNCTION(rbgsl_tensor,free), tnew);
    return rb_ary_new3(2, tt, obj);
    break;
  default:
    rb_raise(rb_eRuntimeError, "undefined operation with %s", 
	     rb_class2name(CLASS_OF(other)));
    break;
  }
}

static VALUE FUNCTION(rb_gsl_tensor,info)(VALUE obj)
{
  GSL_TYPE(rbgsl_tensor) *t;
  char buf[256];
  Data_Get_Struct(obj, GSL_TYPE(rbgsl_tensor), t);
  sprintf(buf, "Class:      %s\n", rb_class2name(CLASS_OF(obj)));
  sprintf(buf, "%sSuperClass: %s\n", buf, rb_class2name(RCLASS(CLASS_OF(obj))->super));
  sprintf(buf, "%sRank:       %d\n", buf, (int) t->tensor->rank);
  sprintf(buf, "%sDimension:  %d\n", buf, (int) t->tensor->dimension);
  sprintf(buf, "%sSize:       %d\n", buf, (int) t->tensor->size);
  return rb_str_new2(buf);
}

void FUNCTION(Init_gsl_tensor,init)(VALUE module)
{
#ifdef BASE_DOUBLE
  cgsl_tensor = rb_define_class_under(module, "Tensor", cGSL_Object);
  cgsl_tensor_int = rb_define_class_under(cgsl_tensor, "Int", cGSL_Object);
  cgsl_tensor_view = rb_define_class_under(cgsl_tensor, "View", cgsl_tensor);
  cgsl_tensor_int_view = rb_define_class_under(cgsl_tensor_int, "View",
					       cgsl_tensor_int);
  /*
  cgsl_index = rb_define_class_under(cgsl_tensor, "Index", 
  cgsl_permutation);*/
#endif

  rb_define_singleton_method(GSL_TYPE(cgsl_tensor), "new",
			     FUNCTION(rb_gsl_tensor,new), -1);
  rb_define_singleton_method(GSL_TYPE(cgsl_tensor), "[]",
			     FUNCTION(rb_gsl_tensor,new), -1);
  rb_define_singleton_method(GSL_TYPE(cgsl_tensor), "alloc",
			     FUNCTION(rb_gsl_tensor,new), -1);
  rb_define_singleton_method(GSL_TYPE(cgsl_tensor), "calloc",
			     FUNCTION(rb_gsl_tensor,calloc), 2);
  rb_define_singleton_method(GSL_TYPE(cgsl_tensor), "copy",
			     FUNCTION(rb_gsl_tensor,copy_singleton), 1);
  rb_define_singleton_method(GSL_TYPE(cgsl_tensor), "memcpy",
			     FUNCTION(rb_gsl_tensor,memcpy_singleton), 2);
  rb_define_singleton_method(GSL_TYPE(cgsl_tensor), "swap",
			     FUNCTION(rb_gsl_tensor,swap_singleton), 2);

  /*****/

  rb_define_method(GSL_TYPE(cgsl_tensor), "copy",
			     FUNCTION(rb_gsl_tensor,copy), 0);
  rb_define_alias(GSL_TYPE(cgsl_tensor), "clone", "copy");
  rb_define_alias(GSL_TYPE(cgsl_tensor), "duplicate", "copy");
  rb_define_method(GSL_TYPE(cgsl_tensor), "set_zero",
			     FUNCTION(rb_gsl_tensor,set_zero), 0);
  rb_define_method(GSL_TYPE(cgsl_tensor), "set_all",
			     FUNCTION(rb_gsl_tensor,set_all), 1);

  rb_define_method(GSL_TYPE(cgsl_tensor), "position",
		   FUNCTION(rb_gsl_tensor,position), -1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "get",
			     FUNCTION(rb_gsl_tensor,get), -1);
  rb_define_alias(GSL_TYPE(cgsl_tensor), "[]", "get");
  rb_define_method(GSL_TYPE(cgsl_tensor), "set",
			     FUNCTION(rb_gsl_tensor,set), -1);
  rb_define_alias(GSL_TYPE(cgsl_tensor), "[]=", "set");

  rb_define_method(GSL_TYPE(cgsl_tensor), "fread",
			     FUNCTION(rb_gsl_tensor,fread), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "fwrite",
			     FUNCTION(rb_gsl_tensor,fwrite), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "fprintf",
			     FUNCTION(rb_gsl_tensor,fprintf), -1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "printf",
			     FUNCTION(rb_gsl_tensor,printf), -1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "fscanf",
			     FUNCTION(rb_gsl_tensor,fscanf), 1);

  rb_define_method(GSL_TYPE(cgsl_tensor), "swap_indices",
			     FUNCTION(rb_gsl_tensor,swap_indices), 2);

  rb_define_method(GSL_TYPE(cgsl_tensor), "max",
			     FUNCTION(rb_gsl_tensor,max), 0);
  rb_define_method(GSL_TYPE(cgsl_tensor), "min",
			     FUNCTION(rb_gsl_tensor,min), 0);
  rb_define_method(GSL_TYPE(cgsl_tensor), "minmax",
			     FUNCTION(rb_gsl_tensor,minmax), 0);
  rb_define_method(GSL_TYPE(cgsl_tensor), "max_index",
			     FUNCTION(rb_gsl_tensor,max_index), 0);
  rb_define_method(GSL_TYPE(cgsl_tensor), "min_index",
			     FUNCTION(rb_gsl_tensor,min_index), 0);
  rb_define_method(GSL_TYPE(cgsl_tensor), "minmax_index",
			     FUNCTION(rb_gsl_tensor,minmax_index), 0);

  rb_define_method(GSL_TYPE(cgsl_tensor), "isnull",
			     FUNCTION(rb_gsl_tensor,isnull), 0);
  rb_define_method(GSL_TYPE(cgsl_tensor), "isnull?",
			     FUNCTION(rb_gsl_tensor,isnull2), 0);

  rb_define_method(GSL_TYPE(cgsl_tensor), "add",
			     FUNCTION(rb_gsl_tensor,add), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "sub",
			     FUNCTION(rb_gsl_tensor,sub), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "mul_elements",
			     FUNCTION(rb_gsl_tensor,mul_elements), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "div_elements",
			     FUNCTION(rb_gsl_tensor,div_elements), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "add_constant",
			     FUNCTION(rb_gsl_tensor,add_constant), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "add_diagonal",
			     FUNCTION(rb_gsl_tensor,add_diagonal), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "scale",
			     FUNCTION(rb_gsl_tensor,scale), 1);
  rb_define_singleton_method(GSL_TYPE(cgsl_tensor), "product",
			     FUNCTION(rb_gsl_tensor,product_singleton), 2);
  rb_define_method(GSL_TYPE(cgsl_tensor), "product",
			     FUNCTION(rb_gsl_tensor,product), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "contract",
			     FUNCTION(rb_gsl_tensor,contract), 2);

  rb_define_alias(GSL_TYPE(cgsl_tensor), "+", "add");
  rb_define_alias(GSL_TYPE(cgsl_tensor), "-", "sub");
  /*  rb_define_alias(GSL_TYPE(cgsl_tensor), "*", "mul_elements");*/
  rb_define_alias(GSL_TYPE(cgsl_tensor), "/", "div_elements");
  rb_define_alias(GSL_TYPE(cgsl_tensor), "*", "product");

  rb_define_method(GSL_TYPE(cgsl_tensor), "add!",
			     FUNCTION(rb_gsl_tensor,add_bang), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "sub!",
			     FUNCTION(rb_gsl_tensor,sub_bang), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "mul_elements!",
			     FUNCTION(rb_gsl_tensor,mul_elements_bang), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "div_elements!",
			     FUNCTION(rb_gsl_tensor,div_elements_bang), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "add_constant!",
			     FUNCTION(rb_gsl_tensor,add_constant_bang), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "add_diagonal!",
			     FUNCTION(rb_gsl_tensor,add_diagonal_bang), 1);
  rb_define_method(GSL_TYPE(cgsl_tensor), "scale!",
			     FUNCTION(rb_gsl_tensor,scale_bang), 1);

  rb_define_method(GSL_TYPE(cgsl_tensor), "+@",
			     FUNCTION(rb_gsl_tensor,uplus), 0);
  rb_define_method(GSL_TYPE(cgsl_tensor), "-@",
			     FUNCTION(rb_gsl_tensor,uminus), 0);

  rb_define_method(GSL_TYPE(cgsl_tensor), "size",
			     FUNCTION(rb_gsl_tensor,size), 0);
  rb_define_method(GSL_TYPE(cgsl_tensor), "rank",
			     FUNCTION(rb_gsl_tensor,rank), 0);
  rb_define_method(GSL_TYPE(cgsl_tensor), "dimension",
			     FUNCTION(rb_gsl_tensor,dimension), 0);
  rb_define_alias(GSL_TYPE(cgsl_tensor), "dim", "dimension");
  rb_define_method(GSL_TYPE(cgsl_tensor), "data",
			     FUNCTION(rb_gsl_tensor,data), 0);

  rb_define_method(GSL_TYPE(cgsl_tensor), "to_v",
			     FUNCTION(rb_gsl_tensor,to_v), 0);
  rb_define_alias(GSL_TYPE(cgsl_tensor), "to_gv", "to_v");

  rb_define_method(GSL_TYPE(cgsl_tensor), "to_vector",
			     FUNCTION(rb_gsl_tensor,2vector), 0);
  rb_define_method(GSL_TYPE(cgsl_tensor), "to_matrix",
			     FUNCTION(rb_gsl_tensor,2matrix), 0);

  rb_define_method(GSL_TYPE(cgsl_tensor), "subtensor",
			     FUNCTION(rb_gsl_tensor,subtensor), -1);
  rb_define_alias(GSL_TYPE(cgsl_tensor), "view", "subtensor");

  rb_define_method(GSL_TYPE(cgsl_tensor), "to_s",
			     FUNCTION(rb_gsl_tensor,to_s), 0);
  rb_define_method(GSL_TYPE(cgsl_tensor), "inspect",
			     FUNCTION(rb_gsl_tensor,inspect), 0);

  rb_define_method(GSL_TYPE(cgsl_tensor), "equal?",
			     FUNCTION(rb_gsl_tensor,equal), -1);
  rb_define_alias(GSL_TYPE(cgsl_tensor), "==", "equal?");

  rb_define_method(GSL_TYPE(cgsl_tensor), "coerce",
			     FUNCTION(rb_gsl_tensor,coerce), 1);

  rb_define_method(GSL_TYPE(cgsl_tensor), "info", 
		   FUNCTION(rb_gsl_tensor,info), 0);
}

#undef NUMCONV
#undef C_TO_VALUE
#undef CHECK_TEN
#undef TEN_P
#undef VEC_P
#undef MAT_P

#endif