linalg.c

/*
  linalg.c
  Ruby/GSL: Ruby extension library for GSL (GNU Scientific Library)
    (C) Copyright 2001-2006 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.
*/

#include <gsl/gsl_math.h>
#include "include/rb_gsl_array.h"
#include "include/rb_gsl_common.h"
#include "include/rb_gsl_linalg.h"

static VALUE cgsl_matrix_LU;
static VALUE cgsl_matrix_QR;
static VALUE cgsl_matrix_QRPT;
static VALUE cgsl_vector_tau;
static VALUE cgsl_matrix_Q;
static VALUE cgsl_matrix_R;

static VALUE cgsl_matrix_LQ;
static VALUE cgsl_matrix_PTLQ;
static VALUE cgsl_matrix_L;

static VALUE cgsl_matrix_SV;
static VALUE cgsl_matrix_U;
static VALUE cgsl_matrix_V;
static VALUE cgsl_vector_S;
static VALUE cgsl_matrix_C;

enum {
  LINALG_DECOMP,
  LINALG_DECOMP_BANG,
};

#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_linalg_LU_decomp_narray(int argc, VALUE *argv, VALUE obj,
                                            int flag)
{
  struct NARRAY *na, *na2;
  VALUE m;
  gsl_matrix_view mv;
  gsl_permutation *p;
  int signum;

  if (argc != 1)
    rb_raise(rb_eArgError, "wrong number of arguments (%d for 1)", argc);
  GetNArray(argv[0], na);
  if (na->rank < 2) rb_raise(rb_eRuntimeError, "rank >= 2 required");
  if (na->shape[0] != na->shape[1])
    rb_raise(rb_eRuntimeError, "square matrix required");
  if (flag == LINALG_DECOMP) {
    m = na_make_object(NA_DFLOAT, 2, na->shape, CLASS_OF(argv[0]));
    GetNArray(m, na2);
    memcpy((double*)na2->ptr, (double*)na->ptr, sizeof(double)*na2->total);
    mv = gsl_matrix_view_array((double*)na2->ptr, na->shape[1], na->shape[0]);
  } else {
    mv = gsl_matrix_view_array((double*)na->ptr, na->shape[1], na->shape[0]);
  }
  p = gsl_permutation_alloc(mv.matrix.size1);
  gsl_linalg_LU_decomp(&mv.matrix, p, &signum);
  if (flag == LINALG_DECOMP) {
    return rb_ary_new3(3, m,
                       Data_Wrap_Struct(cgsl_permutation, 0, gsl_permutation_free, p),
                       INT2FIX(signum));
  } else {
    return rb_ary_new3(3, argv[0],
                       Data_Wrap_Struct(cgsl_permutation, 0, gsl_permutation_free, p),
                       INT2FIX(signum));
  }
}
#endif

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_LU_decomp_nmatrix(int argc, VALUE *argv, VALUE obj,
                                            int flag)
{
  NM_DENSE_STORAGE *input_nmatrix, *temp_nmatrix;
  VALUE m;
  gsl_matrix_view mv;
  gsl_permutation *p;
  int signum;

  if (argc != 1)
    rb_raise(rb_eArgError, "wrong number of arguments (%d for 1)", argc);

  input_nmatrix = NM_STORAGE_DENSE(argv[0]);
  if (input_nmatrix->shape[0] != input_nmatrix->shape[1])
    rb_raise(rb_eRuntimeError, "square matrix required");

  if (flag == LINALG_DECOMP) {
    m = rb_nmatrix_dense_create(FLOAT64, input_nmatrix->shape, 2, 
      input_nmatrix->elements, input_nmatrix->shape[0] * input_nmatrix->shape[1]);
    temp_nmatrix = NM_STORAGE_DENSE(m);
    mv = gsl_matrix_view_array((double*)temp_nmatrix->elements, 
      temp_nmatrix->shape[0], temp_nmatrix->shape[1]);
  } else {
    mv = gsl_matrix_view_array((double*)input_nmatrix->elements, 
      input_nmatrix->shape[0], input_nmatrix->shape[1]);
  }
  p = gsl_permutation_alloc(mv.matrix.size1);
  gsl_linalg_LU_decomp(&mv.matrix, p, &signum);
  if (flag == LINALG_DECOMP) {
    return rb_ary_new3(3, m,
                       Data_Wrap_Struct(cgsl_permutation, 0, gsl_permutation_free, p),
                       INT2FIX(signum));
  } else {
    return rb_ary_new3(3, argv[0],
                       Data_Wrap_Struct(cgsl_permutation, 0, gsl_permutation_free, p),
                       INT2FIX(signum));
  }
}
#endif

static VALUE rb_gsl_linalg_LU_decomposition(int argc, VALUE *argv, VALUE obj, int flag)
{
  gsl_matrix *mtmp = NULL, *m = NULL;
  gsl_permutation *p = NULL;
  int signum, itmp;
  size_t size;
  VALUE objp, objm, omatrix;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0])) {
      return rb_gsl_linalg_LU_decomp_narray(argc, argv, obj, flag);
    }
#endif

#ifdef HAVE_NMATRIX_H
    if (NM_IsNMatrix(argv[0])) {
      return rb_gsl_linalg_LU_decomp_nmatrix(argc, argv, obj, flag);
    }
#endif
    if (MATRIX_COMPLEX_P(argv[0]))
      return rb_gsl_linalg_complex_LU_decomp2(argc, argv, obj);
    omatrix = argv[0];
    itmp = 1;
    break;
  default:
    if (MATRIX_COMPLEX_P(obj))
      return rb_gsl_linalg_complex_LU_decomp2(argc, argv, obj);
    omatrix = obj;
    itmp = 0;
    break;
  }
  CHECK_MATRIX(omatrix);
  Data_Get_Struct(omatrix, gsl_matrix, mtmp);
  if (flag == LINALG_DECOMP_BANG) {
    m = mtmp;
    RBGSL_SET_CLASS(omatrix, cgsl_matrix_LU);
    objm = omatrix;
  } else {
    m = make_matrix_clone(mtmp);
    objm = Data_Wrap_Struct(cgsl_matrix_LU, 0, gsl_matrix_free, m);
  }
  size = m->size1;
  switch (argc-itmp) {
  case 0:
    p = gsl_permutation_alloc(size);
    gsl_linalg_LU_decomp(m, p, &signum);
    objp = Data_Wrap_Struct(cgsl_permutation, 0, gsl_permutation_free, p);
    if (flag == LINALG_DECOMP_BANG) return rb_ary_new3(2, objp, INT2FIX(signum));
    else return rb_ary_new3(3, objm, objp, INT2FIX(signum));
    break;
  case 1:
    CHECK_PERMUTATION(argv[itmp]);
    Data_Get_Struct(argv[itmp], gsl_permutation, p);
    gsl_linalg_LU_decomp(m, p, &signum);
    if (flag == LINALG_DECOMP_BANG) return INT2FIX(signum);
    else return rb_ary_new3(2, objm, INT2FIX(signum));
    break;
  default:
    rb_raise(rb_eArgError, "Usage: LU_decomp() or LU_decomp(permutation)");
    break;
  }
  return Qnil; /* never reach here */
}

static VALUE rb_gsl_linalg_LU_decomp(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_LU_decomposition(argc, argv, obj, LINALG_DECOMP);
}

static VALUE rb_gsl_linalg_LU_decomp_bang(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_LU_decomposition(argc, argv, obj, LINALG_DECOMP_BANG);
}

static gsl_matrix* get_matrix(VALUE obj, VALUE klass,int *flagm);
static gsl_permutation* get_permutation(VALUE obj,  size_t size, int *flagp);
static gsl_vector* get_vector2(VALUE obj,  int *flagv);

static gsl_matrix* get_matrix(VALUE obj, VALUE klass, int *flagm)
{
  gsl_matrix *mtmp = NULL, *m = NULL;
  if (CLASS_OF(obj) == klass) {
    Data_Get_Struct(obj, gsl_matrix, m);
    *flagm = 0;
#ifdef HAVE_NARRAY_H
  } else if (NA_IsNArray(obj)) {
    gsl_matrix_view mv;
    struct NARRAY *na;
    GetNArray(obj, na);
    mv = gsl_matrix_view_array((double*)na->ptr, na->shape[1], na->shape[0]);
    m = &mv.matrix;
    *flagm = -1;
#endif
  } else {
    CHECK_MATRIX(obj);
    Data_Get_Struct(obj, gsl_matrix, mtmp);
    m = make_matrix_clone(mtmp);
    *flagm = 1;
  }
  return m;
}

static gsl_permutation* get_permutation(VALUE obj, size_t size, int *flagp)
{
  gsl_permutation *p = NULL;
  if (CLASS_OF(obj) == cgsl_permutation) {
    Data_Get_Struct(obj, gsl_permutation, p);
    *flagp = 0;
  } else {
    p = gsl_permutation_alloc(size);
    *flagp = 1;
  }
  return p;
}

static gsl_vector* get_vector2(VALUE obj, int *flagv)
{
  gsl_vector *v = NULL;
  if (TYPE(obj) == T_ARRAY) {
    v = make_cvector_from_rarray(obj);
    *flagv = 1;
#ifdef HAVE_NARRAY_H
  } else if (NA_IsNArray(obj)) {
    gsl_vector_view vv;
    struct NARRAY *na;
    GetNArray(obj, na);
    vv = gsl_vector_view_array((double*) na->ptr, na->total);
    v = &vv.vector;
    *flagv = -1;
#endif
  } else {
    CHECK_VECTOR(obj);
    Data_Get_Struct(obj, gsl_vector, v);
    *flagv = 0;
  }
  return v;
}

#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_linalg_LU_solve_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *na, *b;
  VALUE ret;
  gsl_permutation *p;
  gsl_matrix_view mv;
  gsl_vector_view bv, xv;
  double *x;
  int shape[1];
  if (argc < 3)
    rb_raise(rb_eArgError,
             "wrong number of arguments %d(NArray, GSL::Permutation and NArray expected",
             argc);
  GetNArray(argv[0], na);
  mv = gsl_matrix_view_array((double*) na->ptr, na->shape[1], na->shape[0]);
  CHECK_PERMUTATION(argv[1]);
  Data_Get_Struct(argv[1], gsl_permutation, p);
  GetNArray(argv[2], b);
  bv = gsl_vector_view_array((double*) b->ptr, b->total);
  if (argc == 3) {
    shape[0] = b->total;
    ret = na_make_object(NA_DFLOAT, 1, shape, CLASS_OF(argv[0]));
  } else {
    ret = argv[3];
  }
  x = NA_PTR_TYPE(ret,double*);
  xv = gsl_vector_view_array(x, b->total);
  gsl_linalg_LU_solve(&mv.matrix, p, &bv.vector, &xv.vector);
  return ret;
}
#endif

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_LU_solve_nmatrix(int argc, VALUE *argv, VALUE obj)
{
  VALUE ret;
  NM_DENSE_STORAGE *input_nmatrix, *b;
  gsl_permutation *p;
  gsl_matrix_view mv;
  gsl_vector_view bv, xv;
  double *x;
  unsigned int shape[1];

  if (argc < 3) {
    rb_raise(rb_eArgError,
             "wrong number of arguments %d(NMatrix, GSL::Permutation and NMatrix expected)",
             argc);
  }
  input_nmatrix = NM_STORAGE_DENSE(argv[0]);
  mv = gsl_matrix_view_array((double*) input_nmatrix->elements, 
    input_nmatrix->shape[0], input_nmatrix->shape[1]);
  CHECK_PERMUTATION(argv[1]);
  Data_Get_Struct(argv[1], gsl_permutation, p);
  b = NM_STORAGE_DENSE(argv[2]);
  bv = gsl_vector_view_array((double*) b->elements, b->shape[0]);
  if (argc == 3) {
    shape[0] = b->shape[0];
    ret = rb_nmatrix_dense_create(FLOAT64, shape, 1, b->elements, shape[0]);
  } else {
    ret = argv[3];
  }
  x = (double*)NM_DENSE_ELEMENTS(ret);
  xv = gsl_vector_view_array(x, b->shape[0]);
  gsl_linalg_LU_solve(&mv.matrix, p, &bv.vector, &xv.vector);
  return ret;
}
#endif

VALUE rb_gsl_linalg_LU_solve(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *m = NULL;
  gsl_permutation *p = NULL;
  gsl_vector *b = NULL, *x = NULL;
  int signum, flagm = 0, flagp = 0, flagb = 0, flagx = 0, itmp;
  size_t size;
  VALUE bb;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc < 2 || argc > 4)
      rb_raise(rb_eArgError, "Usage: solve(m, b), solve(m, b, x), solve(lu, p, b), solve(lu, p, b, x)");
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0]))
      return rb_gsl_linalg_LU_solve_narray(argc, argv, obj);
#endif

#ifdef HAVE_NMATRIX_H
    if (NM_IsNMatrix(argv[0]))
      return rb_gsl_linalg_LU_solve_nmatrix(argc, argv, obj);
#endif
    m = get_matrix(argv[0], cgsl_matrix_LU, &flagm);
    itmp = 1;
    break;
  default:
    if (argc < 1 || argc > 3)
      rb_raise(rb_eArgError, "Usage: LU_solve(b), LU_solve(p, b), LU_solve(b, x), solve(p, b, x)");
    m = get_matrix(obj, cgsl_matrix_LU, &flagm);
    itmp = 0;
    break;
  }
  size = m->size1;

  p = get_permutation(argv[itmp], size, &flagp);
  if (flagp == 1 && flagm == 0) rb_raise(rb_eArgError, "permutation must be given");
  if (flagp == 0) itmp++;
  bb = argv[itmp];
  b = get_vector2(argv[itmp], &flagb);
  itmp++;

  if (itmp == argc) {
    x = gsl_vector_alloc(size);
    flagx = 1;
  } else {
    CHECK_VECTOR(argv[itmp]);
    Data_Get_Struct(argv[itmp], gsl_vector, x);
  }
  if (flagm == 1) gsl_linalg_LU_decomp(m, p, &signum);
  gsl_linalg_LU_solve(m, p, b, x);
  if (flagm == 1) gsl_matrix_free(m);
  if (flagp == 1) gsl_permutation_free(p);
  if (flagb == 1) gsl_vector_free(b);
  if (flagx == 1) return Data_Wrap_Struct(VECTOR_ROW_COL(bb), 0, gsl_vector_free, x);
  else return argv[argc-1];
}

#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_linalg_LU_svx_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *na, *b;
  gsl_permutation *p;
  gsl_matrix_view mv;
  gsl_vector_view bv;
  if (argc != 3)
    rb_raise(rb_eArgError,
             "wrong number of arguments %d(NArray, GSL::Permutation and NArray expected",
             argc);
  GetNArray(argv[0], na);
  mv = gsl_matrix_view_array((double*) na->ptr, na->shape[1], na->shape[0]);
  CHECK_PERMUTATION(argv[1]);
  Data_Get_Struct(argv[1], gsl_permutation, p);
  GetNArray(argv[2], b);
  bv = gsl_vector_view_array((double*) b->ptr, b->total);
  gsl_linalg_LU_svx(&mv.matrix, p, &bv.vector);
  return argv[2];
}
#endif

/* bb must be Vector, it is replaced by the root of the system */
static VALUE rb_gsl_linalg_LU_svx(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *m = NULL;
  gsl_permutation *p = NULL;
  gsl_vector *b = NULL;
  int signum, flagm = 0, flagp = 0, flagb = 0, itmp;
  size_t size;

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc < 2 || argc > 3)
      rb_raise(rb_eArgError, "Usage: svx(m, b), svx(lu, p, b)");
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0]))
      return rb_gsl_linalg_LU_svx_narray(argc, argv, obj);
#endif
    m = get_matrix(argv[0], cgsl_matrix_LU, &flagm);
    itmp = 1;
    break;
  default:
    if (argc < 1 || argc > 2)
      rb_raise(rb_eArgError, "Usage: LU_svx(b), LU_svx(p, b)");
    m = get_matrix(obj, cgsl_matrix_LU, &flagm);
    itmp = 0;
    break;
  }
  size = m->size1;
  p = get_permutation(argv[itmp], size, &flagp);
  if (flagp == 1 && flagm == 0) rb_raise(rb_eArgError, "permutation must be given");
  if (flagp == 0) itmp++;
  CHECK_VECTOR(argv[itmp]);
  b = get_vector2(argv[itmp], &flagb);
  if (flagm == 1) gsl_linalg_LU_decomp(m, p, &signum);
  gsl_linalg_LU_svx(m, p, b);
  if (flagm == 1) gsl_matrix_free(m);
  if (flagp == 1) gsl_permutation_free(p);
  return argv[itmp];
}

/* singleton */
static VALUE rb_gsl_linalg_LU_refine(VALUE obj, VALUE vm,
                                     VALUE lu, VALUE pp, VALUE bb,
                                     VALUE xx)
{
  gsl_matrix *m = NULL, *mlu = NULL;
  gsl_permutation *p = NULL;
  gsl_vector *b = NULL, *x = NULL, *r = NULL;
  int flagb = 0;
  VALUE vr;
  CHECK_MATRIX(vm);  CHECK_MATRIX(lu);
  CHECK_PERMUTATION(pp);  CHECK_VECTOR(xx);
  Data_Get_Struct(vm, gsl_matrix, m);
  Data_Get_Struct(lu, gsl_matrix, mlu);
  Data_Get_Struct(pp, gsl_permutation, p);
  if (TYPE(bb) == T_ARRAY) {
    b = make_cvector_from_rarray(bb);
    flagb = 1;
  } else {
    CHECK_VECTOR(bb);
    Data_Get_Struct(bb, gsl_vector, b);
  }
  Data_Get_Struct(xx, gsl_vector, x);
  r = gsl_vector_alloc(m->size1);
  gsl_linalg_LU_refine(m, mlu, p, b, x, r);
  vr = Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, r);
  if (flagb == 1) gsl_vector_free(b);
  return rb_ary_new3(2, xx, vr);
}

#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_linalg_LU_invert_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *na;
  VALUE inv;
  gsl_permutation *p;
  gsl_matrix_view mv1, mv2;
  if (argc != 2) {
    rb_raise(rb_eArgError, "Usage: LU.invert(lu, perm)");
  }
  CHECK_PERMUTATION(argv[1]);
  GetNArray(argv[0], na);
  inv = na_make_object(NA_DFLOAT, 2, na->shape, CLASS_OF(argv[0]));
  mv1 = gsl_matrix_view_array((double*)na->ptr, na->shape[1], na->shape[0]);
  mv2 = gsl_matrix_view_array(NA_PTR_TYPE(inv, double*), na->shape[1], na->shape[0]);
  CHECK_PERMUTATION(argv[1]);
  Data_Get_Struct(argv[1], gsl_permutation, p);
  gsl_linalg_LU_invert(&mv1.matrix, p, &mv2.matrix);
  return inv;
}
static VALUE rb_gsl_linalg_LU_det_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *na;
  gsl_matrix_view mv;
  int signum = 1;
  switch (argc) {
  case 2:
    signum = FIX2INT(argv[1]);
  /* no break */
  case 1:
    GetNArray(argv[0], na);
    mv = gsl_matrix_view_array((double*)na->ptr, na->shape[1], na->shape[0]);
    break;
  default:
    rb_raise(rb_eArgError, "Usage: LU.det(lu, perm)");
    break;
  }
  return rb_float_new(gsl_linalg_LU_det(&mv.matrix, signum));
}
static VALUE rb_gsl_linalg_LU_lndet_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *na;
  gsl_matrix_view mv;
  switch (argc) {
  case 1:
    GetNArray(argv[0], na);
    mv = gsl_matrix_view_array((double*)na->ptr, na->shape[1], na->shape[0]);
    break;
  default:
    rb_raise(rb_eArgError, "Usage: LU.lndet(lu)");
    break;
  }
  return rb_float_new(gsl_linalg_LU_lndet(&mv.matrix));
}

#endif

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_LU_invert_nmatrix(int argc, VALUE *argv, VALUE obj)
{
  NM_DENSE_STORAGE* lu_nmatrix;
  VALUE inv;
  gsl_permutation *p;
  gsl_matrix_view mv1, mv2;

  if (argc != 2) {
    rb_raise(rb_eArgError, "Usage: LU.invert(lu, perm)");
  }

  CHECK_PERMUTATION(argv[1]);
  lu_nmatrix = NM_STORAGE_DENSE(argv[0]);
  inv = rb_nmatrix_dense_create(FLOAT64, lu_nmatrix->shape, 2, 
    lu_nmatrix->elements, NM_DENSE_COUNT(argv[0]));
  mv1 = gsl_matrix_view_array((double*)lu_nmatrix->elements, 
    lu_nmatrix->shape[0], lu_nmatrix->shape[1]);
  mv2 = gsl_matrix_view_array((double*)NM_DENSE_ELEMENTS(inv), 
    lu_nmatrix->shape[0], lu_nmatrix->shape[1]);

  Data_Get_Struct(argv[1], gsl_permutation, p);
  gsl_linalg_LU_invert(&mv1.matrix, p, &mv2.matrix);
  return inv;
}
#endif

static VALUE rb_gsl_linalg_LU_invert(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *m = NULL, *inverse = NULL;
  gsl_permutation *p = NULL;
  int signum, flagm = 0, flagp = 0, itmp;
  size_t size;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0]))
      return rb_gsl_linalg_LU_invert_narray(argc, argv, obj);
#endif

#ifdef HAVE_NMATRIX_H
    if (NM_IsNMatrix(argv[0])) {
      return rb_gsl_linalg_LU_invert_nmatrix(argc, argv, obj);
    }
#endif
    m = get_matrix(argv[0], cgsl_matrix_LU, &flagm);
    itmp = 1;
    break;
  default:
    m = get_matrix(obj, cgsl_matrix_LU, &flagm);
    itmp = 0;
  }
  size = m->size1;

  if (argc == itmp) {
    p = gsl_permutation_alloc(size);
    flagp = 1;
  } else {
    CHECK_PERMUTATION(argv[itmp]);
    p = get_permutation(argv[itmp], size, &flagp);
  }
  if (flagp == 1 && flagm == 0) rb_raise(rb_eArgError, "permutation must be given");
  if (flagp == 0) itmp++;
  if (flagm == 1 || flagp == 1) {
    gsl_linalg_LU_decomp(m, p, &signum);
  }
  if (argc-1 == itmp) {
    CHECK_MATRIX(argv[itmp]);
    Data_Get_Struct(argv[itmp], gsl_matrix, inverse);
  } else {
    inverse = gsl_matrix_alloc(size, size);
  }
  gsl_linalg_LU_invert(m, p, inverse);
  if (flagm == 1) gsl_matrix_free(m);
  if (flagp == 1) gsl_permutation_free(p);
  if (argc-1 == itmp) return argv[itmp];
  else return Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, inverse);
}

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_LU_det_nmatrix(int argc, VALUE *argv, VALUE obj)
{
  NM_DENSE_STORAGE *input_nmatrix;
  gsl_matrix_view mv;
  int signum = 1;

  switch (argc) {
  case 2:
    signum = FIX2INT(argv[1]);
  /* no break */
  case 1:
    input_nmatrix = NM_STORAGE_DENSE(argv[0]);
    mv = gsl_matrix_view_array((double*)input_nmatrix->elements, 
      input_nmatrix->shape[0], input_nmatrix->shape[1]);
    break;
  default:
    rb_raise(rb_eArgError, "Usage: LU.det(lu, perm)");
    break;
  }
  return rb_float_new(gsl_linalg_LU_det(&mv.matrix, signum));
}
#endif

static VALUE rb_gsl_linalg_LU_det(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *m = NULL;
  gsl_permutation *p = NULL;
  int flagm = 0, flagp = 0, itmp, sign;
  size_t size;
  double det;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc < 1) rb_raise(rb_eArgError, "wrong number of arguments (%d for 1)",
                           argc);
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0]))
      return rb_gsl_linalg_LU_det_narray(argc, argv, obj);
#endif

#ifdef HAVE_NMATRIX_H
    if (NM_IsNMatrix(argv[0])) {
      return rb_gsl_linalg_LU_det_nmatrix(argc, argv, obj);
    }
#endif
    m = get_matrix(argv[0], cgsl_matrix_LU, &flagm);
    itmp = 1;
    break;
  default:
    m = get_matrix(obj, cgsl_matrix_LU, &flagm);
    itmp = 0;
    break;
  }
  size = m->size1;
  if (flagm == 0) {
    if (argc-itmp == 1)  sign = FIX2INT(argv[itmp]);
    else sign = 1;
  } else {
    if (argc-itmp >= 2) {
      get_permutation(argv[itmp], size, &flagp);
    } else {
      p = gsl_permutation_alloc(size);
      flagp = 1;
    }
  }
  if (flagm == 1) gsl_linalg_LU_decomp(m, p, &sign);
  det = gsl_linalg_LU_det(m, sign);
  if (flagm == 1) gsl_matrix_free(m);
  if (flagp == 1) gsl_permutation_free(p);
  return rb_float_new(det);
}

static VALUE rb_gsl_linalg_LU_lndet(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *m = NULL;
  gsl_permutation *p = NULL;
  int flagm = 0, sign;
  double lndet;

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc < 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                           argc);
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0]))
      return rb_gsl_linalg_LU_lndet_narray(argc, argv, obj);
#endif

    m = get_matrix(argv[0], cgsl_matrix_LU, &flagm);
    break;
  default:
    m = get_matrix(obj, cgsl_matrix_LU, &flagm);
    break;
  }
  if (flagm == 1) {
    p = gsl_permutation_alloc(m->size1);
    gsl_linalg_LU_decomp(m, p, &sign);
  }
  lndet = gsl_linalg_LU_lndet(m);
  if (flagm == 1) {
    gsl_matrix_free(m);
    gsl_permutation_free(p);
  }
  return rb_float_new(lndet);
}

static VALUE rb_gsl_linalg_LU_sgndet(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *m = NULL;
  gsl_permutation *p = NULL;
  int flagm = 0, sign, signdet, itmp;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    m = get_matrix(argv[0], cgsl_matrix_LU, &flagm);
    itmp = 1;
    break;
  default:
    m = get_matrix(obj, cgsl_matrix_LU, &flagm);
    itmp = 0;
    break;
  }
  if (flagm == 1) {
    p = gsl_permutation_alloc(m->size1);
    gsl_linalg_LU_decomp(m, p, &sign);
  } else {
    if (argc-itmp != 1) rb_raise(rb_eArgError, "sign must be given");
    sign = FIX2INT(argv[itmp]);
  }
  signdet = gsl_linalg_LU_sgndet(m, sign);
  if (flagm == 1) {
    gsl_matrix_free(m);
    gsl_permutation_free(p);
  }
  return INT2FIX(signdet);
}

int gsl_linalg_LQ_solve_T(const gsl_matrix*, const gsl_vector*, const gsl_vector*, gsl_vector*);
int gsl_linalg_LQ_svx_T (const gsl_matrix*, const gsl_vector*, gsl_vector*);
int gsl_linalg_LQ_lssolve_T(const gsl_matrix * LQ, const gsl_vector * tau,
                            const gsl_vector * b, gsl_vector * x,
                            gsl_vector * residual);
int
gsl_linalg_LQ_Lsolve_T (const gsl_matrix * LQ, const gsl_vector * b, gsl_vector* x);
int
gsl_linalg_LQ_Lsvx_T (const gsl_matrix * LQ, gsl_vector * x);
int
gsl_linalg_L_solve_T (const gsl_matrix * L, const gsl_vector * b, gsl_vector * x);

enum {
  LINALG_QR_DECOMP,
  LINALG_QR_DECOMP_BANG,
  LINALG_LQ_DECOMP,
  LINALG_LQ_DECOMP_BANG,
  LINALG_QR_SOLVE,
  LINALG_LQ_SOLVE,
  LINALG_QR_QTvec,
  LINALG_QR_Qvec,
  LINALG_LQ_vecQ,
  LINALG_LQ_vecQT,
  LINALG_QR_RSOLVE,
  LINALG_LQ_LSOLVE,
  LINALG_QR_RSVX,
  LINALG_LQ_LSVX,
  LINALG_R_SOLVE,
  LINALG_R_SVX,
  LINALG_L_SOLVE,
  LINALG_L_SVX,
  LINALG_QR_UNPACK,
  LINALG_LQ_UNPACK,
};

static VALUE rb_gsl_linalg_QR_LQ_decomposition(int argc, VALUE *argv, VALUE obj,
                                               int flag)
{
  gsl_matrix *m = NULL, *mtmp = NULL;
  gsl_vector *tau = NULL;
  int (*fdecomp)(gsl_matrix *, gsl_vector *);
  int itmp, status;
  VALUE vtau, mdecomp, omatrix;

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc < 1) rb_raise(rb_eArgError, "too few arguments.");
    omatrix = argv[0];
    itmp = 1;
    break;
  default:
    omatrix = obj;
    itmp = 0;
    break;
  }
  CHECK_MATRIX(omatrix);
  Data_Get_Struct(omatrix, gsl_matrix, mtmp);

  switch (flag) {
  case LINALG_QR_DECOMP:
    fdecomp = &gsl_linalg_QR_decomp;
    m = make_matrix_clone(mtmp);
    mdecomp = Data_Wrap_Struct(cgsl_matrix_QR, 0, gsl_matrix_free, m);
    break;
  case LINALG_QR_DECOMP_BANG:
    fdecomp = &gsl_linalg_QR_decomp;
    m = mtmp;
    mdecomp = omatrix;
    RBGSL_SET_CLASS(mdecomp, cgsl_matrix_QR);
    break;
  case LINALG_LQ_DECOMP:
    fdecomp = &gsl_linalg_LQ_decomp;
    m = make_matrix_clone(mtmp);
    mdecomp = Data_Wrap_Struct(cgsl_matrix_LQ, 0, gsl_matrix_free, m);
    break;
  case LINALG_LQ_DECOMP_BANG:
    fdecomp = &gsl_linalg_LQ_decomp;
    m = mtmp;
    mdecomp = omatrix;
    RBGSL_SET_CLASS(mdecomp, cgsl_matrix_LQ);
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  switch (argc - itmp) {
  case 0:
    tau = gsl_vector_alloc(GSL_MIN(mtmp->size1, mtmp->size2));
    break;
  case 1:
    CHECK_VECTOR(argv[itmp]);
    Data_Get_Struct(argv[itmp], gsl_vector, tau);
    break;
  default:
    rb_raise(rb_eArgError, "wrong number of arguments");
    break;
  }
  status = (*fdecomp)(m, tau);
  switch (flag) {
  case LINALG_QR_DECOMP:
  case LINALG_LQ_DECOMP:
    if (argc == itmp) {
      vtau = Data_Wrap_Struct(cgsl_vector_tau, 0, gsl_vector_free, tau);
      return rb_ary_new3(2, mdecomp, vtau);
    } else {
      RBGSL_SET_CLASS(argv[itmp], cgsl_vector_tau);
      return mdecomp;
    }
    break;
  case LINALG_QR_DECOMP_BANG:
  case LINALG_LQ_DECOMP_BANG:
    if (argc == itmp) {
      return Data_Wrap_Struct(cgsl_vector_tau, 0, gsl_vector_free, tau);
    } else {
      RBGSL_SET_CLASS(argv[itmp], cgsl_vector_tau);
      return INT2FIX(status);
    }
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  return Qnil;
}

#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_linalg_QR_unpack_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *m, *tau;
  gsl_matrix_view mv, mq, mr;
  gsl_vector_view vv;
  int shape[2];
  VALUE q, r;
  if (argc != 2) rb_raise(rb_eArgError, "wrong number of arguments (%d for 2)",
                          argc);
  GetNArray(argv[0], m);
  GetNArray(argv[1], tau);
  mv = gsl_matrix_view_array((double*)m->ptr, m->shape[1], m->shape[0]);
  vv = gsl_vector_view_array((double*)tau->ptr, tau->shape[0]);
  shape[0] = m->shape[1];
  shape[1] = m->shape[1];
  q = na_make_object(NA_DFLOAT, 2, shape, CLASS_OF(argv[0]));
  shape[0] = m->shape[1];
  shape[1] = m->shape[0];
  r = na_make_object(NA_DFLOAT, 2, shape, CLASS_OF(argv[0]));
  mq = gsl_matrix_view_array(NA_PTR_TYPE(q,double*), m->shape[1], m->shape[1]);
  mr = gsl_matrix_view_array(NA_PTR_TYPE(r,double*), m->shape[1], m->shape[0]);
  //  printf("OK 4 %d %d\n", mq.matrix.size1, mr.matrix.size2);
  gsl_linalg_QR_unpack(&mv.matrix, &vv.vector, &mq.matrix, &mr.matrix);
  //  printf("OK 5\n");
  return rb_ary_new3(2, q, r);
}

static VALUE rb_gsl_linalg_QR_solve_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *qr, *tau, *b;
  VALUE x;
  gsl_matrix_view mv;
  gsl_vector_view tv, bv, xv;
  if (argc != 3) rb_raise(rb_eArgError, "Usage: QR.solve(qr, tau, b)");
  GetNArray(argv[0], qr);
  GetNArray(argv[1], tau);
  GetNArray(argv[2], b);
  x = na_make_object(NA_DFLOAT, 1, b->shape, CLASS_OF(argv[2]));
  mv = gsl_matrix_view_array((double*)qr->ptr, qr->shape[1], qr->shape[0]);
  tv = gsl_vector_view_array((double*)tau->ptr, tau->shape[0]);
  bv = gsl_vector_view_array((double*)b->ptr, b->shape[0]);
  xv = gsl_vector_view_array(NA_PTR_TYPE(x,double*), b->shape[0]);
  gsl_linalg_QR_solve(&mv.matrix, &tv.vector, &bv.vector, &xv.vector);
  return x;
}
static VALUE rb_gsl_linalg_QR_svx_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *qr, *tau, *b;
  gsl_matrix_view mv;
  gsl_vector_view tv, bv;
  if (argc != 3) rb_raise(rb_eArgError, "Usage: QR.solve(qr, tau, b)");
  GetNArray(argv[0], qr);
  GetNArray(argv[1], tau);
  GetNArray(argv[2], b);
  mv = gsl_matrix_view_array((double*)qr->ptr, qr->shape[1], qr->shape[0]);
  tv = gsl_vector_view_array((double*)tau->ptr, tau->shape[0]);
  bv = gsl_vector_view_array((double*)b->ptr, b->shape[0]);
  gsl_linalg_QR_svx(&mv.matrix, &tv.vector, &bv.vector);
  return argv[2];
}
static VALUE rb_gsl_linalg_QR_decomp_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *na;
  gsl_matrix_view mv;
  gsl_vector_view vv;
  int shapem[2], shapev[1];
  VALUE qr, tau;
  if (argc < 1) rb_raise(rb_eArgError, "too few arguments.");
  GetNArray(argv[0], na);
  shapem[0] = na->shape[1];
  shapem[1] = na->shape[1];
  shapev[0] = shapem[0];
  qr = na_make_object(NA_DFLOAT, 2, shapem, CLASS_OF(argv[0]));
  tau = na_make_object(NA_DFLOAT, 1, shapev, cNVector);
  memcpy(NA_PTR_TYPE(qr,double*),na->ptr,sizeof(double)*shapem[0]*shapem[1]);
  mv = gsl_matrix_view_array(NA_PTR_TYPE(qr,double*), shapem[0], shapem[1]);
  vv = gsl_vector_view_array(NA_PTR_TYPE(tau,double*), shapev[0]);
  gsl_linalg_QR_decomp(&mv.matrix, &vv.vector);
  return rb_ary_new3(2, qr, tau);
}
#endif

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_QR_decomp_nmatrix(int argc, VALUE *argv, VALUE obj)
{
  NM_DENSE_STORAGE *nm;
  gsl_matrix_view mv;
  gsl_vector_view vv;
  unsigned shapem[2], shapev[1];
  VALUE qr, tau;

  if (argc < 1) rb_raise(rb_eArgError, "too few arguments.");
  nm = NM_STORAGE_DENSE(argv[0]);
  shapem[0] = nm->shape[1];
  shapem[1] = nm->shape[1];
  shapev[0] = shapem[0];
  qr = rb_nmatrix_dense_create(FLOAT64, shapem, 2, nm->elements, 
    shapem[0]*shapem[1]);
  tau = rb_nmatrix_dense_create(FLOAT64, shapev, 1, nm->elements, shapev[0]);
  mv = gsl_matrix_view_array((double*)NM_DENSE_ELEMENTS(qr), shapem[0], shapem[1]);
  vv = gsl_vector_view_array((double*)NM_DENSE_ELEMENTS(tau), shapev[0]);
  gsl_linalg_QR_decomp(&mv.matrix, &vv.vector);
  return rb_ary_new3(2, qr, tau);
}
#endif

static VALUE rb_gsl_linalg_QR_decomp(int argc, VALUE *argv, VALUE obj)
{
#ifdef HAVE_NARRAY_H
  if (argc >= 1 && NA_IsNArray(argv[0]))
    return rb_gsl_linalg_QR_decomp_narray(argc, argv, obj);
#endif

#ifdef HAVE_NMATRIX_H
  if (argc >= 1 && NM_IsNMatrix(argv[0]))
    return rb_gsl_linalg_QR_decomp_nmatrix(argc, argv, obj);
#endif

  return rb_gsl_linalg_QR_LQ_decomposition(argc, argv, obj, LINALG_QR_DECOMP);
}

static VALUE rb_gsl_linalg_QR_decomp_bang(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QR_LQ_decomposition(argc, argv, obj, LINALG_QR_DECOMP_BANG);
}

static VALUE rb_gsl_linalg_LQ_decomp(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QR_LQ_decomposition(argc, argv, obj, LINALG_LQ_DECOMP);
}

static VALUE rb_gsl_linalg_LQ_decomp_bang(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QR_LQ_decomposition(argc, argv, obj, LINALG_LQ_DECOMP_BANG);
}

static VALUE rb_gsl_linalg_QR_LQ_solve(int argc, VALUE *argv, VALUE obj, int flag)
{
  gsl_matrix *m = NULL;
  gsl_vector *b = NULL, *x = NULL, *tau = NULL;
  VALUE omatrix;
  int flagm = 0, flagt = 0, flagb = 0, flagx = 0, itmp;
  size_t size;
  int (*fdecomp)(gsl_matrix*, gsl_vector*);
  int (*fsolve)(const gsl_matrix*, const gsl_vector*, const gsl_vector*, gsl_vector*);

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc < 1) rb_raise(rb_eArgError, "too few arguments.");
    omatrix = argv[0];
    itmp = 1;
    break;
  default:
    omatrix = obj;
    itmp = 0;
    break;
  }
  if (argc-itmp < 1 || argc-itmp > 3)
    rb_raise(rb_eArgError,  "wrong number of arguments");
  CHECK_MATRIX(omatrix);
  switch (flag) {
  case LINALG_QR_SOLVE:
    m = get_matrix(omatrix, cgsl_matrix_QR, &flagm);
    fdecomp = &gsl_linalg_QR_decomp;
    fsolve = &gsl_linalg_QR_solve;
    break;
  case LINALG_LQ_SOLVE:
    m = get_matrix(omatrix, cgsl_matrix_LQ, &flagm);
    fdecomp = &gsl_linalg_LQ_decomp;
    fsolve = &gsl_linalg_LQ_solve_T;
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operatioin");
    break;
  }
  size = m->size1;
  if (flagm == 0) { /* the matrix given is already decomped */
    if (CLASS_OF(argv[itmp]) != cgsl_vector_tau)
      rb_raise(rb_eArgError, "tau vector must be given");
    Data_Get_Struct(argv[itmp], gsl_vector, tau);
    flagt = 0;
    itmp++;
  } else {
    if (CLASS_OF(argv[itmp]) == cgsl_vector_tau) {
      Data_Get_Struct(argv[itmp], gsl_vector, tau);
      flagt = 0;
      itmp++;
    } else {
      tau = gsl_vector_alloc(size);
      flagt = 1;
    }
  }
  b = get_vector2(argv[itmp], &flagb);
  itmp++;
  if (itmp == argc) {
    x = gsl_vector_alloc(m->size1);
    flagx = 1;
  } else {
    CHECK_VECTOR(argv[itmp]);
    Data_Get_Struct(argv[itmp], gsl_vector, x);
    flagx = 0;
  }
  if (flagm == 1) (*fdecomp)(m, tau);
  (*fsolve)(m, tau, b, x);
  if (flagm == 1) gsl_matrix_free(m);
  if (flagt == 1) gsl_vector_free(tau);
  if (flagb == 1) gsl_vector_free(b);
  if (flagx == 1) return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
  else return argv[itmp];
}

static VALUE rb_gsl_linalg_QR_LQ_svx(int argc, VALUE *argv, VALUE obj, int flag)
{
  gsl_matrix *m = NULL;
  gsl_vector *b = NULL, *tau = NULL;
  VALUE omatrix;
  int flagm = 0, flagt = 0, flagb = 0, itmp;
  size_t size;
  int (*fdecomp)(gsl_matrix*, gsl_vector*);
  int (*fsvx)(const gsl_matrix*, const gsl_vector*, gsl_vector*);

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc < 1) rb_raise(rb_eArgError, "too few arguments.");
    omatrix = argv[0];
    itmp = 1;
    break;
  default:
    omatrix = obj;
    itmp = 0;
    break;
  }
  if (argc-itmp < 1 || argc-itmp > 2)
    rb_raise(rb_eArgError,  "wrong number of arguments");
  CHECK_MATRIX(omatrix);
  switch (flag) {
  case LINALG_QR_SOLVE:
    m = get_matrix(omatrix, cgsl_matrix_QR, &flagm);
    fdecomp = &gsl_linalg_QR_decomp;
    fsvx = &gsl_linalg_QR_svx;
    break;
  case LINALG_LQ_SOLVE:
    m = get_matrix(omatrix, cgsl_matrix_LQ, &flagm);
    fdecomp = &gsl_linalg_LQ_decomp;
    fsvx = &gsl_linalg_LQ_svx_T;
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operatioin");
    break;
  }
  size = m->size1;
  if (flagm == 0) { /* the matrix given is already decomped */
    if (CLASS_OF(argv[itmp]) != cgsl_vector_tau)
      rb_raise(rb_eArgError, "tau vector must be given");
    Data_Get_Struct(argv[itmp], gsl_vector, tau);
    flagt = 0;
    itmp++;
  } else {
    if (CLASS_OF(argv[itmp]) == cgsl_vector_tau) {
      Data_Get_Struct(argv[itmp], gsl_vector, tau);
      flagt = 0;
      itmp++;
    } else {
      tau = gsl_vector_alloc(size);
      flagt = 1;
    }
  }
  b = get_vector2(argv[itmp], &flagb);
  if (flagm == 1 && flagt == 1) (*fdecomp)(m, tau);
  (*fsvx)(m, tau, b);
  if (flagm == 1) gsl_matrix_free(m);
  if (flagt == 1) gsl_vector_free(tau);
  return argv[itmp];
}

static VALUE rb_gsl_linalg_QR_LQ_lssolve(int argc, VALUE *argv, VALUE obj, int flag)
{
  gsl_matrix *m = NULL;
  gsl_vector *b = NULL, *x = NULL, *tau = NULL, *r = NULL;
  VALUE omatrix;
  int flagm = 0, flagt = 0, flagb = 0, itmp, status;
  size_t size;
  int (*fdecomp)(gsl_matrix*, gsl_vector*);
  int (*flssolve)(const gsl_matrix*, const gsl_vector*, const gsl_vector*, gsl_vector*,
                  gsl_vector*);

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc < 1) rb_raise(rb_eArgError, "too few arguments.");
    omatrix = argv[0];
    itmp = 1;
    break;
  default:
    omatrix = obj;
    itmp = 0;
    break;
  }
  if (argc-itmp < 1 || argc-itmp > 4)
    rb_raise(rb_eArgError,  "wrong number of arguments");
  CHECK_MATRIX(omatrix);
  switch (flag) {
  case LINALG_QR_SOLVE:
    m = get_matrix(omatrix, cgsl_matrix_QR, &flagm);
    fdecomp = &gsl_linalg_QR_decomp;
    flssolve = &gsl_linalg_QR_lssolve;
    break;
  case LINALG_LQ_SOLVE:
    m = get_matrix(omatrix, cgsl_matrix_LQ, &flagm);
    fdecomp = &gsl_linalg_LQ_decomp;
    flssolve = &gsl_linalg_LQ_lssolve_T;
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operatioin");
    break;
  }
  size = m->size1;
  if (flagm == 0) { /* the matrix given is already decomped */
    if (CLASS_OF(argv[itmp]) != cgsl_vector_tau)
      rb_raise(rb_eArgError, "tau vector must be given");
    Data_Get_Struct(argv[itmp], gsl_vector, tau);
    flagt = 0;
    itmp++;
  } else {
    if (CLASS_OF(argv[itmp]) == cgsl_vector_tau) {
      Data_Get_Struct(argv[itmp], gsl_vector, tau);
      flagt = 0;
      itmp++;
    } else {
      tau = gsl_vector_alloc(size);
      flagt = 1;
    }
  }
  b = get_vector2(argv[itmp], &flagb);
  itmp++;
  switch (argc - itmp) {
  case 2:
    CHECK_VECTOR(argv[argc-2]);
    Data_Get_Struct(argv[argc-2], gsl_vector, x);
    CHECK_VECTOR(argv[argc-1]);
    Data_Get_Struct(argv[argc-1], gsl_vector, r);
    break;
  case 1:
    CHECK_VECTOR(argv[argc-1]);
    Data_Get_Struct(argv[argc-1], gsl_vector, x);
    r = gsl_vector_alloc(x->size);
    break;
  case 0:
    x = gsl_vector_alloc(m->size1);
    r = gsl_vector_alloc(m->size1);
    break;
  default:
    rb_raise(rb_eArgError, "wrong number of arguments");
    break;
  }
  if (flagm == 1) (*fdecomp)(m, tau);
  status = (*flssolve)(m, tau, b, x, r);
  if (flagm == 1) gsl_matrix_free(m);
  if (flagt == 1) gsl_vector_free(tau);
  if (flagb == 1) gsl_vector_free(b);
  switch (argc - itmp) {
  case 2:
    return INT2FIX(status);
    break;
  case 1:
    return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, r);
    break;
  default:
    return rb_ary_new3(2, Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x),
                       Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, r));
  }
  return Qnil;
}

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_QR_solve_nmatrix(int argc, VALUE *argv, VALUE obj)
{
  NM_DENSE_STORAGE *qr, *tau, *b;
  VALUE x;
  gsl_matrix_view mv;
  gsl_vector_view tv, bv, xv;

  if (argc != 3) rb_raise(rb_eArgError, "Usage: QR.solve(qr, tau, b)");
  qr = NM_STORAGE_DENSE(argv[0]);
  tau = NM_STORAGE_DENSE(argv[1]);
  b = NM_STORAGE_DENSE(argv[2]);

  x = rb_nmatrix_dense_create(FLOAT64, b->shape, 1, b->elements, b->shape[0]);
  mv = gsl_matrix_view_array((double*)qr->elements, qr->shape[0], qr->shape[1]);
  tv = gsl_vector_view_array((double*)tau->elements, tau->shape[0]);
  bv = gsl_vector_view_array((double*)b->elements, b->shape[0]);
  xv = gsl_vector_view_array((double*)NM_DENSE_ELEMENTS(x), b->shape[0]);
  gsl_linalg_QR_solve(&mv.matrix, &tv.vector, &bv.vector, &xv.vector);
  return x;
}
#endif

static VALUE rb_gsl_linalg_QR_solve(int argc, VALUE *argv, VALUE obj)
{
#ifdef HAVE_NARRAY_H
  if (argc == 3 && NA_IsNArray(argv[0]))
    return rb_gsl_linalg_QR_solve_narray(argc, argv, obj);
#endif

#ifdef HAVE_NMATRIX_H
  if (argc == 3 && NM_IsNMatrix(argv[0])) {
    return rb_gsl_linalg_QR_solve_nmatrix(argc, argv, obj);
  }
#endif
  return rb_gsl_linalg_QR_LQ_solve(argc, argv, obj, LINALG_QR_SOLVE);
}

static VALUE rb_gsl_linalg_QR_svx(int argc, VALUE *argv, VALUE obj)
{
#ifdef HAVE_NARRAY_H
  if (argc == 2 && NA_IsNArray(argv[0]))
    return rb_gsl_linalg_QR_svx_narray(argc, argv, obj);
#endif
  return rb_gsl_linalg_QR_LQ_svx(argc, argv, obj, LINALG_QR_SOLVE);
}

static VALUE rb_gsl_linalg_QR_lssolve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QR_LQ_lssolve(argc, argv, obj, LINALG_QR_SOLVE);
}

static VALUE rb_gsl_linalg_LQ_solve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QR_LQ_solve(argc, argv, obj, LINALG_LQ_SOLVE);
}

static VALUE rb_gsl_linalg_LQ_svx(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QR_LQ_svx(argc, argv, obj, LINALG_LQ_SOLVE);
}

static VALUE rb_gsl_linalg_LQ_lssolve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QR_LQ_lssolve(argc, argv, obj, LINALG_LQ_SOLVE);
}

static VALUE rb_gsl_linalg_QRLQ_QTvec(int argc, VALUE *argv, VALUE obj,
                                      int flag)
{
  gsl_matrix *QR = NULL;
  gsl_vector *tau = NULL, *v = NULL;
  VALUE ret;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 3) rb_raise(rb_eArgError,
                            "wrong number of arguments (%d for 3)", argc);
    CHECK_MATRIX(argv[0]); CHECK_VECTOR(argv[1]); CHECK_VECTOR(argv[2]);
    Data_Get_Struct(argv[0], gsl_matrix, QR);
    Data_Get_Struct(argv[1], gsl_vector, tau);
    Data_Get_Struct(argv[2], gsl_vector, v);
    ret = argv[2];
    break;
  default:
    if (argc != 2) rb_raise(rb_eArgError,
                            "wrong number of arguments (%d for 2)", argc);
    CHECK_VECTOR(argv[2]); CHECK_VECTOR(argv[1]);
    Data_Get_Struct(obj, gsl_matrix, QR);
    Data_Get_Struct(argv[0], gsl_vector, tau);
    Data_Get_Struct(argv[1], gsl_vector, v);
    ret = argv[1];
    break;
  }
  switch (flag) {
  case LINALG_QR_QTvec:
    gsl_linalg_QR_QTvec(QR, tau, v);
    break;
  case LINALG_QR_Qvec:
    gsl_linalg_QR_Qvec(QR, tau, v);
    break;
  case LINALG_LQ_vecQ:
    gsl_linalg_LQ_vecQ(QR, tau, v);
    break;
  case LINALG_LQ_vecQT:
    gsl_linalg_LQ_vecQT(QR, tau, v);
    break;
  default:
    break;
  }
  return ret;
}

static VALUE rb_gsl_linalg_QR_QTvec(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_QTvec(argc, argv, obj, LINALG_QR_QTvec);
}

static VALUE rb_gsl_linalg_QR_Qvec(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_QTvec(argc, argv, obj, LINALG_QR_Qvec);
}

static VALUE rb_gsl_linalg_LQ_vecQT(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_QTvec(argc, argv, obj, LINALG_LQ_vecQT);
}

static VALUE rb_gsl_linalg_LQ_vecQ(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_QTvec(argc, argv, obj, LINALG_LQ_vecQ);
}

static VALUE rb_gsl_linalg_QRLQ_unpack(int argc, VALUE *argv, VALUE obj,
                                       int flag)
{
  gsl_matrix *QR = NULL, *Q = NULL, *R = NULL;
  gsl_vector *tau = NULL;
  int itmp;
  VALUE vtmp, vQ, vR, klass;
  switch (flag) {
  case LINALG_QR_UNPACK:
    klass = cgsl_matrix_QR;
    break;
  case LINALG_LQ_UNPACK:
    klass = cgsl_matrix_LQ;
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 2) rb_raise(rb_eArgError,
                            "wrong number of arguments (%d for 2)", argc);
    vtmp = argv[0];
    itmp = 1;
    break;
  default:
    if (argc != 1) rb_raise(rb_eArgError,
                            "wrong number of arguments (%d for 1)", argc);
    vtmp = obj;
    itmp = 0;
    break;
  }
  CHECK_MATRIX(vtmp);
  if (CLASS_OF(vtmp) != klass) {
    rb_raise(rb_eTypeError, "not a QR matrix");
  }
  Data_Get_Struct(vtmp, gsl_matrix, QR);
  if (CLASS_OF(argv[itmp]) != cgsl_vector_tau)
    rb_raise(rb_eTypeError, "tau vector must be given.");
  Data_Get_Struct(argv[itmp], gsl_vector, tau);
  Q = gsl_matrix_alloc(QR->size1, QR->size1);
  R = gsl_matrix_alloc(QR->size1, QR->size2);
  switch (flag) {
  case LINALG_QR_UNPACK:
    gsl_linalg_QR_unpack(QR, tau, Q, R);
    vQ = Data_Wrap_Struct(cgsl_matrix_Q, 0, gsl_matrix_free, Q);
    vR = Data_Wrap_Struct(cgsl_matrix_R, 0, gsl_matrix_free, R);
    break;
  case LINALG_LQ_UNPACK:
    gsl_linalg_LQ_unpack(QR, tau, Q, R);
    vQ = Data_Wrap_Struct(cgsl_matrix_L, 0, gsl_matrix_free, Q);
    vR = Data_Wrap_Struct(cgsl_matrix_Q, 0, gsl_matrix_free, R);
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  return rb_ary_new3(2, vQ, vR);
}

static VALUE rb_gsl_linalg_QR_unpack(int argc, VALUE *argv, VALUE obj)
{
#ifdef HAVE_NARRAY_H
  if (argc == 2 && NA_IsNArray(argv[0]))
    return rb_gsl_linalg_QR_unpack_narray(argc, argv, obj);
#endif
  return rb_gsl_linalg_QRLQ_unpack(argc, argv, obj, LINALG_QR_UNPACK);
}

static VALUE rb_gsl_linalg_LQ_unpack(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_unpack(argc, argv, obj, LINALG_LQ_UNPACK);
}

/* singleton */
static VALUE rb_gsl_linalg_QRLQ_QRLQsolve(int argc, VALUE *argv, VALUE obj,
                                          int flag)
{
  gsl_matrix *Q = NULL, *R = NULL;
  gsl_vector *b = NULL, *x = NULL;
  int (*fsolve)(gsl_matrix*, gsl_matrix *, const gsl_vector*, gsl_vector *);
  int flagb = 0;
  VALUE retval;
  switch (argc) {
  case 3:
    CHECK_MATRIX(argv[0]); CHECK_MATRIX(argv[1]);
    Data_Get_Struct(argv[0], gsl_matrix, Q);
    Data_Get_Struct(argv[1], gsl_matrix, R);
    x = gsl_vector_alloc(Q->size1);
    retval = Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
    break;
  case 4:
    CHECK_MATRIX(argv[0]); CHECK_MATRIX(argv[1]);
    CHECK_VECTOR(argv[3]);
    Data_Get_Struct(argv[0], gsl_matrix, Q);
    Data_Get_Struct(argv[1], gsl_matrix, R);
    Data_Get_Struct(argv[3], gsl_vector, x);
    retval = argv[3];
    break;
  default:
    rb_raise(rb_eArgError, "wrong number of arguments (%d for 3 or 4)", argc);
    break;
  }
  switch (flag) {
  case LINALG_QR_DECOMP:
    if (CLASS_OF(argv[0]) != cgsl_matrix_Q)
      rb_raise(rb_eTypeError, "not a Q matrix");
    if (CLASS_OF(argv[1]) != cgsl_matrix_R)
      rb_raise(rb_eTypeError, "not a R matrix");
    fsolve = &gsl_linalg_QR_QRsolve;
    break;
  case LINALG_LQ_DECOMP:
    /*    if (CLASS_OF(argv[0]) != cgsl_matrix_L)
      rb_raise(rb_eTypeError, "not a L matrix");
    if (CLASS_OF(argv[1]) != cgsl_matrix_Q)
    rb_raise(rb_eTypeError, "not a Q matrix");*/
    fsolve = &gsl_linalg_LQ_LQsolve;
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  if (TYPE(argv[2]) == T_ARRAY) {
    b = make_cvector_from_rarray(argv[2]);
    flagb = 1;
  } else {
    CHECK_VECTOR(argv[2]);
    Data_Get_Struct(argv[2], gsl_vector, b);
  }
  (*fsolve)(Q, R, b, x);
  if (flagb == 1) gsl_vector_free(b);
  return retval;
}

/*****/
static VALUE rb_gsl_linalg_QRLQ_RLsolve(int argc, VALUE *argv, VALUE obj,
                                        int flag)
{
  gsl_matrix *QR = NULL, *mtmp;
  gsl_vector *b = NULL, *x = NULL, *tau = NULL;
  size_t istart;
  int (*fsolve)(const gsl_matrix*, const gsl_vector*, gsl_vector *);
  int flagb = 0, flagq = 0;
  VALUE omatrix,retval;
  switch (TYPE(obj)) {
  case T_MODULE: case T_CLASS: case T_OBJECT:
    if (argc < 1) rb_raise(rb_eArgError,  "too few arguments");
    omatrix = argv[0];
    istart = 1;
    break;
  default:
    omatrix = obj;
    istart = 0;
    break;
  }
  CHECK_MATRIX(omatrix);
  Data_Get_Struct(omatrix, gsl_matrix, mtmp);
  switch (argc - istart) {
  case 1:
    x = gsl_vector_alloc(mtmp->size1);
    retval = Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
    break;
  case 2:
    Data_Get_Struct(argv[istart+1], gsl_vector, x);
    retval = argv[istart+1];
    break;
  default:
    rb_raise(rb_eArgError, "wrong number of arguments (%d for 3 or 4)", argc);
    break;
  }
  QR = mtmp; flagq = 0;
  switch (flag) {
  case LINALG_QR_RSOLVE:
    if (CLASS_OF(omatrix) != cgsl_matrix_QR) {
      QR = make_matrix_clone(mtmp);
      tau = gsl_vector_alloc(QR->size1);
      gsl_linalg_QR_decomp(QR, tau);
      flagq = 1;
    }
    fsolve = &gsl_linalg_QR_Rsolve;
    break;
  case LINALG_R_SOLVE:
    if (CLASS_OF(omatrix) != cgsl_matrix_QR) {
      QR = make_matrix_clone(mtmp);
      tau = gsl_vector_alloc(QR->size1);
      gsl_linalg_QR_decomp(QR, tau);
      flagq = 1;
    }
    fsolve = &gsl_linalg_R_solve;
    break;
  case LINALG_LQ_LSOLVE:
    if (CLASS_OF(omatrix) != cgsl_matrix_LQ) {
      QR = make_matrix_clone(mtmp);
      tau = gsl_vector_alloc(QR->size1);
      gsl_linalg_LQ_decomp(QR, tau);
      flagq = 1;
    }
    fsolve = &gsl_linalg_LQ_Lsolve_T;
    break;
  case LINALG_L_SOLVE:
    if (CLASS_OF(omatrix) != cgsl_matrix_LQ) {
      QR = make_matrix_clone(mtmp);
      tau = gsl_vector_alloc(QR->size1);
      gsl_linalg_LQ_decomp(QR, tau);
      flagq = 1;
    }
    fsolve = &gsl_linalg_L_solve_T;
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  if (TYPE(argv[istart]) == T_ARRAY) {
    b = make_cvector_from_rarray(argv[istart]);
    flagb = 1;
  } else {
    CHECK_VECTOR(argv[istart]);
    Data_Get_Struct(argv[istart], gsl_vector, b);
  }
  (*fsolve)(QR, b, x);
  if (flagb == 1) gsl_vector_free(b);
  if (flagq == 1) {
    gsl_matrix_free(QR);
    gsl_vector_free(tau);
  }
  return retval;
}

static VALUE rb_gsl_linalg_QRLQ_RLsvx(int argc, VALUE *argv, VALUE obj,
                                      int flag)
{
  gsl_matrix *QR = NULL, *mtmp;
  gsl_vector *x = NULL, *tau = NULL;
  size_t istart;
  int (*fsolve)(const gsl_matrix*, gsl_vector *);
  int flagq = 0;
  VALUE omatrix,retval;
  switch (TYPE(obj)) {
  case T_MODULE: case T_CLASS: case T_OBJECT:
    if (argc < 1) rb_raise(rb_eArgError,  "too few arguments");
    omatrix = argv[0];
    istart = 1;
    break;
  default:
    omatrix = obj;
    istart = 0;
    break;
  }
  CHECK_MATRIX(omatrix);
  Data_Get_Struct(omatrix, gsl_matrix, mtmp);
  switch (argc - istart) {
  case 0:
    x = gsl_vector_alloc(mtmp->size1);
    retval = Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
    break;
  case 1:
    Data_Get_Struct(argv[istart+1], gsl_vector, x);
    retval = argv[istart+1];
    break;
  default:
    rb_raise(rb_eArgError, "wrong number of arguments (%d for 3 or 4)", argc);
    break;
  }
  QR = mtmp; flagq = 0;
  switch (flag) {
  case LINALG_QR_RSVX:
    if (CLASS_OF(omatrix) != cgsl_matrix_QR) {
      QR = make_matrix_clone(mtmp);
      tau = gsl_vector_alloc(QR->size1);
      gsl_linalg_QR_decomp(QR, tau);
      flagq = 1;
    }
    fsolve = &gsl_linalg_QR_Rsvx;
    break;
  /*
case LINALG_R_SVX:
  if (CLASS_OF(omatrix) != cgsl_matrix_QR) {
    QR = make_matrix_clone(mtmp);
    tau = gsl_vector_alloc(QR->size1);
    gsl_linalg_QR_decomp(QR, tau);
    flagq = 1;
  }
  fsolve = &gsl_linalg_R_svx;
  break;
  */
  case LINALG_LQ_LSVX:
    if (CLASS_OF(omatrix) != cgsl_matrix_LQ) {
      QR = make_matrix_clone(mtmp);
      tau = gsl_vector_alloc(QR->size1);
      gsl_linalg_LQ_decomp(QR, tau);
      flagq = 1;
    }
    fsolve = &gsl_linalg_LQ_Lsvx_T;
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  (*fsolve)(QR, x);
  if (flagq == 1) {
    gsl_matrix_free(QR);
    gsl_vector_free(tau);
  }
  return retval;
}

static VALUE rb_gsl_linalg_QR_Rsolve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_RLsolve(argc, argv, obj, LINALG_QR_RSOLVE);
}

static VALUE rb_gsl_linalg_QR_Rsvx(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_RLsvx(argc, argv, obj, LINALG_QR_RSVX);
}

static VALUE rb_gsl_linalg_R_solve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_RLsolve(argc, argv, obj, LINALG_R_SOLVE);
}

/* singleton */
static VALUE rb_gsl_linalg_QR_QRsolve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_QRLQsolve(argc, argv, obj, LINALG_QR_DECOMP);
}

static VALUE rb_gsl_linalg_LQ_Lsolve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_RLsolve(argc, argv, obj, LINALG_LQ_LSOLVE);
}

static VALUE rb_gsl_linalg_LQ_Lsvx(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_RLsvx(argc, argv, obj, LINALG_LQ_LSVX);
}

static VALUE rb_gsl_linalg_L_solve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_RLsolve(argc, argv, obj, LINALG_L_SOLVE);
}

/* singleton */
static VALUE rb_gsl_linalg_LQ_LQsolve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQ_QRLQsolve(argc, argv, obj, LINALG_LQ_DECOMP);
}

static VALUE rb_gsl_linalg_QRLQ_update(VALUE obj, VALUE qq, VALUE rr, VALUE ww,
                                       VALUE vv, int flag)
{
  gsl_matrix *Q = NULL, *R = NULL;
  gsl_vector *w = NULL, *v = NULL;
  int status;
  CHECK_MATRIX(qq); CHECK_MATRIX(rr);
  CHECK_VECTOR(ww); CHECK_VECTOR(vv);
  Data_Get_Struct(qq, gsl_matrix, Q);
  Data_Get_Struct(rr, gsl_matrix, R);
  Data_Get_Struct(ww, gsl_vector, w);
  Data_Get_Struct(vv, gsl_vector, v);
  switch (flag) {
  case LINALG_QR_DECOMP:
    status = gsl_linalg_QR_update(Q, R, w, v);
    break;
  case LINALG_LQ_DECOMP:
    status = gsl_linalg_LQ_update(Q, R, w, v);
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  return INT2FIX(status);
}

/* singleton */
static VALUE rb_gsl_linalg_QR_update(VALUE obj, VALUE qq, VALUE rr, VALUE ww,
                                     VALUE vv)
{
  return rb_gsl_linalg_QRLQ_update(obj, qq, rr, ww, vv, LINALG_QR_DECOMP);
}

static VALUE rb_gsl_linalg_LQ_update(VALUE obj, VALUE qq, VALUE rr, VALUE ww,
                                     VALUE vv)
{
  return rb_gsl_linalg_QRLQ_update(obj, qq, rr, ww, vv, LINALG_LQ_DECOMP);
}

/******/
enum {
  LINALG_QRPT,
  LINALG_PTLQ,
};

static VALUE rb_gsl_linalg_QRLQPT_decomp(int argc, VALUE *argv, VALUE obj, int flag)
{
  gsl_matrix *A = NULL, *QR = NULL;
  gsl_vector *tau = NULL, *norm = NULL;
  gsl_permutation *p = NULL;
  int signum;
  size_t size0;
  VALUE vtau, vp, vA, vQR;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of arguments (%d for 1)",
                            argc);
    vA = argv[0];
    break;
  default:
    vA = obj;
    break;
  }
  CHECK_MATRIX(vA);
  Data_Get_Struct(vA, gsl_matrix, A);
  QR = make_matrix_clone(A);
  size0 = GSL_MIN(A->size1, A->size2);
  tau = gsl_vector_alloc(size0);
  p = gsl_permutation_alloc(size0);
  norm = gsl_vector_alloc(size0);
  switch (flag) {
  case LINALG_QRPT:
    vQR = Data_Wrap_Struct(cgsl_matrix_QRPT, 0, gsl_matrix_free, QR);
    vtau = Data_Wrap_Struct(cgsl_vector_tau, 0, gsl_vector_free, tau);
    vp = Data_Wrap_Struct(cgsl_permutation, 0, gsl_permutation_free, p);
    gsl_linalg_QRPT_decomp(QR, tau, p, &signum, norm);
    break;
  case LINALG_PTLQ:
    vQR = Data_Wrap_Struct(cgsl_matrix_PTLQ, 0, gsl_matrix_free, QR);
    vtau = Data_Wrap_Struct(cgsl_vector_tau, 0, gsl_vector_free, tau);
    vp = Data_Wrap_Struct(cgsl_permutation, 0, gsl_permutation_free, p);
    gsl_linalg_PTLQ_decomp(QR, tau, p, &signum, norm);
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  gsl_vector_free(norm);
  return rb_ary_new3(4, vQR, vtau, vp, INT2FIX(signum));
}

static VALUE rb_gsl_linalg_QRLQPT_decomp_bang(int argc, VALUE *argv, VALUE obj, int flag)
{
  gsl_matrix *A = NULL;
  gsl_vector *tau = NULL, *norm = NULL;
  gsl_permutation *p = NULL;
  int signum;
  size_t size0;
  VALUE vtau, vp, vA;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of arguments (%d for 1)",
                            argc);
    vA = argv[0];
    break;
  default:
    vA = obj;
    break;
  }
  CHECK_MATRIX(vA);
  Data_Get_Struct(vA, gsl_matrix, A);
  size0 = GSL_MIN(A->size1, A->size2);
  tau = gsl_vector_alloc(size0);
  p = gsl_permutation_alloc(size0);
  norm = gsl_vector_alloc(size0);
  switch (flag) {
  case LINALG_QRPT:
    RBGSL_SET_CLASS(vA, cgsl_matrix_QRPT);
    vtau = Data_Wrap_Struct(cgsl_vector_tau, 0, gsl_vector_free, tau);
    vp = Data_Wrap_Struct(cgsl_permutation, 0, gsl_permutation_free, p);
    gsl_linalg_QRPT_decomp(A, tau, p, &signum, norm);
    break;
  case LINALG_PTLQ:
    RBGSL_SET_CLASS(vA, cgsl_matrix_PTLQ);
    vtau = Data_Wrap_Struct(cgsl_vector_tau, 0, gsl_vector_free, tau);
    vp = Data_Wrap_Struct(cgsl_permutation, 0, gsl_permutation_free, p);
    gsl_linalg_PTLQ_decomp(A, tau, p, &signum, norm);
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  gsl_vector_free(norm);
  return rb_ary_new3(3, vtau, vp, INT2FIX(signum));
}

static VALUE rb_gsl_linalg_QRPT_decomp(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_decomp(argc, argv, obj, LINALG_QRPT);
}

static VALUE rb_gsl_linalg_QRPT_decomp_bang(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_decomp_bang(argc, argv, obj, LINALG_QRPT);
}

static VALUE rb_gsl_linalg_PTLQ_decomp(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_decomp(argc, argv, obj, LINALG_PTLQ);
}

static VALUE rb_gsl_linalg_PTLQ_decomp_bang(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_decomp_bang(argc, argv, obj, LINALG_PTLQ);
}

static VALUE rb_gsl_linalg_QRLQPT_decomp2(int argc, VALUE *argv, VALUE obj,int flag)
{
  gsl_matrix *A = NULL, *Q = NULL, *R = NULL;
  gsl_vector *tau = NULL, *norm = NULL;
  gsl_permutation *p = NULL;
  int signum;
  size_t size0;
  VALUE vtau, vp, vA, vQ, vR;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of arguments");
    vA = argv[0];
    break;
  default:
    if (argc != 0) rb_raise(rb_eArgError, "wrong number of arguments");
    vA = obj;
    break;
  }
  CHECK_MATRIX(vA);
  Data_Get_Struct(vA, gsl_matrix, A);
  Q = gsl_matrix_alloc(A->size1, A->size2);
  R = gsl_matrix_alloc(A->size1, A->size2);
  size0 = GSL_MIN(A->size1, A->size2);
  tau = gsl_vector_alloc(size0);
  p = gsl_permutation_alloc(size0);
  norm = gsl_vector_alloc(size0);
  /*  vQ = Data_Wrap_Struct(cgsl_matrix_Q, 0, gsl_matrix_free, Q);
      vR = Data_Wrap_Struct(cgsl_matrix_R, 0, gsl_matrix_free, R);*/
  vtau = Data_Wrap_Struct(cgsl_vector_tau, 0, gsl_vector_free, tau);
  vp = Data_Wrap_Struct(cgsl_permutation, 0, gsl_permutation_free, p);
  switch (flag) {
  case LINALG_QRPT:
    vQ = Data_Wrap_Struct(cgsl_matrix_Q, 0, gsl_matrix_free, Q);
    vR = Data_Wrap_Struct(cgsl_matrix_R, 0, gsl_matrix_free, R);
    gsl_linalg_QRPT_decomp2(A, Q, R, tau, p, &signum, norm);
    break;
  case LINALG_PTLQ:
    vR = Data_Wrap_Struct(cgsl_matrix_L, 0, gsl_matrix_free, R);
    vQ = Data_Wrap_Struct(cgsl_matrix_Q, 0, gsl_matrix_free, Q);
    gsl_linalg_PTLQ_decomp2(A, Q, R, tau, p, &signum, norm);
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
  }
  gsl_vector_free(norm);
  return rb_ary_new3(5, vQ, vR, vtau, vp, INT2FIX(signum));
}

static VALUE rb_gsl_linalg_QRPT_decomp2(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_decomp2(argc, argv, obj, LINALG_QRPT);
}

static VALUE rb_gsl_linalg_PTLQ_decomp2(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_decomp2(argc, argv, obj, LINALG_PTLQ);
}

int gsl_linalg_PTLQ_solve_T(const gsl_matrix * QR, const gsl_vector * tau,
                            const gsl_permutation * p, const gsl_vector * b,
                            gsl_vector * x);
int gsl_linalg_PTLQ_svx_T(const gsl_matrix * LQ,
                          const gsl_vector * tau,
                          const gsl_permutation * p,
                          gsl_vector * x);
int gsl_linalg_PTLQ_LQsolve_T (const gsl_matrix * Q, const gsl_matrix * L,
                               const gsl_permutation * p,
                               const gsl_vector * b,
                               gsl_vector * x);
int gsl_linalg_PTLQ_Lsolve_T (const gsl_matrix * LQ,
                              const gsl_permutation * p,
                              const gsl_vector * b,
                              gsl_vector * x);
int gsl_linalg_PTLQ_Lsvx_T (const gsl_matrix * LQ,
                            const gsl_permutation * p,
                            gsl_vector * x);

static VALUE rb_gsl_linalg_QRLQPT_solve(int argc, VALUE *argv, VALUE obj, int flag)
{
  gsl_matrix *QR = NULL, *A = NULL;
  gsl_vector *tau = NULL, *b = NULL, *x = NULL, *norm = NULL;
  gsl_permutation *p = NULL;
  int signum, itmp, flagb = 0, flagq = 0;
  VALUE vtmp, klass;
  size_t size0;
  int (*fdecomp)(gsl_matrix*, gsl_vector*, gsl_permutation*, int *, gsl_vector*);
  int (*fsolve)(const gsl_matrix*, const gsl_vector*, const gsl_permutation*,
                const gsl_vector*, gsl_vector *);
  switch (flag) {
  case LINALG_QRPT:
    klass = cgsl_matrix_QRPT;
    fdecomp = &gsl_linalg_QRPT_decomp;
    fsolve = &gsl_linalg_QRPT_solve;
    break;
  case LINALG_PTLQ:
    klass = cgsl_matrix_PTLQ;
    fdecomp = &gsl_linalg_PTLQ_decomp;
    fsolve = &gsl_linalg_PTLQ_solve_T;
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc < 1) rb_raise(rb_eArgError, "too few arguments.");
    vtmp = argv[0];
    itmp = 1;
    break;
  default:
    vtmp = obj;
    itmp = 0;
    break;
  }
  CHECK_MATRIX(vtmp);
  if (CLASS_OF(vtmp) == klass) {
    if (argc-itmp != 3) rb_raise(rb_eArgError,
                                 "wrong number of arguments (%d for %d)",
                                 argc, 4-itmp);
    CHECK_VECTOR(argv[itmp]);
    if (CLASS_OF(argv[itmp]) != cgsl_vector_tau)
      rb_raise(rb_eTypeError, "not a tau vector");
    CHECK_PERMUTATION(argv[itmp+1]);
    Data_Get_Struct(argv[itmp], gsl_vector, tau);
    Data_Get_Struct(argv[itmp+1], gsl_permutation, p);
    Data_Get_Struct(vtmp, gsl_matrix, QR);
    size0 = GSL_MIN(QR->size1, QR->size2);
    itmp += 2;
  } else {
    if (argc-itmp != 1) rb_raise(rb_eArgError,
                                 "wrong number of arguments (%d for %d)", argc, 2-itmp);
    Data_Get_Struct(vtmp, gsl_matrix, A);
    QR = make_matrix_clone(A);
    size0 = GSL_MIN(QR->size1, QR->size2);
    flagq = 1;
    p = gsl_permutation_alloc(size0);
    tau = gsl_vector_alloc(size0);
  }
  norm = gsl_vector_alloc(size0);
  if (TYPE(argv[itmp]) == T_ARRAY) {
    b = make_cvector_from_rarray(argv[itmp]);
    flagb = 1;
  } else {
    CHECK_VECTOR(argv[itmp]);
    Data_Get_Struct(argv[itmp], gsl_vector, b);
  }
  x = gsl_vector_alloc(b->size);
  if (flagq == 1) (*fdecomp)(QR, tau, p, &signum, norm);
  (*fsolve)(QR, tau, p, b, x);
  if (flagb == 1) gsl_vector_free(b);
  if (flagq == 1) {
    gsl_matrix_free(QR);
    gsl_permutation_free(p);
    gsl_vector_free(tau);
    gsl_vector_free(norm);
  }
  return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
}

static VALUE rb_gsl_linalg_QRPT_solve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_solve(argc, argv, obj, LINALG_QRPT);
}

static VALUE rb_gsl_linalg_PTLQ_solve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_solve(argc, argv, obj, LINALG_PTLQ);
}

static VALUE rb_gsl_linalg_QRLQPT_svx(int argc, VALUE *argv, VALUE obj, int flag)
{
  gsl_matrix *QR = NULL, *A = NULL;
  gsl_vector *tau = NULL, *b = NULL, *norm = NULL;
  gsl_permutation *p = NULL;
  int signum, itmp, flagq = 0;
  VALUE vtmp, klass;
  size_t size0;
  int (*fdecomp)(gsl_matrix*, gsl_vector*, gsl_permutation*, int *, gsl_vector*);
  int (*fsvx)(const gsl_matrix*, const gsl_vector*, const gsl_permutation*,
              gsl_vector *);
  switch (flag) {
  case LINALG_QRPT:
    klass = cgsl_matrix_QRPT;
    fdecomp = &gsl_linalg_QRPT_decomp;
    fsvx = &gsl_linalg_QRPT_svx;
    break;
  case LINALG_PTLQ:
    klass = cgsl_matrix_PTLQ;
    fdecomp = &gsl_linalg_PTLQ_decomp;
    fsvx = &gsl_linalg_PTLQ_svx_T;
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of arguments (%d for 1)",
                            argc);
    vtmp = argv[0];
    itmp = 1;
    break;
  default:
    vtmp = obj;
    itmp = 0;
    break;
  }
  CHECK_MATRIX(vtmp);
  if (CLASS_OF(vtmp) == klass) {
    if (argc-itmp != 3) rb_raise(rb_eArgError,
                                 "wrong number of arguments (%d for %d)",
                                 argc, 3+itmp);
    CHECK_VECTOR(argv[itmp]);
    if (CLASS_OF(argv[itmp]) != cgsl_vector_tau)
      rb_raise(rb_eTypeError, "not a tau vector");
    CHECK_PERMUTATION(argv[itmp+1]);
    Data_Get_Struct(argv[itmp], gsl_vector, tau);
    Data_Get_Struct(argv[itmp+1], gsl_permutation, p);
    Data_Get_Struct(vtmp, gsl_matrix, QR);
    size0 = GSL_MIN(QR->size1, QR->size2);
    itmp += 2;
  } else {
    if (argc-itmp != 1) rb_raise(rb_eArgError,
                                 "wrong number of arguments (%d for %d)", argc, 2+itmp);
    Data_Get_Struct(vtmp, gsl_matrix, A);
    QR = make_matrix_clone(A);
    size0 = GSL_MIN(QR->size1, QR->size2);
    flagq = 1;
    p = gsl_permutation_alloc(size0);
    tau = gsl_vector_alloc(size0);
  }
  norm = gsl_vector_alloc(size0);
  CHECK_VECTOR(argv[itmp]);
  Data_Get_Struct(argv[itmp], gsl_vector, b);
  if (flagq == 1) (*fdecomp)(QR, tau, p, &signum, norm);
  (*fsvx)(QR, tau, p, b);
  if (flagq == 1) {
    gsl_matrix_free(QR);
    gsl_permutation_free(p);
    gsl_vector_free(tau);
    gsl_vector_free(norm);
  }
  return argv[itmp];
}

static VALUE rb_gsl_linalg_QRPT_svx(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_svx(argc, argv, obj, LINALG_QRPT);
}

static VALUE rb_gsl_linalg_PTLQ_svx(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_svx(argc, argv, obj, LINALG_PTLQ);
}

/* singleton */
static VALUE rb_gsl_linalg_QRLQPT_QRLQsolve(VALUE obj, VALUE qq, VALUE rr,
                                            VALUE pp, VALUE bb, int flag)
{
  gsl_matrix *Q = NULL, *R = NULL;
  gsl_vector *b = NULL, *x = NULL;
  gsl_permutation *p = NULL;
  int flagb = 0;
  int (*fsolve)(const gsl_matrix*, const gsl_matrix*, const gsl_permutation*,
                const gsl_vector*, gsl_vector*);
  switch (flag) {
  case LINALG_QRPT:
    if (CLASS_OF(qq) != cgsl_matrix_Q) rb_raise(rb_eTypeError, "not a Q matrix");
    if (CLASS_OF(rr) != cgsl_matrix_R) rb_raise(rb_eTypeError, "not a R matrix");
    fsolve = &gsl_linalg_QRPT_QRsolve;
    break;
  case LINALG_PTLQ:
    if (CLASS_OF(qq) != cgsl_matrix_Q) rb_raise(rb_eTypeError, "not a Q matrix");
    if (CLASS_OF(rr) != cgsl_matrix_L) rb_raise(rb_eTypeError, "not a L matrix");
    fsolve = &gsl_linalg_PTLQ_LQsolve_T;
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  if (TYPE(bb) == T_ARRAY) {
    b = make_cvector_from_rarray(bb);
    flagb = 1;
  } else {
    CHECK_VECTOR(bb);
    Data_Get_Struct(bb, gsl_vector, b);
  }
  CHECK_PERMUTATION(pp);
  Data_Get_Struct(qq, gsl_matrix, Q);
  Data_Get_Struct(rr, gsl_matrix, R);
  Data_Get_Struct(pp, gsl_permutation, p);
  x = gsl_vector_alloc(b->size);
  (*fsolve)(Q, R, p, b, x);
  if (flagb == 1) gsl_vector_free(b);
  return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
}

static VALUE rb_gsl_linalg_QRPT_QRsolve(VALUE obj, VALUE qq, VALUE rr,
                                        VALUE pp, VALUE bb)
{
  return rb_gsl_linalg_QRLQPT_QRLQsolve(obj, qq, rr, pp, bb, LINALG_QRPT);
}

static VALUE rb_gsl_linalg_PTLQ_LQsolve(VALUE obj, VALUE qq, VALUE rr,
                                        VALUE pp, VALUE bb)
{
  return rb_gsl_linalg_QRLQPT_QRLQsolve(obj, qq, rr, pp, bb, LINALG_PTLQ);
}

/* singleton */
static VALUE rb_gsl_linalg_QRLQPT_update(VALUE obj, VALUE qq, VALUE rr,
                                         VALUE pp, VALUE ww, VALUE vv, int flag)
{
  gsl_matrix *Q = NULL, *R = NULL;
  gsl_vector *w = NULL, *v = NULL;
  gsl_permutation *p = NULL;
  switch (flag) {
  case LINALG_QRPT:
    if (CLASS_OF(qq) != cgsl_matrix_Q) rb_raise(rb_eTypeError, "not a Q matrix");
    if (CLASS_OF(rr) != cgsl_matrix_R) rb_raise(rb_eTypeError, "not a R matrix");
    break;
  case LINALG_PTLQ:
    if (CLASS_OF(qq) != cgsl_matrix_Q) rb_raise(rb_eTypeError, "not a Q matrix");
    if (CLASS_OF(rr) != cgsl_matrix_L) rb_raise(rb_eTypeError, "not a L matrix");
    break;
  }
  CHECK_PERMUTATION(pp);
  Data_Get_Struct(qq, gsl_matrix, Q);
  Data_Get_Struct(rr, gsl_matrix, R);
  Data_Get_Struct(pp, gsl_permutation, p);
  Data_Get_Struct(ww, gsl_vector, w);
  Data_Get_Struct(vv, gsl_vector, v);
  switch (flag) {
  case LINALG_QRPT:
    gsl_linalg_QRPT_update(Q, R, p, w, v);
    break;
  case LINALG_PTLQ:
    gsl_linalg_PTLQ_update(Q, R, p, w, v);
    break;
  }
  return obj;
}

static VALUE rb_gsl_linalg_QRPT_update(VALUE obj, VALUE qq, VALUE rr,
                                       VALUE pp, VALUE ww, VALUE vv)
{
  return rb_gsl_linalg_QRLQPT_update(obj, qq, rr, pp, ww, vv, LINALG_QRPT);
}

static VALUE rb_gsl_linalg_PTLQ_update(VALUE obj, VALUE qq, VALUE rr,
                                       VALUE pp, VALUE ww, VALUE vv)
{
  return rb_gsl_linalg_QRLQPT_update(obj, qq, rr, pp, ww, vv, LINALG_PTLQ);
}

static VALUE rb_gsl_linalg_QRLQPT_RLsolve(int argc, VALUE *argv, VALUE obj, int flag)
{
  gsl_matrix *QR = NULL;
  gsl_vector *b = NULL, *x = NULL;
  gsl_permutation *p = NULL;
  int itmp, flagb = 0;
  VALUE vtmp, klass;
  int (*fsolve)(const gsl_matrix*, const gsl_permutation*, const gsl_vector*,
                gsl_vector*);
  switch (flag) {
  case LINALG_QRPT:
    klass = cgsl_matrix_QRPT;
    fsolve = &gsl_linalg_QRPT_Rsolve;
    break;
  case LINALG_PTLQ:
    klass = cgsl_matrix_PTLQ;
    fsolve = &gsl_linalg_PTLQ_Lsolve_T;
    break;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of arguments (%d for 1)",
                            argc);
    vtmp = argv[0];
    itmp = 1;
    break;
  default:
    vtmp = obj;
    itmp = 0;
    break;
  }
  if (argc-itmp != 2)
    rb_raise(rb_eArgError, "wrong number of argument (%d for %d)", argc, 2+itmp);
  CHECK_MATRIX(vtmp);
  if (CLASS_OF(vtmp) != klass) {
    rb_raise(rb_eArgError, "not a QR matrix");
  }
  CHECK_PERMUTATION(argv[itmp]);
  Data_Get_Struct(argv[itmp], gsl_permutation, p);
  Data_Get_Struct(vtmp, gsl_matrix, QR);
  itmp++;
  if (TYPE(argv[itmp]) == T_ARRAY) {
    b = make_cvector_from_rarray(argv[itmp]);
    flagb = 1;
  } else {
    CHECK_VECTOR(argv[itmp]);
    Data_Get_Struct(argv[itmp], gsl_vector, b);
  }
  x = gsl_vector_alloc(b->size);
  (*fsolve)(QR, p, b, x);
  if (flagb == 1) gsl_vector_free(b);
  return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
}

static VALUE rb_gsl_linalg_QRPT_Rsolve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_RLsolve(argc, argv, obj, LINALG_QRPT);
}

static VALUE rb_gsl_linalg_PTLQ_Lsolve(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_RLsolve(argc, argv, obj, LINALG_PTLQ);
}

static VALUE rb_gsl_linalg_QRLQPT_RLsvx(int argc, VALUE *argv, VALUE obj, int flag)
{
  gsl_matrix *QR = NULL;
  gsl_vector *b = NULL;
  gsl_permutation *p = NULL;
  int itmp;
  VALUE vtmp, klass;
  int (*fsvx)(const gsl_matrix*, const gsl_permutation*, gsl_vector*);
  switch (flag) {
  case LINALG_QRPT:
    klass = cgsl_matrix_QRPT;
    fsvx = &gsl_linalg_QRPT_Rsvx;
    break;
  case LINALG_PTLQ:
    klass = cgsl_matrix_PTLQ;
    fsvx = &gsl_linalg_PTLQ_Lsvx_T;
  default:
    rb_raise(rb_eRuntimeError, "unknown operation");
    break;
  }
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of arguments (%d for 1)",
                            argc);
    vtmp = argv[0];
    itmp = 1;
    break;
  default:
    vtmp = obj;
    itmp = 0;
    break;
  }
  if (argc-itmp != 2)
    rb_raise(rb_eArgError, "wrong number of argument (%d for %d)", argc, 2+itmp);
  CHECK_MATRIX(vtmp);
  if (CLASS_OF(vtmp) != klass) {
    rb_raise(rb_eArgError, "not a QR matrix");
  }
  CHECK_PERMUTATION(argv[itmp]);
  Data_Get_Struct(argv[itmp], gsl_permutation, p);
  Data_Get_Struct(vtmp, gsl_matrix, QR);
  itmp++;
  if (TYPE(argv[itmp]) == T_ARRAY) {
    b = make_cvector_from_rarray(argv[itmp]);
  } else {
    CHECK_VECTOR(argv[itmp]);
    Data_Get_Struct(argv[itmp], gsl_vector, b);
  }
  (*fsvx)(QR, p, b);
  return argv[itmp];
}

static VALUE rb_gsl_linalg_QRPT_Rsvx(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_RLsvx(argc, argv, obj, LINALG_QRPT);
}

static VALUE rb_gsl_linalg_PTLQ_Lsvx(int argc, VALUE *argv, VALUE obj)
{
  return rb_gsl_linalg_QRLQPT_RLsvx(argc, argv, obj, LINALG_PTLQ);
}

/*******/
#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_linalg_SV_decomp_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *A;
  gsl_matrix_view uv, vv;
  gsl_vector_view sv;
  gsl_vector *work;
  VALUE u, s, v;
  int shape[2];
  GetNArray(argv[0], A);
  shape[0] = A->shape[0];
  shape[1] = A->shape[0];
  u = na_make_object(NA_DFLOAT, 2, A->shape, CLASS_OF(argv[0]));
  v = na_make_object(NA_DFLOAT, 2, shape, CLASS_OF(argv[0]));
  s = na_make_object(NA_DFLOAT, 1, &(shape[0]), cNVector);
  uv = gsl_matrix_view_array(NA_PTR_TYPE(u,double*), A->shape[1], A->shape[0]);
  vv = gsl_matrix_view_array(NA_PTR_TYPE(v,double*), shape[1], shape[0]);
  sv = gsl_vector_view_array(NA_PTR_TYPE(s,double*), shape[0]);
  work = gsl_vector_alloc(shape[0]);
  memcpy(NA_PTR_TYPE(u,double*), (double*)A->ptr, sizeof(double)*A->total);
  gsl_linalg_SV_decomp(&uv.matrix, &vv.matrix, &sv.vector, work);
  gsl_vector_free(work);
  return rb_ary_new3(3, u, v, s);
}

static VALUE rb_gsl_linalg_SV_decomp_jacobi_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *A;
  gsl_matrix_view uv, vv;
  gsl_vector_view sv;
  VALUE u, s, v;
  int shape[2];
  GetNArray(argv[0], A);
  shape[0] = A->shape[0];
  shape[1] = A->shape[0];
  u = na_make_object(NA_DFLOAT, 2, A->shape, CLASS_OF(argv[0]));
  v = na_make_object(NA_DFLOAT, 2, shape, CLASS_OF(argv[0]));
  s = na_make_object(NA_DFLOAT, 1, &(shape[0]), cNVector);
  uv = gsl_matrix_view_array(NA_PTR_TYPE(u,double*), A->shape[1], A->shape[0]);
  vv = gsl_matrix_view_array(NA_PTR_TYPE(v,double*), shape[1], shape[0]);
  sv = gsl_vector_view_array(NA_PTR_TYPE(s,double*), shape[0]);
  memcpy(NA_PTR_TYPE(u,double*), (double*)A->ptr, sizeof(double)*A->total);
  gsl_linalg_SV_decomp_jacobi(&uv.matrix, &vv.matrix, &sv.vector);
  return rb_ary_new3(3, u, v, s);
}

static VALUE rb_gsl_linalg_SV_solve_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *A;
  gsl_matrix_view uv, vv;
  gsl_vector_view sv, bv, xv;
  VALUE x;
  if (argc != 4)
    rb_raise(rb_eArgError, "Usage: SV.solve(u, v, s, b)");
  GetNArray(argv[0], A);
  uv = gsl_matrix_view_array(NA_PTR_TYPE(argv[0],double*), A->shape[1], A->shape[0]);
  vv = gsl_matrix_view_array(NA_PTR_TYPE(argv[1],double*), A->shape[0], A->shape[0]);
  sv = gsl_vector_view_array(NA_PTR_TYPE(argv[2],double*), A->shape[0]);
  bv = gsl_vector_view_array(NA_PTR_TYPE(argv[3],double*), A->shape[0]);
  x = na_make_object(NA_DFLOAT, 1, &(A->shape[0]), CLASS_OF(argv[3]));
  xv = gsl_vector_view_array(NA_PTR_TYPE(x,double*), A->shape[0]);
  gsl_linalg_SV_solve(&uv.matrix, &vv.matrix, &sv.vector, &bv.vector, &xv.vector);
  return x;
}

#endif

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_SV_decomp_nmatrix(int argc, VALUE *argv, VALUE obj)
{
  NM_DENSE_STORAGE *A;
  gsl_matrix_view uv, vv;
  gsl_vector_view sv;
  gsl_vector *work;
  VALUE u, s, v;
  unsigned shape[2];

  A = NM_STORAGE_DENSE(argv[0]);
  shape[0] = A->shape[0];
  shape[1] = A->shape[1];
  u = rb_nmatrix_dense_create(FLOAT64, A->shape, 2, A->elements, shape[0]*shape[1]);
  v = rb_nmatrix_dense_create(FLOAT64, shape, 2, A->elements, shape[0]*shape[0]);
  s = rb_nvector_dense_create(FLOAT64, A->elements, shape[0]);
  uv = gsl_matrix_view_array((double*)NM_DENSE_ELEMENTS(u), shape[1], shape[0]);
  vv = gsl_matrix_view_array((double*)NM_DENSE_ELEMENTS(v), shape[0], shape[0]);
  sv = gsl_vector_view_array((double*)NM_DENSE_ELEMENTS(s), shape[0]);
  work = gsl_vector_alloc(shape[0]);
  gsl_linalg_SV_decomp(&uv.matrix, &vv.matrix, &sv.vector, work);
  gsl_vector_free(work);
  return rb_ary_new3(3, u, v, s);
}
#endif

static VALUE rb_gsl_linalg_SV_decomp(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *V = NULL, *U = NULL;
  gsl_vector *w = NULL, *S = NULL;
  int flag = 1;
  VALUE vs, vv, vu;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    switch (argc) {
    case 2:
      CHECK_VECTOR(argv[1]);
      Data_Get_Struct(argv[1], gsl_vector, w);
      flag = 0;
    /* no break, do next */
    case 1:
#ifdef HAVE_NARRAY_H
      if (NA_IsNArray(argv[0]))
        return rb_gsl_linalg_SV_decomp_narray(argc, argv, obj);
#endif

#ifdef HAVE_NMATRIX_H
      if (NM_IsNMatrix(argv[0])) {
        return rb_gsl_linalg_SV_decomp_nmatrix(argc, argv, obj);
      }
#endif
      CHECK_MATRIX(argv[0]);
      Data_Get_Struct(argv[0], gsl_matrix, A);
      break;
    default:
      rb_raise(rb_eArgError, "wrong number of arguments (%d for 1 or 2)", argc);
      break;
    }
    break;
  default:
    switch (argc) {
    case 0:
      /* do nothing */
      break;
    case 1:
      CHECK_VECTOR(argv[0]);
      Data_Get_Struct(argv[0], gsl_vector, w);
      flag = 0;
      break;
    default:
      rb_raise(rb_eArgError, "wrong number of arguments (%d for 0 or 1)", argc);
      break;
    }
    Data_Get_Struct(obj, gsl_matrix, A);
    break;
  }
  U = make_matrix_clone(A);
  S = gsl_vector_alloc(A->size2);   /* see manual p 123 */
  V = gsl_matrix_alloc(A->size2, A->size2);
  if (flag == 1) w = gsl_vector_alloc(A->size2);
  gsl_linalg_SV_decomp(U, V, S, w);
  if (flag == 1) gsl_vector_free(w);
  vu = Data_Wrap_Struct(cgsl_matrix_U, 0, gsl_matrix_free, U);
  vv = Data_Wrap_Struct(cgsl_matrix_V, 0, gsl_matrix_free, V);
  vs = Data_Wrap_Struct(cgsl_vector_S, 0, gsl_vector_free, S);
  return rb_ary_new3(3, vu, vv, vs);
}

static VALUE rb_gsl_linalg_SV_decomp_mod(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *V = NULL, *U = NULL, *X = NULL;
  gsl_vector *w = NULL, *S = NULL;
  VALUE vs, vv, vu;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError,
                            "wrong number of argument (%d for 1)", argc);
    CHECK_MATRIX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    break;
  default:
    Data_Get_Struct(obj, gsl_matrix, A);
    break;
  }
  U = make_matrix_clone(A);
  S = gsl_vector_alloc(A->size2);   /* see manual p 123 */
  V = gsl_matrix_alloc(A->size2, A->size2);
  X = gsl_matrix_alloc(A->size2, A->size2);
  w = gsl_vector_alloc(A->size2);
  gsl_linalg_SV_decomp_mod(U, X, V, S, w);
  gsl_vector_free(w);
  gsl_matrix_free(X);
  vu = Data_Wrap_Struct(cgsl_matrix_U, 0, gsl_matrix_free, U);
  vv = Data_Wrap_Struct(cgsl_matrix_V, 0, gsl_matrix_free, V);
  vs = Data_Wrap_Struct(cgsl_vector_S, 0, gsl_vector_free, S);
  return rb_ary_new3(3, vu, vv, vs);
}

static VALUE rb_gsl_linalg_SV_decomp_jacobi(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *V = NULL, *U = NULL;
  gsl_vector *S = NULL;
  VALUE vs, vv, vu;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError,
                            "wrong number of argument (%d for 1)", argc);
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0]))
      return rb_gsl_linalg_SV_decomp_jacobi_narray(argc, argv, obj);
#endif
    CHECK_MATRIX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    break;
  default:
    Data_Get_Struct(obj, gsl_matrix, A);
    break;
  }
  U = make_matrix_clone(A);
  S = gsl_vector_alloc(A->size2);   /* see manual p 123 */
  V = gsl_matrix_alloc(A->size2, A->size2);
  gsl_linalg_SV_decomp_jacobi(U, V, S);
  vu = Data_Wrap_Struct(cgsl_matrix_U, 0, gsl_matrix_free, U);
  vv = Data_Wrap_Struct(cgsl_matrix_V, 0, gsl_matrix_free, V);
  vs = Data_Wrap_Struct(cgsl_vector_S, 0, gsl_vector_free, S);
  return rb_ary_new3(3, vu, vv, vs);
}

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_SV_solve_nmatrix(int argc, VALUE *argv, VALUE obj)
{
  NM_DENSE_STORAGE *A;
  gsl_matrix_view uv, vv;
  gsl_vector_view sv, bv, xv;
  VALUE x;
  if (argc != 4){
    rb_raise(rb_eArgError, "Usage: SV.solve(u, v, s, b)");
  }
  A = NM_STORAGE_DENSE(argv[0]);
  uv = gsl_matrix_view_array((double*)NM_DENSE_ELEMENTS(argv[0]), A->shape[0], A->shape[1]);
  vv = gsl_matrix_view_array((double*)NM_DENSE_ELEMENTS(argv[1]), A->shape[0], A->shape[0]);
  sv = gsl_vector_view_array((double*)NM_DENSE_ELEMENTS(argv[2]), A->shape[0]);
  bv = gsl_vector_view_array((double*)NM_DENSE_ELEMENTS(argv[3]), A->shape[0]);
  x = rb_nvector_dense_create(FLOAT64, (double*)NM_DENSE_ELEMENTS(argv[3]), A->shape[0]);
  xv = gsl_vector_view_array((double*)NM_DENSE_ELEMENTS(x), A->shape[0]);
  gsl_linalg_SV_solve(&uv.matrix, &vv.matrix, &sv.vector, &bv.vector, &xv.vector);
  return x;
}
#endif

static VALUE rb_gsl_linalg_SV_solve(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *U = NULL, *V = NULL;
  gsl_vector *S = NULL, *b = NULL, *x = NULL;
  int flagb = 0, flagv = 0;

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc < 1) rb_raise(rb_eArgError, "too few arguments.");
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0]))
      return rb_gsl_linalg_SV_solve_narray(argc, argv, obj);
#endif

#ifdef HAVE_NMATRIX_H
    if (NM_IsNMatrix(argv[0]))
      return rb_gsl_linalg_SV_solve_nmatrix(argc, argv, obj);
#endif

    CHECK_MATRIX(argv[0]);
    if (CLASS_OF(argv[0]) == cgsl_matrix_U) {
      if (argc != 4) rb_raise(rb_eArgError,
                              "wrong number of arguments (%d for 4)", argc);
      Data_Get_Struct(argv[0], gsl_matrix, U);
      CHECK_MATRIX(argv[1]);
      if (CLASS_OF(argv[1]) != cgsl_matrix_V)
        rb_raise(rb_eTypeError, "not a V matrix");
      Data_Get_Struct(argv[1], gsl_matrix, V);
      CHECK_VECTOR(argv[2]);
      if (CLASS_OF(argv[2]) != cgsl_vector_S)
        rb_raise(rb_eTypeError, "not a S vector");
      Data_Get_Struct(argv[2], gsl_vector, S);
      if (TYPE(argv[3]) == T_ARRAY) {
        b = make_cvector_from_rarray(argv[3]);
        flagb = 1;
      } else {
        CHECK_VECTOR(argv[3]);
        Data_Get_Struct(argv[3], gsl_vector, b);
      }
    } else {
      if (argc != 2) rb_raise(rb_eArgError,
                              "wrong number of arguments (%d for 2)", argc);
      Data_Get_Struct(argv[0], gsl_matrix, A);
      U = make_matrix_clone(A);
      if (TYPE(argv[1]) == T_ARRAY) {
        b = make_cvector_from_rarray(argv[1]);
        flagb = 1;
      } else {
        CHECK_VECTOR(argv[1]);
        Data_Get_Struct(argv[1], gsl_vector, b);
      }
      S = gsl_vector_alloc(A->size2);   /* see manual p 123 */
      V = gsl_matrix_alloc(A->size2, A->size2);
      gsl_linalg_SV_decomp_jacobi(U, V, S);
      flagv = 1;
    }
    break;
  default:
    if (argc != 1) rb_raise(rb_eArgError,
                            "wrong number of arguments (%d for 1)", argc);
    Data_Get_Struct(obj, gsl_matrix, A);
    U = make_matrix_clone(A);
    if (TYPE(argv[0]) == T_ARRAY) {
      b = make_cvector_from_rarray(argv[0]);
      flagb = 1;
    } else {
      CHECK_VECTOR(argv[0]);
      Data_Get_Struct(argv[0], gsl_vector, b);
    }
    S = gsl_vector_alloc(A->size2);   /* see manual p 123 */
    V = gsl_matrix_alloc(A->size2, A->size2);
    gsl_linalg_SV_decomp_jacobi(U, V, S);
    flagv = 1;
    break;
  }
  //  x = gsl_vector_alloc(b->size);
  // Bug report #25842
  x = gsl_vector_alloc(S->size);
  gsl_linalg_SV_solve(U, V, S, b, x);
  if (flagv == 1) {
    gsl_matrix_free(U);
    gsl_matrix_free(V);
    gsl_vector_free(S);
  }
  if (flagb == 1) gsl_vector_free(b);
  return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
}

/*****/

#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_linalg_cholesky_decomp_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *na;
  VALUE chol;
  gsl_matrix_view mv;
  GetNArray(argv[0], na);
  chol = na_make_object(NA_DFLOAT, 2, na->shape, CLASS_OF(argv[0]));
  memcpy(NA_PTR_TYPE(chol,double*), (double*)na->ptr, sizeof(double)*na->total);
  mv = gsl_matrix_view_array(NA_PTR_TYPE(chol,double*), na->shape[1], na->shape[0]);
  gsl_linalg_cholesky_decomp(&mv.matrix);
  return chol;
}

static VALUE rb_gsl_linalg_cholesky_solve_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *nm, *nb;
  VALUE x;
  gsl_matrix_view mv;
  gsl_vector_view bv, xv;
  switch (argc) {
  case 2:
    GetNArray(argv[0], nm);
    GetNArray(argv[1], nb);
    x = na_make_object(NA_DFLOAT, 1, nb->shape, CLASS_OF(argv[1]));
    break;
  case 3:
    GetNArray(argv[0], nm);
    GetNArray(argv[1], nb);
    x = argv[2];
    break;
  default:
    rb_raise(rb_eArgError,
             "Usage: Cholesky.solve(chol, b) or Cholesky.solve(chol, b, x)");
    break;
  }
  mv = gsl_matrix_view_array((double*)nm->ptr, nm->shape[1], nm->shape[0]);
  bv = gsl_vector_view_array((double*)nb->ptr, nb->shape[0]);
  xv = gsl_vector_view_array(NA_PTR_TYPE(x,double*), nb->shape[0]);
  gsl_linalg_cholesky_solve(&mv.matrix, &bv.vector, &xv.vector);
  return x;
}

static VALUE rb_gsl_linalg_cholesky_svx_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *nm, *nb;
  gsl_matrix_view mv;
  gsl_vector_view bv;
  GetNArray(argv[0], nm);  GetNArray(argv[1], nb);
  mv = gsl_matrix_view_array((double*)nm->ptr, nm->shape[1], nm->shape[0]);
  bv = gsl_vector_view_array((double*)nb->ptr, nb->shape[0]);
  gsl_linalg_cholesky_svx(&mv.matrix, &bv.vector);
  return argv[1];
}
#endif

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_cholesky_decomp_nmatrix(int argc, VALUE *argv, VALUE obj)
{
  NM_DENSE_STORAGE *nm;
  VALUE chol;
  gsl_matrix_view mv;
  nm = NM_STORAGE_DENSE(argv[0]);

  chol = rb_nmatrix_dense_create(FLOAT64, nm->shape, 2, nm->elements,
    nm->shape[0]*nm->shape[1]);
  mv   = gsl_matrix_view_array((double*)NM_DENSE_ELEMENTS(chol), nm->shape[0],
    nm->shape[1]);
  gsl_linalg_cholesky_decomp(&mv.matrix);
  return chol;
}
#endif

static VALUE rb_gsl_linalg_cholesky_decomp(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *Atmp = NULL;
  switch(TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0]))
      return rb_gsl_linalg_cholesky_decomp_narray(argc, argv, obj);
#endif

#ifdef HAVE_NMATRIX_H
    if (NM_IsNMatrix(argv[0]))
      return rb_gsl_linalg_cholesky_decomp_nmatrix(argc, argv, obj);
#endif
    CHECK_MATRIX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix, Atmp);
    break;
  default:
    CHECK_MATRIX(obj);
    Data_Get_Struct(obj, gsl_matrix, Atmp);
    break;
  }
  A = make_matrix_clone(Atmp);
  gsl_linalg_cholesky_decomp(A);
  return Data_Wrap_Struct(cgsl_matrix_C, 0, gsl_matrix_free, A);
}

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_cholesky_solve_nmatrix(int argc, VALUE *argv, VALUE obj)
{
  NM_DENSE_STORAGE *nm, *nb;
  VALUE x;
  gsl_matrix_view mv;
  gsl_vector_view bv, xv;

  nm = NM_STORAGE_DENSE(argv[0]);
  nb = NM_STORAGE_DENSE(argv[1]);
  switch (argc) {
  case 2:
    x = rb_nvector_dense_create(FLOAT64, nb->elements, nb->shape[0]);
    break;
  case 3:
    x = argv[2];
    break;
  default:
    rb_raise(rb_eArgError,
             "Usage: Cholesky.solve(chol, b) or Cholesky.solve(chol, b, x)");
    break;
  }
  mv = gsl_matrix_view_array((double*)nm->elements, nm->shape[0], nm->shape[1]);
  bv = gsl_vector_view_array((double*)nb->elements, nb->shape[0]);
  xv = gsl_vector_view_array((double*)NM_DENSE_ELEMENTS(x), nb->shape[0]);
  gsl_linalg_cholesky_solve(&mv.matrix, &bv.vector, &xv.vector);
  return x;
}
#endif

static VALUE rb_gsl_linalg_cholesky_solve(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *Atmp = NULL;
  gsl_vector *b = NULL, *x = NULL;
  int flagb = 0, flaga = 0;
  VALUE vA, vb;
  switch(TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 2) rb_raise(rb_eArgError, "wrong number of argument (%d for 2)",
                            argc);
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0]))
      return rb_gsl_linalg_cholesky_solve_narray(argc, argv, obj);
#endif

#ifdef HAVE_NMATRIX_H
    if (NM_IsNMatrix(argv[0]))
      return rb_gsl_linalg_cholesky_solve_nmatrix(argc, argv, obj);
#endif
    vA = argv[0];
    vb = argv[1];
    break;
  default:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    vA = obj;
    vb = argv[0];
    break;
  }
  CHECK_MATRIX(vA);
  Data_Get_Struct(vA, gsl_matrix, Atmp);
  if (TYPE(vb) == T_ARRAY) {
    b = make_cvector_from_rarray(vb);
    flagb = 1;
  } else {
    CHECK_VECTOR(vb);
    Data_Get_Struct(vb, gsl_vector, b);
  }
  if (CLASS_OF(vA) == cgsl_matrix_C) {
    A = Atmp;
  } else {
    A = make_matrix_clone(Atmp);
    flaga = 1;
    gsl_linalg_cholesky_decomp(A);
  }
  x = gsl_vector_alloc(b->size);
  gsl_linalg_cholesky_solve(A, b, x);
  if (flaga == 1) gsl_matrix_free(A);
  if (flagb == 1) gsl_vector_free(b);
  return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
}

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_cholesky_svx_nmatrix(int argc, VALUE *argv, VALUE obj)
{
  NM_DENSE_STORAGE *nm, *nb;
  gsl_matrix_view mv;
  gsl_vector_view bv;
  nm = NM_STORAGE_DENSE(argv[0]);
  nb = NM_STORAGE_DENSE(argv[1]);

  mv = gsl_matrix_view_array((double*)nm->elements, nm->shape[0], nm->shape[1]);
  bv = gsl_vector_view_array((double*)nb->elements, nb->shape[0]);
  gsl_linalg_cholesky_svx(&mv.matrix, &bv.vector);

  return argv[1];
}
#endif

static VALUE rb_gsl_linalg_cholesky_svx(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *Atmp = NULL;
  gsl_vector *b = NULL;
  int flaga = 0;
  VALUE vA, vb;
  switch(TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 2) rb_raise(rb_eArgError, "wrong number of argument (%d for 2)",
                            argc);
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0]))
      return rb_gsl_linalg_cholesky_svx_narray(argc, argv, obj);
#endif

#ifdef HAVE_NMATRIX_H
    if (NM_IsNMatrix(argv[0]))
      return rb_gsl_linalg_cholesky_svx_nmatrix(argc, argv, obj);
#endif
    vA = argv[0];
    vb = argv[1];
    break;
  default:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    vA = obj;
    vb = argv[0];
    break;
  }
  CHECK_MATRIX(vA);
  Data_Get_Struct(vA, gsl_matrix, Atmp);
  CHECK_VECTOR(vb);
  Data_Get_Struct(vb, gsl_vector, b);
  if (CLASS_OF(vA) == cgsl_matrix_C) {
    A = Atmp;
  } else {
    A = make_matrix_clone(Atmp);
    flaga = 1;
    gsl_linalg_cholesky_decomp(A);
  }
  gsl_linalg_cholesky_svx(A, b);
  if (flaga == 1) gsl_matrix_free(A);
  return vb;
}

static VALUE rb_gsl_linalg_symmtd_decomp(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *Atmp = NULL;
  gsl_vector *tau = NULL;
  VALUE vQ, vtau;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    CHECK_MATRIX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix, Atmp);
    break;
  default:
    CHECK_MATRIX(obj);
    Data_Get_Struct(obj, gsl_matrix, Atmp);
    break;
  }
  A = make_matrix_clone(Atmp);
  tau = gsl_vector_alloc(A->size1);
  gsl_linalg_symmtd_decomp(A, tau);
  vQ = Data_Wrap_Struct(cgsl_matrix_Q, 0, gsl_matrix_free, A);
  vtau = Data_Wrap_Struct(cgsl_vector_tau, 0, gsl_vector_free, tau);
  return rb_ary_new3(2, vQ, vtau);
}


static VALUE rb_gsl_linalg_symmtd_decomp2(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL;
  gsl_vector *tau = NULL;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    CHECK_MATRIX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    break;
  default:
    CHECK_MATRIX(obj);
    Data_Get_Struct(obj, gsl_matrix, A);
    break;
  }
  tau = gsl_vector_alloc(A->size1);
  gsl_linalg_symmtd_decomp(A, tau);
  return Data_Wrap_Struct(cgsl_vector_tau, 0, gsl_vector_free, tau);
}

static VALUE rb_gsl_linalg_symmtd_unpack(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *Q = NULL;
  gsl_vector *tau = NULL, *d = NULL, *sd = NULL;
  VALUE vq, vd, vsd;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 2) rb_raise(rb_eArgError, "wrong number of argument (%d for 2)",
                            argc);
    CHECK_MATRIX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    Data_Get_Struct(argv[1], gsl_vector, tau);
    break;
  default:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    CHECK_MATRIX(obj);
    Data_Get_Struct(obj, gsl_matrix, A);
    Data_Get_Struct(argv[0], gsl_vector, tau);
    break;
  }
  Q = gsl_matrix_alloc(A->size1, A->size2);
  d = gsl_vector_alloc(tau->size);
  sd = gsl_vector_alloc(tau->size);
  gsl_linalg_symmtd_unpack(A, tau, Q, d, sd);

  vq = Data_Wrap_Struct(cgsl_matrix_Q, 0, gsl_matrix_free, Q);
  vd = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, d);
  vsd = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, sd);

  return rb_ary_new3(3, vq, vd, vsd);
}

static VALUE rb_gsl_linalg_symmtd_unpack_T(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL;
  gsl_vector *d = NULL, *sd = NULL;
  VALUE vd, vsd;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 2)",
                            argc);
    CHECK_MATRIX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    break;
  default:
    Data_Get_Struct(obj, gsl_matrix, A);
    break;
  }
  d = gsl_vector_alloc(A->size1);
  sd = gsl_vector_alloc(A->size1);
  gsl_linalg_symmtd_unpack_T(A, d, sd);

  vd = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, d);
  vsd = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, sd);

  return rb_ary_new3(2, vd, vsd);
}

/*****/

static VALUE rb_gsl_linalg_hermtd_decomp(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix_complex *A = NULL, *Atmp = NULL;
  gsl_vector_complex *tau = NULL;
  VALUE vQ, vtau;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    CHECK_MATRIX_COMPLEX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix_complex, Atmp);
    break;
  default:
    CHECK_MATRIX_COMPLEX(obj);
    Data_Get_Struct(obj, gsl_matrix_complex, Atmp);
    break;
  }
  A = make_matrix_complex_clone(Atmp);
  tau = gsl_vector_complex_alloc(A->size1);
  gsl_linalg_hermtd_decomp(A, tau);
  vQ = Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free, A);
  vtau = Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, tau);
  return rb_ary_new3(2, vQ, vtau);
}

static VALUE rb_gsl_linalg_hermtd_decomp2(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix_complex *A = NULL;
  gsl_vector_complex *tau = NULL;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    CHECK_MATRIX_COMPLEX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix_complex, A);
    break;
  default:
    CHECK_MATRIX_COMPLEX(obj);
    Data_Get_Struct(obj, gsl_matrix_complex, A);
    break;
  }
  tau = gsl_vector_complex_alloc(A->size1);
  gsl_linalg_hermtd_decomp(A, tau);
  return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, tau);
}

static VALUE rb_gsl_linalg_hermtd_unpack(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix_complex *A = NULL, *Q = NULL;
  gsl_vector_complex *tau = NULL;
  gsl_vector *d = NULL, *sd = NULL;
  VALUE vq, vd, vsd;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 2) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    CHECK_MATRIX_COMPLEX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix_complex, A);
    Data_Get_Struct(argv[1], gsl_vector_complex, tau);
    break;
  default:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    CHECK_MATRIX_COMPLEX(obj);
    Data_Get_Struct(obj, gsl_matrix_complex, A);
    Data_Get_Struct(argv[0], gsl_vector_complex, tau);
    break;
  }
  Q = gsl_matrix_complex_alloc(A->size1, A->size2);
  d = gsl_vector_alloc(tau->size);
  sd = gsl_vector_alloc(tau->size);
  gsl_linalg_hermtd_unpack(A, tau, Q, d, sd);

  vq = Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free, Q);
  vd = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, d);
  vsd = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, sd);

  return rb_ary_new3(3, vq, vd, vsd);
}

static VALUE rb_gsl_linalg_hermtd_unpack_T(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix_complex *A = NULL;
  gsl_vector *d = NULL, *sd = NULL;
  VALUE vd, vsd;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    CHECK_MATRIX_COMPLEX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix_complex, A);
    break;
  default:
    Data_Get_Struct(obj, gsl_matrix_complex, A);
    break;
  }
  d = gsl_vector_alloc(A->size1);
  sd = gsl_vector_alloc(A->size1);
  gsl_linalg_hermtd_unpack_T(A, d, sd);

  vd = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, d);
  vsd = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, sd);

  return rb_ary_new3(2, vd, vsd);
}

/******/

static VALUE rb_gsl_linalg_bidiag_decomp(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *Atmp = NULL;
  gsl_vector *tau_U = NULL, *tau_V = NULL;
  size_t size0;
  // local variable "status" was defined and set, but never used
  //int status;
  VALUE vu, vv, vA;

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of arguments (%d for 1)",
                            argc);
    Data_Get_Struct(argv[0], gsl_matrix, Atmp);
    break;
  default:
    Data_Get_Struct(obj, gsl_matrix, Atmp);
    break;
  }
  A = make_matrix_clone(Atmp);
  size0 = GSL_MIN(A->size1, A->size2);
  tau_U = gsl_vector_alloc(size0);
  tau_V = gsl_vector_alloc(size0-1);
  /*status =*/ gsl_linalg_bidiag_decomp(A, tau_U, tau_V);
  vA = Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, A);
  vu = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, tau_U);
  vv = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, tau_V);
  return rb_ary_new3(3, vA, vu, vv);
}

static VALUE rb_gsl_linalg_bidiag_decomp2(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL;
  gsl_vector *tau_U = NULL, *tau_V = NULL;
  size_t size0;
  VALUE vu, vv;

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of arguments (%d for 1)",
                            argc);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    break;
  default:
    Data_Get_Struct(obj, gsl_matrix, A);
    break;
  }
  size0 = GSL_MIN(A->size1, A->size2);
  tau_U = gsl_vector_alloc(size0);
  tau_V = gsl_vector_alloc(size0-1);
  gsl_linalg_bidiag_decomp(A, tau_U, tau_V);
  vu = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, tau_U);
  vv = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, tau_V);
  return rb_ary_new3(2, vu, vv);
}

static VALUE rb_gsl_linalg_bidiag_unpack(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *U = NULL, *V = NULL;
  gsl_vector *tau_U = NULL, *tau_V = NULL, *d = NULL, *s = NULL;
  size_t size0;
  VALUE vu, vv, vd, vs;

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 3) rb_raise(rb_eArgError, "wrong number of arguments (%d for 3)",
                            argc);
    CHECK_MATRIX(argv[0]);
    CHECK_VECTOR(argv[1]);
    CHECK_VECTOR(argv[2]);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    Data_Get_Struct(argv[1], gsl_vector, tau_U);
    Data_Get_Struct(argv[2], gsl_vector, tau_V);
    break;
  default:
    if (argc != 2) rb_raise(rb_eArgError, "wrong number of arguments (%d for 2)",
                            argc);
    CHECK_MATRIX(obj);
    CHECK_VECTOR(argv[0]);
    CHECK_VECTOR(argv[1]);
    Data_Get_Struct(obj, gsl_matrix, A);
    Data_Get_Struct(argv[0], gsl_vector, tau_U);
    Data_Get_Struct(argv[1], gsl_vector, tau_V);
    break;
  }
  size0 = GSL_MIN(A->size1, A->size2);
  U = gsl_matrix_alloc(A->size1, A->size2);
  V = gsl_matrix_alloc(size0, size0);

  d = gsl_vector_alloc(size0);
  s = gsl_vector_alloc(size0-1);
  gsl_linalg_bidiag_unpack(A, tau_U, U, tau_V, V, d, s);
  vu = Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, U);
  vv = Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, V);
  vd = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, d);
  vs = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, s);
  return rb_ary_new3(4, vu, vv, vd, vs);
}

static VALUE rb_gsl_linalg_bidiag_unpack2(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *V = NULL;
  gsl_vector *tau_V = NULL, *tau_U = NULL;
  VALUE vv;

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 3) rb_raise(rb_eArgError, "wrong number of arguments (%d for 3)",
                            argc);
    CHECK_MATRIX(argv[0]);
    CHECK_VECTOR(argv[1]);
    CHECK_VECTOR(argv[2]);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    Data_Get_Struct(argv[1], gsl_vector, tau_U);
    Data_Get_Struct(argv[2], gsl_vector, tau_V);
    break;
  default:
    if (argc != 2) rb_raise(rb_eArgError, "wrong number of arguments (%d for 2)",
                            argc);
    CHECK_MATRIX(obj);
    CHECK_VECTOR(argv[0]);
    CHECK_VECTOR(argv[1]);
    Data_Get_Struct(obj, gsl_matrix, A);
    Data_Get_Struct(argv[0], gsl_vector, tau_U);
    Data_Get_Struct(argv[1], gsl_vector, tau_V);
    break;
  }
  V = gsl_matrix_alloc(A->size2, A->size2);
  gsl_linalg_bidiag_unpack2(A, tau_U, tau_V, V);
  vv = Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, V);
  return vv;
}

static VALUE rb_gsl_linalg_bidiag_unpack_B(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL;
  gsl_vector *d = NULL, *s = NULL;
  size_t size0;
  VALUE vd, vs;

  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of arguments (%d for 3)",
                            argc);
    CHECK_MATRIX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    break;
  default:
    CHECK_MATRIX(obj);
    Data_Get_Struct(obj, gsl_matrix, A);
    break;
  }
  size0 = GSL_MIN(A->size1, A->size2);
  d = gsl_vector_alloc(size0);
  s = gsl_vector_alloc(size0);
  gsl_linalg_bidiag_unpack_B(A, d, s);
  vd = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, d);
  vs = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, s);
  return rb_ary_new3(2, vd, vs);
}

/* Householder Transformations 11.Jul.2004 */
static VALUE rb_gsl_linalg_householder_transform(int argc, VALUE *argv, VALUE obj)
{
  gsl_vector *v = NULL;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc < 1) rb_raise(rb_eArgError, "too few arguments.");
    CHECK_VECTOR(argv[0]);
    Data_Get_Struct(argv[0], gsl_vector, v);
    break;
  default:
    Data_Get_Struct(obj, gsl_vector, v);
    break;
  }
  return rb_float_new(gsl_linalg_householder_transform(v));
}

/* singleton */
static VALUE rb_gsl_linalg_householder_hm(VALUE obj, VALUE t, VALUE vv, VALUE aa)
{
  gsl_vector *v = NULL;
  double tau;
  gsl_matrix *A = NULL;
  CHECK_VECTOR(vv);
  CHECK_MATRIX(aa);
  tau = NUM2DBL(t);
  Data_Get_Struct(vv, gsl_vector, v);
  Data_Get_Struct(aa, gsl_matrix, A);
  gsl_linalg_householder_hm(tau, v, A);
  return aa;
}

static VALUE rb_gsl_linalg_householder_mh(VALUE obj, VALUE t, VALUE vv, VALUE aa)
{
  gsl_vector *v = NULL;
  double tau;
  gsl_matrix *A = NULL;
  CHECK_VECTOR(vv);
  CHECK_MATRIX(aa);
  tau = NUM2DBL(t);
  Data_Get_Struct(vv, gsl_vector, v);
  Data_Get_Struct(aa, gsl_matrix, A);
  gsl_linalg_householder_mh(tau, v, A);
  return aa;
}

static VALUE rb_gsl_linalg_householder_hv(VALUE obj, VALUE t, VALUE vv, VALUE ww)
{
  gsl_vector *v = NULL, *w = NULL;
  double tau;
  CHECK_VECTOR(vv);
  CHECK_VECTOR(ww);
  tau = NUM2DBL(t);
  Data_Get_Struct(vv, gsl_vector, v);
  Data_Get_Struct(ww, gsl_vector, w);
  gsl_linalg_householder_hv(tau, v, w);
  return ww;
}

#ifdef HAVE_NARRAY_H
static VALUE rb_gsl_linalg_HH_solve_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *na;
  gsl_vector_view bv, xv;
  VALUE x;
  gsl_matrix *mtmp;
  GetNArray(argv[0], na);
  bv = gsl_vector_view_array(NA_PTR_TYPE(argv[1],double*), na->shape[1]);
  x = na_make_object(NA_DFLOAT, 1, &na->shape[1], CLASS_OF(argv[1]));
  xv = gsl_vector_view_array(NA_PTR_TYPE(x,double*), na->shape[1]);
  mtmp = gsl_matrix_alloc(na->shape[1], na->shape[0]);
  memcpy(mtmp->data, (double*)na->ptr, sizeof(double)*na->total);
  gsl_linalg_HH_solve(mtmp, &bv.vector, &xv.vector);
  gsl_matrix_free(mtmp);
  return x;
}
static VALUE rb_gsl_linalg_HH_svx_narray(int argc, VALUE *argv, VALUE obj)
{
  struct NARRAY *na;
  gsl_matrix *mtmp;
  gsl_vector_view bv;
  GetNArray(argv[0], na);
  bv = gsl_vector_view_array(NA_PTR_TYPE(argv[1],double*), na->shape[1]);
  mtmp = gsl_matrix_alloc(na->shape[1], na->shape[0]);
  memcpy(mtmp->data, (double*)na->ptr, sizeof(double)*na->total);
  gsl_linalg_HH_svx(mtmp, &bv.vector);
  gsl_matrix_free(mtmp);
  return argv[1];
}
#endif

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_HH_solve_nmatrix(int argc, VALUE *argv, VALUE obj)
{
  NM_DENSE_STORAGE *nm;
  gsl_vector_view bv, xv;
  VALUE x;
  gsl_matrix *mtmp;

  nm = NM_STORAGE_DENSE(argv[0]);
  bv = gsl_vector_view_array((double*)NM_DENSE_ELEMENTS(argv[1]), nm->shape[0]);
  x = rb_nvector_dense_create(FLOAT64, nm->elements, nm->shape[0]);
  xv = gsl_vector_view_array((double*)NM_DENSE_ELEMENTS(x), nm->shape[0]);
  mtmp = gsl_matrix_alloc(nm->shape[0], nm->shape[1]);
  memcpy(mtmp->data, (double*)nm->elements, sizeof(double)*nm->shape[0]*nm->shape[1]);
  gsl_linalg_HH_solve(mtmp, &bv.vector, &xv.vector);
  gsl_matrix_free(mtmp);
  return x;
}
#endif

/* 17.Apr.2004 */
static VALUE rb_gsl_linalg_HH_solve(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *Atmp = NULL;
  gsl_vector *b = NULL, *x = NULL;
  int flagb = 0;
  VALUE vA, vb;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 2) rb_raise(rb_eArgError, "wrong number of argument (%d for 2)",
                            argc);
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0]))
      return rb_gsl_linalg_HH_solve_narray(argc, argv, obj);
#endif

#ifdef HAVE_NMATRIX_H
    if (NM_IsNMatrix(argv[0]))
      return rb_gsl_linalg_HH_solve_nmatrix(argc, argv, obj);
#endif
    vA = argv[0];
    vb = argv[1];
    break;
  default:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    vA = obj;
    vb = argv[0];
    break;
  }
  CHECK_MATRIX(vA);
  Data_Get_Struct(vA, gsl_matrix, Atmp);
  if (TYPE(vb) == T_ARRAY) {
    b = make_cvector_from_rarray(vb);
    flagb = 1;
  } else {
    CHECK_VECTOR(vb);
    Data_Get_Struct(vb, gsl_vector, b);
  }
  A = make_matrix_clone(Atmp);
  x = gsl_vector_alloc(b->size);
  gsl_linalg_HH_solve(A, b, x);
  gsl_matrix_free(A);
  if (flagb == 1) gsl_vector_free(b);
  return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
}

static VALUE rb_gsl_linalg_HH_solve_bang(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL;
  gsl_vector *b = NULL, *x = NULL;
  int flagb = 0;
  VALUE vA, vb;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 2) rb_raise(rb_eArgError, "wrong number of argument (%d for 2)",
                            argc);
    vA = argv[0];
    vb = argv[1];
    break;
  default:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    vA = obj;
    vb = argv[0];
    break;
  }
  CHECK_MATRIX(vA);
  Data_Get_Struct(vA, gsl_matrix, A);
  if (TYPE(vb) == T_ARRAY) {
    b = make_cvector_from_rarray(vb);
    flagb = 1;
  } else {
    CHECK_VECTOR(vb);
    Data_Get_Struct(vb, gsl_vector, b);
  }
  x = gsl_vector_alloc(b->size);
  gsl_linalg_HH_solve(A, b, x);
  if (flagb == 1) gsl_vector_free(b);
  return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
}

#ifdef HAVE_NMATRIX_H
static VALUE rb_gsl_linalg_HH_svx_nmatrix(int argc, VALUE *argv, VALUE obj)
{
  NM_DENSE_STORAGE *nm;
  gsl_matrix *mtmp;
  gsl_vector_view bv;

  nm = NM_STORAGE_DENSE(argv[0]);
  bv = gsl_vector_view_array((double*)NM_DENSE_ELEMENTS(argv[1]), nm->shape[0]);
  mtmp = gsl_matrix_alloc(nm->shape[0], nm->shape[1]);
  memcpy(mtmp->data, (double*)nm->elements, sizeof(double)*nm->shape[0]*nm->shape[1]);
  gsl_linalg_HH_svx(mtmp, &bv.vector);
  gsl_matrix_free(mtmp);
  return argv[1];
}
#endif

static VALUE rb_gsl_linalg_HH_svx(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *Atmp = NULL;
  gsl_vector *b = NULL;
  VALUE vA, vb;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 2) rb_raise(rb_eArgError, "wrong number of argument (%d for 2)",
                            argc);
#ifdef HAVE_NARRAY_H
    if (NA_IsNArray(argv[0]))
      return rb_gsl_linalg_HH_svx_narray(argc, argv, obj);
#endif

#ifdef HAVE_NMATRIX_H
    if (NM_IsNMatrix(argv[0]))
      return rb_gsl_linalg_HH_svx_nmatrix(argc, argv, obj);
#endif
    vA = argv[0];
    vb = argv[1];
    break;
  default:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
                            argc);
    vA = obj;
    vb = argv[0];
    break;
  }
  CHECK_MATRIX(vA);
  Data_Get_Struct(vA, gsl_matrix, Atmp);
  CHECK_VECTOR(vb);
  Data_Get_Struct(vb, gsl_vector, b);
  A = make_matrix_clone(Atmp);
  gsl_linalg_HH_svx(A, b);
  gsl_matrix_free(A);
  return vb;
}

static VALUE rb_gsl_linalg_solve_symm_tridiag(VALUE obj, VALUE dd, VALUE ee, VALUE bb)
{
  gsl_vector *b = NULL, *x = NULL, *d = NULL, *e = NULL;

  Data_Get_Struct(dd, gsl_vector, d);
  Data_Get_Struct(ee, gsl_vector, e);
  Data_Get_Struct(bb, gsl_vector, b);
  x = gsl_vector_alloc(b->size);

  gsl_linalg_solve_symm_tridiag(d, e, b, x);
  return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
}

static VALUE rb_gsl_linalg_solve_tridiag(VALUE obj, VALUE dd, VALUE ee, VALUE ff,
                                         VALUE bb)
{
  gsl_vector *b = NULL, *x = NULL, *d = NULL, *e = NULL, *f = NULL;

  Data_Get_Struct(dd, gsl_vector, d);
  Data_Get_Struct(ee, gsl_vector, e);
  Data_Get_Struct(ff, gsl_vector, f);
  Data_Get_Struct(bb, gsl_vector, b);
  x = gsl_vector_alloc(b->size);

  gsl_linalg_solve_tridiag(d, e, f, b, x);
  return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
}

static VALUE rb_gsl_linalg_solve_symm_cyc_tridiag(VALUE obj, VALUE dd, VALUE ee, VALUE bb)
{
  gsl_vector *b = NULL, *x = NULL, *d = NULL, *e = NULL;

  Data_Get_Struct(dd, gsl_vector, d);
  Data_Get_Struct(ee, gsl_vector, e);
  Data_Get_Struct(bb, gsl_vector, b);
  x = gsl_vector_alloc(b->size);

  gsl_linalg_solve_symm_cyc_tridiag(d, e, b, x);
  return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
}

static VALUE rb_gsl_linalg_solve_cyc_tridiag(VALUE obj, VALUE dd, VALUE ee,
                                             VALUE ff, VALUE bb)
{
  gsl_vector *b = NULL, *x = NULL, *d = NULL, *e = NULL, *f = NULL;
  Data_Get_Struct(dd, gsl_vector, d);
  Data_Get_Struct(ee, gsl_vector, e);
  Data_Get_Struct(ff, gsl_vector, f);
  Data_Get_Struct(bb, gsl_vector, b);
  x = gsl_vector_alloc(b->size);
  gsl_linalg_solve_cyc_tridiag(d, e, f, b, x);
  return Data_Wrap_Struct(cgsl_vector_col, 0, gsl_vector_free, x);
}

static void rb_gsl_linalg_balance_columns_init(int argc, VALUE *argv, VALUE obj,
                                               VALUE *mat, VALUE *vec,
                                               gsl_matrix **M, gsl_vector **V)
{
  gsl_matrix *A = NULL;
  gsl_vector *D = NULL;
  switch (TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    switch (argc) {
    case 2:
      CHECK_MATRIX(argv[0]); CHECK_VECTOR(argv[1]);
      Data_Get_Struct(argv[0], gsl_matrix, A);
      Data_Get_Struct(argv[1], gsl_vector, D);
      *vec = argv[1];
      break;
    case 1:
      CHECK_MATRIX(argv[0]);
      Data_Get_Struct(argv[0], gsl_matrix, A);
      D = gsl_vector_alloc(A->size2);
      *vec = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, D);
      break;
    default:
      rb_raise(rb_eArgError, "wrong number of arguments (%d for 1 or 2)", argc);
      break;
    }
    *mat = argv[0];
    break;
  default:
    Data_Get_Struct(obj, gsl_matrix, A);
    switch (argc) {
    case 1:
      CHECK_VECTOR(argv[0]);
      Data_Get_Struct(argv[0], gsl_vector, D);
      *vec = argv[0];
      break;
    case 0:
      D = gsl_vector_alloc(A->size2);
      *vec = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, D);
      break;
    default:
      rb_raise(rb_eArgError, "wrong number of arguments (%d for 0 or 1)", argc);
      break;
    }
    *mat = obj;
    break;
  }
  *M = A;
  *V = D;
}

static VALUE rb_gsl_linalg_balance_columns_bang(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL;
  gsl_vector *D = NULL;
  VALUE mat, vec;
  // local variable "status" was defined and set, but never used
  //int status;
  rb_gsl_linalg_balance_columns_init(argc, argv, obj, &mat, &vec, &A, &D);
  /*status =*/ gsl_linalg_balance_columns(A, D);
  return rb_ary_new3(2, mat, vec);
}

static VALUE rb_gsl_linalg_balance_columns(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix *A = NULL, *Anew;
  gsl_vector *D = NULL;
  VALUE mat, vec;
  // local variable "status" was defined and set, but never used
  //int status;
  rb_gsl_linalg_balance_columns_init(argc, argv, obj, &mat, &vec, &A, &D);
  Anew = make_matrix_clone(A);
  mat = Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, Anew);
  /*status =*/ gsl_linalg_balance_columns(Anew, D);
  return rb_ary_new3(2, mat, vec);
}

static VALUE rb_gsl_linalg_hessenberg_decomp(VALUE module, VALUE AA)
{
  gsl_matrix *A = NULL, *Atmp = NULL;
  gsl_vector *tau = NULL;
  VALUE vH, vtau;
  CHECK_MATRIX(AA);
  Data_Get_Struct(AA, gsl_matrix, Atmp);
  A = make_matrix_clone(Atmp);
  tau = gsl_vector_alloc(A->size1);
  gsl_linalg_hessenberg_decomp(A, tau);
  vH = Data_Wrap_Struct(cgsl_matrix_Q, 0, gsl_matrix_free, A);
  vtau = Data_Wrap_Struct(cgsl_vector_tau, 0, gsl_vector_free, tau);
  return rb_ary_new3(2, vH, vtau);
}

static VALUE rb_gsl_linalg_hessenberg_unpack(VALUE module, VALUE HH, VALUE tt)
{
  gsl_matrix *H = NULL, *U = NULL;
  gsl_vector *tau = NULL;
  CHECK_MATRIX(HH);
  CHECK_VECTOR(tt);
  Data_Get_Struct(HH, gsl_matrix, H);
  Data_Get_Struct(tt, gsl_vector, tau);
  U = gsl_matrix_alloc(H->size1, H->size2);
  gsl_linalg_hessenberg_unpack(H, tau, U);

  return Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, U);
}

static VALUE rb_gsl_linalg_hessenberg_unpack_accum(int argc, VALUE *argv, VALUE module)
{
  gsl_matrix *H = NULL, *V = NULL;
  gsl_vector *tau = NULL;
  size_t i;
  VALUE val = Qnil;
  switch (argc) {
  case 2:
    /* nothing to do */
    break;
  case 3:
    CHECK_MATRIX(argv[2]);
    Data_Get_Struct(argv[2], gsl_matrix, V);
    val = argv[2];
    break;
  default:
    rb_raise(rb_eArgError, "Wrong number of arguments (%d for 2 or 3)", argc);
  }
  CHECK_MATRIX(argv[0]);
  CHECK_VECTOR(argv[1]);
  Data_Get_Struct(argv[0], gsl_matrix, H);
  Data_Get_Struct(argv[1], gsl_vector, tau);
  if (argc == 2) {
    V = gsl_matrix_alloc(H->size1, H->size2);
    val = Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, V);
    for (i = 0; i < V->size1; i++) gsl_matrix_set(V, i, i, 1.0);
  }
  gsl_linalg_hessenberg_unpack_accum(H, tau, V);
  return val;
}
static VALUE rb_gsl_linalg_hessenberg_set_zero(VALUE module, VALUE HH)
{
  gsl_matrix *H;
  CHECK_MATRIX(HH);
  Data_Get_Struct(HH, gsl_matrix, H);
  return INT2FIX(gsl_linalg_hessenberg_set_zero(H));
  /*  gsl_linalg_hessenberg_set_zero(H);
      return INT2FIX(0);*/
}
static VALUE rb_gsl_linalg_hesstri_decomp(int argc, VALUE *argv, VALUE module)
{
  gsl_matrix *A = NULL, *B = NULL, *Anew, *Bnew;
  gsl_matrix *U = NULL, *V = NULL;
  gsl_vector *work = NULL;
  VALUE vH, vR, vU = Qnil, vV = Qnil, ary;
  int flag = 0;
  switch (argc) {
  case 2:
    flag = 1;
    break;
  case 3:
    CHECK_VECTOR(argv[2]);
    Data_Get_Struct(argv[2], gsl_vector, work);
    break;
  case 4:
    CHECK_MATRIX(argv[2]);
    CHECK_MATRIX(argv[3]);
    Data_Get_Struct(argv[2], gsl_matrix, U);
    Data_Get_Struct(argv[3], gsl_matrix, V);
    flag = 1;
    break;
  case 5:
    CHECK_MATRIX(argv[2]);
    CHECK_MATRIX(argv[3]);
    CHECK_VECTOR(argv[4]);
    Data_Get_Struct(argv[2], gsl_matrix, U);
    Data_Get_Struct(argv[3], gsl_matrix, V);
    Data_Get_Struct(argv[4], gsl_vector, work);
    vU = argv[2];
    vV = argv[3];
    break;
  default:
    rb_raise(rb_eArgError, "Wrong number of arguments (%d for 2-55)", argc);
  }
  CHECK_MATRIX(argv[0]);
  CHECK_MATRIX(argv[1]);
  Data_Get_Struct(argv[0], gsl_matrix, A);
  Data_Get_Struct(argv[1], gsl_matrix, B);
  Anew = make_matrix_clone(A);
  Bnew = make_matrix_clone(B);
  if (flag == 1) work = gsl_vector_alloc(A->size1);
  gsl_linalg_hesstri_decomp(Anew, Bnew, U, V, work);
  if (flag == 1) gsl_vector_free(work);
  vH = Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, Anew);
  vR = Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, Bnew);
  if (argc == 2 || argc == 3) {
    ary = rb_ary_new3(2, vH, vR);
  } else {
    ary = rb_ary_new3(4, vH, vR, vU, vV);
  }
  return ary;
}
static VALUE rb_gsl_linalg_hesstri_decomp_bang(int argc, VALUE *argv, VALUE module)
{
  gsl_matrix *A = NULL, *B = NULL;
  gsl_matrix *U = NULL, *V = NULL;
  gsl_vector *work = NULL;
  VALUE vH, vR, vU = Qnil, vV = Qnil, ary;
  int flag = 0;
  switch (argc) {
  case 2:
    flag = 1;
    break;
  case 3:
    CHECK_VECTOR(argv[2]);
    Data_Get_Struct(argv[2], gsl_vector, work);
    break;
  case 4:
    CHECK_MATRIX(argv[2]);
    CHECK_MATRIX(argv[3]);
    Data_Get_Struct(argv[2], gsl_matrix, U);
    Data_Get_Struct(argv[3], gsl_matrix, V);
    flag = 1;
    break;
  case 5:
    CHECK_MATRIX(argv[2]);
    CHECK_MATRIX(argv[3]);
    CHECK_VECTOR(argv[4]);
    Data_Get_Struct(argv[2], gsl_matrix, U);
    Data_Get_Struct(argv[3], gsl_matrix, V);
    Data_Get_Struct(argv[4], gsl_vector, work);
    vU = argv[2];
    vV = argv[3];
    break;
  default:
    rb_raise(rb_eArgError, "Wrong number of arguments (%d for 2-55)", argc);
  }
  CHECK_MATRIX(argv[0]);
  CHECK_MATRIX(argv[1]);
  Data_Get_Struct(argv[0], gsl_matrix, A);
  Data_Get_Struct(argv[1], gsl_matrix, B);
  if (flag == 1) work = gsl_vector_alloc(A->size1);
  gsl_linalg_hesstri_decomp(A, B, U, V, work);
  if (flag == 1) gsl_vector_free(work);
  vH = argv[0];
  vR = argv[1];
  if (argc == 2 || argc == 3) {
    ary = rb_ary_new3(2, vH, vR);
  } else {
    ary = rb_ary_new3(4, vH, vR, vU, vV);
  }
  return ary;
}

static VALUE rb_gsl_linalg_balance_matrix(int argc, VALUE *argv, VALUE module)
{
  gsl_matrix *A, *Anew;
  gsl_vector *D;
  VALUE vA, vD;
  switch (argc) {
  case 1:
    CHECK_MATRIX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    Anew = make_matrix_clone(A);
    D = gsl_vector_alloc(A->size1);
    vD = Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, D);
    break;
  case 2:
    CHECK_MATRIX(argv[0]);
    CHECK_VECTOR(argv[1]);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    Data_Get_Struct(argv[1], gsl_vector, D);
    Anew = make_matrix_clone(A);
    vD = argv[1];
    break;
  default:
    rb_raise(rb_eArgError, "Wrong number of arguments (%d for 1 or 2)", argc);
  }
  gsl_linalg_balance_matrix(Anew, D);
  vA = Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, Anew);
  return rb_ary_new3(2, vA, vD);

}
static VALUE rb_gsl_linalg_balance_matrix2(int argc, VALUE *argv, VALUE module)
{
  gsl_matrix *A;
  gsl_vector *D;
  switch (argc) {
  case 1:
    CHECK_MATRIX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    D = gsl_vector_alloc(A->size1);
    gsl_linalg_balance_matrix(A, D);
    return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, D);
    break;
  case 2:
    CHECK_MATRIX(argv[0]);
    CHECK_VECTOR(argv[1]);
    Data_Get_Struct(argv[0], gsl_matrix, A);
    Data_Get_Struct(argv[1], gsl_vector, D);
    return INT2FIX(gsl_linalg_balance_matrix(A, D));
    break;
  default:
    rb_raise(rb_eArgError, "Wrong number of arguments (%d for 1 or 2)", argc);
  }
  return Qtrue;
}

void Init_gsl_linalg_complex(VALUE module);
void Init_gsl_linalg(VALUE module)
{
  VALUE mgsl_linalg;
  VALUE mgsl_linalg_LU;
  VALUE mgsl_linalg_QR;
  VALUE mgsl_linalg_QRPT;
  VALUE mgsl_linalg_LQ;
  VALUE mgsl_linalg_PTLQ;
  VALUE mgsl_linalg_SV;
  VALUE mgsl_linalg_cholesky;
  VALUE mgsl_linalg_symmtd;
  VALUE mgsl_linalg_hermtd;
  VALUE mgsl_linalg_bidiag;
  VALUE mgsl_linalg_tridiag;
  VALUE mgsl_linalg_HH;
  VALUE mgsl_linalg_Householder;
  VALUE mhessen;

  mgsl_linalg = rb_define_module_under(module, "Linalg");
  mgsl_linalg_LU = rb_define_module_under(mgsl_linalg, "LU");
  cgsl_matrix_LU = rb_define_class_under(mgsl_linalg_LU, "LUMatrix", cgsl_matrix);
  mgsl_linalg_QR = rb_define_module_under(mgsl_linalg, "QR");
  mgsl_linalg_QRPT = rb_define_module_under(mgsl_linalg, "QRPT");
  cgsl_matrix_QR = rb_define_class_under(mgsl_linalg, "QRMatrix", cgsl_matrix);
  cgsl_matrix_QRPT = rb_define_class_under(mgsl_linalg, "QRPTMatrix", cgsl_matrix);
  cgsl_vector_tau = rb_define_class_under(mgsl_linalg, "TauVector", cgsl_vector);
  cgsl_matrix_Q = rb_define_class_under(mgsl_linalg, "QMatrix", cgsl_matrix);
  cgsl_matrix_R = rb_define_class_under(mgsl_linalg, "RMatrix", cgsl_matrix);

  mgsl_linalg_LQ = rb_define_module_under(mgsl_linalg, "LQ");
  mgsl_linalg_PTLQ = rb_define_module_under(mgsl_linalg, "PTLQ");
  cgsl_matrix_LQ = rb_define_class_under(mgsl_linalg, "LQMatrix", cgsl_matrix);
  cgsl_matrix_PTLQ = rb_define_class_under(mgsl_linalg, "PTLQMatrix", cgsl_matrix);
  cgsl_matrix_L = rb_define_class_under(mgsl_linalg, "LMatrix", cgsl_matrix);

  /*****/
  mgsl_linalg_SV = rb_define_module_under(mgsl_linalg, "SV");
  cgsl_matrix_SV = rb_define_class_under(mgsl_linalg_SV, "SVMatrix", cgsl_matrix);
  cgsl_matrix_U = rb_define_class_under(mgsl_linalg_SV, "UMatrix", cgsl_matrix);
  cgsl_matrix_V = rb_define_class_under(mgsl_linalg_SV, "VMatrix", cgsl_matrix);
  cgsl_vector_S = rb_define_class_under(mgsl_linalg_SV, "SingularValues", cgsl_vector);

  /*****/
  mgsl_linalg_cholesky = rb_define_module_under(mgsl_linalg, "Cholesky");
  cgsl_matrix_C = rb_define_class_under(mgsl_linalg_cholesky, "CholeskyMatrix", cgsl_matrix);
  mgsl_linalg_symmtd = rb_define_module_under(mgsl_linalg, "Symmtd");

  mgsl_linalg_hermtd = rb_define_module_under(mgsl_linalg, "Hermtd");
  mgsl_linalg_bidiag = rb_define_module_under(mgsl_linalg, "Bidiag");
  mgsl_linalg_tridiag = rb_define_module_under(mgsl_linalg, "Tridiag");

  mgsl_linalg_HH = rb_define_module_under(mgsl_linalg, "HH");
  mgsl_linalg_Householder = rb_define_module_under(mgsl_linalg, "Householder");

  /*****/
  rb_define_module_function(mgsl_linalg, "LU_decomp!", rb_gsl_linalg_LU_decomp_bang, -1);
  rb_define_module_function(mgsl_linalg_LU, "decomp!", rb_gsl_linalg_LU_decomp_bang, -1);
  rb_define_module_function(mgsl_linalg, "LU_decomp", rb_gsl_linalg_LU_decomp, -1);
  rb_define_module_function(mgsl_linalg_LU, "decomp", rb_gsl_linalg_LU_decomp, -1);
  rb_define_method(cgsl_matrix, "LU_decomp!", rb_gsl_linalg_LU_decomp_bang, -1);
  rb_define_method(cgsl_matrix, "LU_decomp", rb_gsl_linalg_LU_decomp, -1);

  rb_define_module_function(mgsl_linalg, "LU_solve", rb_gsl_linalg_LU_solve, -1);
  rb_define_module_function(mgsl_linalg_LU, "solve", rb_gsl_linalg_LU_solve, -1);
  rb_define_method(cgsl_matrix, "LU_solve", rb_gsl_linalg_LU_solve, -1);
  rb_define_method(cgsl_matrix_LU, "solve", rb_gsl_linalg_LU_solve, -1);

  rb_define_module_function(mgsl_linalg, "LU_svx", rb_gsl_linalg_LU_svx, -1);
  rb_define_module_function(mgsl_linalg_LU, "svx", rb_gsl_linalg_LU_svx, -1);
  rb_define_method(cgsl_matrix, "LU_svx", rb_gsl_linalg_LU_svx, -1);
  rb_define_method(cgsl_matrix_LU, "svx", rb_gsl_linalg_LU_svx, -1);

  rb_define_module_function(mgsl_linalg, "LU_invert", rb_gsl_linalg_LU_invert, -1);
  rb_define_module_function(mgsl_linalg_LU, "invert", rb_gsl_linalg_LU_invert, -1);
  rb_define_module_function(mgsl_linalg_LU, "inv", rb_gsl_linalg_LU_invert, -1);
  rb_define_module_function(mgsl_linalg_LU, "refine", rb_gsl_linalg_LU_refine, 5);

  rb_define_method(cgsl_matrix, "invert", rb_gsl_linalg_LU_invert, -1);
  rb_define_alias(cgsl_matrix, "LU_invert", "invert");
  rb_define_alias(cgsl_matrix, "inv", "invert");

  rb_define_module_function(mgsl_linalg, "LU_det", rb_gsl_linalg_LU_det, -1);
  rb_define_module_function(mgsl_linalg_LU, "det", rb_gsl_linalg_LU_det, -1);
  rb_define_method(cgsl_matrix, "LU_det", rb_gsl_linalg_LU_det, -1);
  rb_define_alias(cgsl_matrix, "det", "LU_det");

  rb_define_module_function(mgsl_linalg, "LU_lndet", rb_gsl_linalg_LU_lndet, -1);
  rb_define_module_function(mgsl_linalg_LU, "lndet", rb_gsl_linalg_LU_lndet, -1);
  rb_define_method(cgsl_matrix, "LU_lndet", rb_gsl_linalg_LU_lndet, -1);
  rb_define_alias(cgsl_matrix, "lndet", "LU_lndet");

  rb_define_module_function(mgsl_linalg, "LU_sgndet", rb_gsl_linalg_LU_sgndet, -1);
  rb_define_module_function(mgsl_linalg_LU, "sgndet", rb_gsl_linalg_LU_sgndet, -1);
  rb_define_method(cgsl_matrix, "LU_sgndet", rb_gsl_linalg_LU_sgndet, -1);
  rb_define_alias(cgsl_matrix, "sgndet", "LU_sgndet");

  /*****/
  rb_define_module_function(mgsl_linalg, "QR_decomp", rb_gsl_linalg_QR_decomp, -1);
  rb_define_module_function(mgsl_linalg_QR, "decomp", rb_gsl_linalg_QR_decomp, -1);
  rb_define_method(cgsl_matrix, "QR_decomp", rb_gsl_linalg_QR_decomp, -1);
  rb_define_module_function(mgsl_linalg, "QR_decomp!", rb_gsl_linalg_QR_decomp_bang, -1);
  rb_define_module_function(mgsl_linalg_QR, "decomp!", rb_gsl_linalg_QR_decomp_bang, -1);
  rb_define_method(cgsl_matrix, "QR_decomp!", rb_gsl_linalg_QR_decomp_bang, -1);

  rb_define_module_function(mgsl_linalg, "QR_solve", rb_gsl_linalg_QR_solve, -1);
  rb_define_module_function(mgsl_linalg_QR, "solve", rb_gsl_linalg_QR_solve, -1);
  rb_define_module_function(mgsl_linalg, "QR_svx", rb_gsl_linalg_QR_svx, -1);
  rb_define_module_function(mgsl_linalg_QR, "svx", rb_gsl_linalg_QR_svx, -1);
  rb_define_method(cgsl_matrix, "QR_solve", rb_gsl_linalg_QR_solve, -1);
  rb_define_method(cgsl_matrix_QR, "solve", rb_gsl_linalg_QR_solve, -1);
  rb_define_method(cgsl_matrix, "QR_svx", rb_gsl_linalg_QR_svx, -1);
  rb_define_method(cgsl_matrix_QR, "svx", rb_gsl_linalg_QR_svx, -1);

  rb_define_module_function(mgsl_linalg_QR, "lssolve", rb_gsl_linalg_QR_lssolve, -1);
  rb_define_method(cgsl_matrix, "QR_lssolve", rb_gsl_linalg_QR_lssolve, -1);
  rb_define_method(cgsl_matrix_QR, "lssolve", rb_gsl_linalg_QR_lssolve, -1);

  rb_define_module_function(mgsl_linalg_QR, "QTvec", rb_gsl_linalg_QR_QTvec, -1);
  rb_define_method(cgsl_matrix_QR, "QTvec", rb_gsl_linalg_QR_QTvec, -1);
  rb_define_module_function(mgsl_linalg_QR, "Qvec", rb_gsl_linalg_QR_Qvec, -1);
  rb_define_method(cgsl_matrix_QR, "Qvec", rb_gsl_linalg_QR_Qvec, -1);

  rb_define_module_function(mgsl_linalg_QR, "Rsolve", rb_gsl_linalg_QR_Rsolve, -1);
  rb_define_method(cgsl_matrix, "QR_Rsolve", rb_gsl_linalg_QR_Rsolve, -1);
  rb_define_method(cgsl_matrix_QR, "Rsolve", rb_gsl_linalg_QR_Rsolve, -1);

  rb_define_module_function(mgsl_linalg_QR, "Rsvx", rb_gsl_linalg_QR_Rsvx, -1);
  rb_define_method(cgsl_matrix_QR, "Rsvx", rb_gsl_linalg_QR_Rsvx, -1);

  rb_define_module_function(mgsl_linalg_QR, "unpack", rb_gsl_linalg_QR_unpack, -1);
  rb_define_method(cgsl_matrix_QR, "unpack", rb_gsl_linalg_QR_unpack, -1);

  rb_define_module_function(mgsl_linalg_QR, "QRsolve", rb_gsl_linalg_QR_QRsolve, -1);
  rb_define_module_function(mgsl_linalg_QR, "update", rb_gsl_linalg_QR_update, 4);

  rb_define_method(mgsl_linalg, "R_solve", rb_gsl_linalg_R_solve, -1);
  rb_define_method(cgsl_matrix_R, "solve", rb_gsl_linalg_R_solve, -1);
  /*
  rb_define_method(cgsl_matrix_R, "svx", rb_gsl_linalg_R_svx, -1);
  */
  rb_define_module_function(mgsl_linalg_QRPT, "decomp", rb_gsl_linalg_QRPT_decomp, -1);
  rb_define_method(cgsl_matrix, "QRPT_decomp", rb_gsl_linalg_QRPT_decomp, -1);
  rb_define_module_function(mgsl_linalg_QRPT, "decomp!", rb_gsl_linalg_QRPT_decomp_bang, -1);
  rb_define_method(cgsl_matrix, "QRPT_decomp!", rb_gsl_linalg_QRPT_decomp_bang, -1);

  rb_define_module_function(mgsl_linalg_QRPT, "decomp2", rb_gsl_linalg_QRPT_decomp2, -1);
  rb_define_method(cgsl_matrix, "QRPT_decomp2", rb_gsl_linalg_QRPT_decomp2, -1);

  rb_define_module_function(mgsl_linalg_QRPT, "solve", rb_gsl_linalg_QRPT_solve, -1);
  rb_define_method(cgsl_matrix, "QRPT_solve", rb_gsl_linalg_QRPT_solve, -1);
  rb_define_method(cgsl_matrix_QRPT, "solve", rb_gsl_linalg_QRPT_solve, -1);

  rb_define_module_function(mgsl_linalg_QRPT, "svx", rb_gsl_linalg_QRPT_svx, -1);
  rb_define_method(cgsl_matrix, "QRPT_svx", rb_gsl_linalg_QRPT_svx, -1);
  rb_define_method(cgsl_matrix_QRPT, "svx", rb_gsl_linalg_QRPT_svx, -1);

  rb_define_module_function(mgsl_linalg_QRPT, "QRsolve", rb_gsl_linalg_QRPT_QRsolve, 4);
  rb_define_module_function(mgsl_linalg_QRPT, "update", rb_gsl_linalg_QRPT_update, 5);

  rb_define_module_function(mgsl_linalg_QRPT, "Rsolve", rb_gsl_linalg_QRPT_Rsolve, -1);
  rb_define_method(cgsl_matrix_QRPT, "Rsolve", rb_gsl_linalg_QRPT_Rsolve, -1);
  rb_define_module_function(mgsl_linalg_QRPT, "Rsvx", rb_gsl_linalg_QRPT_Rsvx, -1);
  rb_define_method(cgsl_matrix_QRPT, "Rsvx", rb_gsl_linalg_QRPT_Rsvx, -1);

  /*****/
  rb_define_module_function(mgsl_linalg_SV, "decomp", rb_gsl_linalg_SV_decomp, -1);
  rb_define_method(cgsl_matrix, "SV_decomp", rb_gsl_linalg_SV_decomp, -1);
  rb_define_alias(cgsl_matrix, "SVD", "SV_decomp");
  rb_define_alias(cgsl_matrix, "svd", "SV_decomp");
  rb_define_module_function(mgsl_linalg_SV, "decomp_mod", rb_gsl_linalg_SV_decomp_mod, -1);
  rb_define_method(cgsl_matrix, "SV_decomp_mod", rb_gsl_linalg_SV_decomp_mod, -1);
  rb_define_module_function(mgsl_linalg_SV, "decomp_jacobi", rb_gsl_linalg_SV_decomp_jacobi, -1);
  rb_define_method(cgsl_matrix, "SV_decomp_jacobi", rb_gsl_linalg_SV_decomp_jacobi, -1);

  rb_define_module_function(mgsl_linalg_SV, "solve", rb_gsl_linalg_SV_solve, -1);

  rb_define_method(cgsl_matrix, "SV_solve", rb_gsl_linalg_SV_solve, -1);

  /*****/
  rb_define_module_function(mgsl_linalg_cholesky, "decomp", rb_gsl_linalg_cholesky_decomp, -1);
  rb_define_method(cgsl_matrix, "cholesky_decomp", rb_gsl_linalg_cholesky_decomp, -1);

  rb_define_module_function(mgsl_linalg_cholesky, "solve", rb_gsl_linalg_cholesky_solve, -1);
  rb_define_method(cgsl_matrix, "cholesky_solve", rb_gsl_linalg_cholesky_solve, -1);
  rb_define_method(cgsl_matrix_C, "solve", rb_gsl_linalg_cholesky_solve, -1);
  rb_define_module_function(mgsl_linalg_cholesky, "svx", rb_gsl_linalg_cholesky_svx, -1);
  rb_define_method(cgsl_matrix, "cholesky_svx", rb_gsl_linalg_cholesky_svx, -1);
  rb_define_method(cgsl_matrix_C, "svx", rb_gsl_linalg_cholesky_svx, -1);

  /*****/

  rb_define_module_function(mgsl_linalg_symmtd, "decomp", rb_gsl_linalg_symmtd_decomp, -1);
  rb_define_method(cgsl_matrix, "symmtd_decomp", rb_gsl_linalg_symmtd_decomp, -1);
  rb_define_module_function(mgsl_linalg_symmtd, "decomp!", rb_gsl_linalg_symmtd_decomp2, -1);
  rb_define_method(cgsl_matrix, "symmtd_decomp!", rb_gsl_linalg_symmtd_decomp2, -1);

  rb_define_method(cgsl_matrix, "symmtd_unpack", rb_gsl_linalg_symmtd_unpack, -1);
  rb_define_method(cgsl_matrix, "symmtd_unpack_T", rb_gsl_linalg_symmtd_unpack_T, -1);

  rb_define_module_function(mgsl_linalg_symmtd, "unpack", rb_gsl_linalg_symmtd_unpack, -1);
  rb_define_module_function(mgsl_linalg_symmtd, "unpack_T", rb_gsl_linalg_symmtd_unpack_T, -1);
  /*****/
  rb_define_module_function(mgsl_linalg_hermtd, "decomp", rb_gsl_linalg_hermtd_decomp, -1);
  rb_define_method(cgsl_matrix, "hermtd_decomp", rb_gsl_linalg_hermtd_decomp, -1);
  rb_define_module_function(mgsl_linalg_hermtd, "decomp!", rb_gsl_linalg_hermtd_decomp2, -1);
  rb_define_method(cgsl_matrix, "hermtd_decomp!", rb_gsl_linalg_hermtd_decomp2, -1);

  rb_define_method(cgsl_matrix_complex, "hermtd_unpack", rb_gsl_linalg_hermtd_unpack, -1);
  rb_define_module_function(mgsl_linalg_hermtd, "unpack", rb_gsl_linalg_hermtd_unpack, -1);
  rb_define_method(cgsl_matrix_complex, "hermtd_unpack_T", rb_gsl_linalg_hermtd_unpack_T, -1);
  rb_define_module_function(mgsl_linalg_hermtd, "unpack_T", rb_gsl_linalg_hermtd_unpack_T, -1);

  /*****/
  rb_define_method(cgsl_matrix, "bidiag_decomp", rb_gsl_linalg_bidiag_decomp, -1);
  rb_define_method(cgsl_matrix, "bidiag_decomp!", rb_gsl_linalg_bidiag_decomp2, -1);

  rb_define_module_function(mgsl_linalg, "bidiag_decomp", rb_gsl_linalg_bidiag_decomp, -1);
  rb_define_module_function(mgsl_linalg, "bidiag_decomp!", rb_gsl_linalg_bidiag_decomp2, -1);
  rb_define_module_function(mgsl_linalg_bidiag, "decomp", rb_gsl_linalg_bidiag_decomp, -1);
  rb_define_module_function(mgsl_linalg_bidiag, "decomp!", rb_gsl_linalg_bidiag_decomp2, -1);

  rb_define_method(cgsl_matrix, "bidiag_unpack", rb_gsl_linalg_bidiag_unpack, -1);
  rb_define_method(cgsl_matrix, "bidiag_unpack2", rb_gsl_linalg_bidiag_unpack2, -1);
  rb_define_module_function(mgsl_linalg, "bidiag_unpack", rb_gsl_linalg_bidiag_unpack, -1);
  rb_define_module_function(mgsl_linalg, "bidiag_unpack2", rb_gsl_linalg_bidiag_unpack2, -1);
  rb_define_module_function(mgsl_linalg_bidiag, "unpack", rb_gsl_linalg_bidiag_unpack, -1);
  rb_define_module_function(mgsl_linalg_bidiag, "unpack2", rb_gsl_linalg_bidiag_unpack2, -1);

  rb_define_method(cgsl_matrix, "bidiag_unpack_B", rb_gsl_linalg_bidiag_unpack_B, -1);
  rb_define_module_function(mgsl_linalg, "bidiag_unpack_B", rb_gsl_linalg_bidiag_unpack_B, -1);
  rb_define_module_function(mgsl_linalg_bidiag, "unpack_B", rb_gsl_linalg_bidiag_unpack_B, -1);
  /*****/
  rb_define_module_function(mgsl_linalg, "householder_transform",
                            rb_gsl_linalg_householder_transform, -1);
  rb_define_module_function(mgsl_linalg_Householder, "transform",
                            rb_gsl_linalg_householder_transform, -1);
  rb_define_module_function(mgsl_linalg_HH, "transform",
                            rb_gsl_linalg_householder_transform, -1);
  rb_define_method(cgsl_vector, "householder_transform",
                   rb_gsl_linalg_householder_transform, -1);

  rb_define_module_function(mgsl_linalg, "householder_hm",
                            rb_gsl_linalg_householder_hm, 3);
  rb_define_module_function(mgsl_linalg_Householder, "hm",
                            rb_gsl_linalg_householder_hm, 3);
  rb_define_module_function(mgsl_linalg_HH, "hm",
                            rb_gsl_linalg_householder_hm, 3);

  rb_define_module_function(mgsl_linalg, "householder_mh",
                            rb_gsl_linalg_householder_mh, 3);
  rb_define_module_function(mgsl_linalg_Householder, "mh",
                            rb_gsl_linalg_householder_mh, 3);
  rb_define_module_function(mgsl_linalg_HH, "mh",
                            rb_gsl_linalg_householder_mh, 3);

  rb_define_module_function(mgsl_linalg, "householder_hv",
                            rb_gsl_linalg_householder_hv, 3);
  rb_define_module_function(mgsl_linalg_Householder, "hv",
                            rb_gsl_linalg_householder_hv, 3);
  rb_define_module_function(mgsl_linalg_HH, "hv",
                            rb_gsl_linalg_householder_hv, 3);

  rb_define_module_function(mgsl_linalg_HH, "solve", rb_gsl_linalg_HH_solve, -1);
  rb_define_module_function(mgsl_linalg_HH, "solve!", rb_gsl_linalg_HH_solve_bang, -1);
  rb_define_method(cgsl_matrix, "HH_solve", rb_gsl_linalg_HH_solve, -1);
  rb_define_method(cgsl_matrix, "HH_solve!", rb_gsl_linalg_HH_solve_bang, -1);
  rb_define_module_function(mgsl_linalg_HH, "svx", rb_gsl_linalg_HH_svx, -1);
  rb_define_method(cgsl_matrix, "HH_svx", rb_gsl_linalg_HH_svx, -1);

  /*****/

  rb_define_module_function(mgsl_linalg, "solve_symm_tridiag", rb_gsl_linalg_solve_symm_tridiag, 3);

  rb_define_module_function(mgsl_linalg_tridiag, "solve_symm", rb_gsl_linalg_solve_symm_tridiag, 3);

  rb_define_module_function(mgsl_linalg, "solve_tridiag", rb_gsl_linalg_solve_tridiag, 4);
  rb_define_module_function(mgsl_linalg_tridiag, "solve", rb_gsl_linalg_solve_tridiag, 4);
  rb_define_module_function(mgsl_linalg, "solve_symm_cyc_tridiag", rb_gsl_linalg_solve_symm_cyc_tridiag, 3);
  rb_define_module_function(mgsl_linalg, "solve_cyc_tridiag", rb_gsl_linalg_solve_cyc_tridiag, 4);
  rb_define_module_function(mgsl_linalg_tridiag, "solve_symm_cyc", rb_gsl_linalg_solve_symm_cyc_tridiag, 3);
  rb_define_module_function(mgsl_linalg_tridiag, "solve_cyc", rb_gsl_linalg_solve_cyc_tridiag, 4);

  /*****/
  rb_define_module_function(mgsl_linalg, "balance_columns!",
                            rb_gsl_linalg_balance_columns_bang, -1);
  rb_define_method(cgsl_matrix, "balance_columns!",
                   rb_gsl_linalg_balance_columns_bang, -1);
  rb_define_module_function(mgsl_linalg, "balance_columns",
                            rb_gsl_linalg_balance_columns, -1);
  rb_define_method(cgsl_matrix, "balance_columns",
                   rb_gsl_linalg_balance_columns, -1);
  rb_define_alias(cgsl_matrix, "balance", "balance_columns");
  rb_define_alias(cgsl_matrix, "balanc", "balance_columns");
  /*****/

  Init_gsl_linalg_complex(mgsl_linalg);

  /** GSL-1.6 **/
  rb_define_module_function(mgsl_linalg, "LQ_decomp", rb_gsl_linalg_LQ_decomp, -1);
  rb_define_module_function(mgsl_linalg_LQ, "decomp", rb_gsl_linalg_LQ_decomp, -1);
  rb_define_method(cgsl_matrix, "LQ_decomp", rb_gsl_linalg_LQ_decomp, -1);
  rb_define_module_function(mgsl_linalg, "LQ_decomp!", rb_gsl_linalg_LQ_decomp_bang, -1);
  rb_define_module_function(mgsl_linalg_LQ, "decomp!", rb_gsl_linalg_LQ_decomp_bang, -1);
  rb_define_method(cgsl_matrix, "LQ_decomp!", rb_gsl_linalg_LQ_decomp_bang, -1);

  rb_define_module_function(mgsl_linalg, "LQ_solve_T", rb_gsl_linalg_LQ_solve, -1);
  rb_define_module_function(mgsl_linalg_LQ, "solve_T", rb_gsl_linalg_LQ_solve, -1);
  rb_define_module_function(mgsl_linalg, "LQ_svx_T", rb_gsl_linalg_LQ_svx, -1);
  rb_define_module_function(mgsl_linalg_LQ, "svx_T", rb_gsl_linalg_LQ_svx, -1);
  rb_define_method(cgsl_matrix, "LQ_solve_T", rb_gsl_linalg_LQ_solve, -1);
  rb_define_method(cgsl_matrix_LQ, "solve_T", rb_gsl_linalg_LQ_solve, -1);
  rb_define_method(cgsl_matrix, "LQ_svx_T", rb_gsl_linalg_LQ_svx, -1);
  rb_define_method(cgsl_matrix_LQ, "svx_T", rb_gsl_linalg_LQ_svx, -1);

  rb_define_module_function(mgsl_linalg_LQ, "lssolve_T", rb_gsl_linalg_LQ_lssolve, -1);
  rb_define_method(cgsl_matrix, "LQ_lssolve_T", rb_gsl_linalg_LQ_lssolve, -1);
  rb_define_method(cgsl_matrix_LQ, "lssolve_T", rb_gsl_linalg_LQ_lssolve, -1);

  rb_define_module_function(mgsl_linalg_LQ, "vecQT", rb_gsl_linalg_LQ_vecQT, -1);
  rb_define_method(cgsl_matrix_LQ, "vecQT", rb_gsl_linalg_LQ_vecQT, -1);
  rb_define_module_function(mgsl_linalg_LQ, "vecQ", rb_gsl_linalg_LQ_vecQ, -1);
  rb_define_method(cgsl_matrix_LQ, "vecQ", rb_gsl_linalg_LQ_vecQ, -1);

  rb_define_module_function(mgsl_linalg_LQ, "Lsolve_T", rb_gsl_linalg_LQ_Lsolve, -1);
  rb_define_method(cgsl_matrix, "LQ_Lsolve_T", rb_gsl_linalg_LQ_Lsolve, -1);
  rb_define_method(cgsl_matrix_LQ, "Lsolve_T", rb_gsl_linalg_LQ_Lsolve, -1);

  rb_define_module_function(mgsl_linalg_LQ, "Lsvx_T", rb_gsl_linalg_LQ_Lsvx, -1);
  rb_define_method(cgsl_matrix_LQ, "Lsvx_T", rb_gsl_linalg_LQ_Lsvx, -1);

  rb_define_module_function(mgsl_linalg_LQ, "unpack", rb_gsl_linalg_LQ_unpack, -1);
  rb_define_method(cgsl_matrix_LQ, "unpack", rb_gsl_linalg_LQ_unpack, -1);

  rb_define_module_function(mgsl_linalg_LQ, "LQsolve_T", rb_gsl_linalg_LQ_LQsolve, -1);
  rb_define_module_function(mgsl_linalg_LQ, "update", rb_gsl_linalg_LQ_update, 4);

  rb_define_method(mgsl_linalg, "L_solve_T", rb_gsl_linalg_L_solve, -1);
  rb_define_method(cgsl_matrix_L, "solve_T", rb_gsl_linalg_L_solve, -1);
  /*
  rb_define_method(cgsl_matrix_R, "svx", rb_gsl_linalg_R_svx, -1);
  */
  rb_define_module_function(mgsl_linalg_PTLQ, "decomp", rb_gsl_linalg_PTLQ_decomp, -1);
  rb_define_method(cgsl_matrix, "PTLQ_decomp", rb_gsl_linalg_PTLQ_decomp, -1);
  rb_define_module_function(mgsl_linalg_PTLQ, "decomp!", rb_gsl_linalg_PTLQ_decomp_bang, -1);
  rb_define_method(cgsl_matrix, "PTLQ_decomp!", rb_gsl_linalg_PTLQ_decomp_bang, -1);

  rb_define_module_function(mgsl_linalg_PTLQ, "decomp2", rb_gsl_linalg_PTLQ_decomp2, -1);
  rb_define_method(cgsl_matrix, "PTLQ_decomp2", rb_gsl_linalg_PTLQ_decomp2, -1);

  rb_define_module_function(mgsl_linalg_PTLQ, "solve_T", rb_gsl_linalg_PTLQ_solve, -1);
  rb_define_method(cgsl_matrix, "PTLQ_solve_T", rb_gsl_linalg_PTLQ_solve, -1);
  rb_define_method(cgsl_matrix_PTLQ, "solve_T", rb_gsl_linalg_PTLQ_solve, -1);

  rb_define_module_function(mgsl_linalg_PTLQ, "svx_T", rb_gsl_linalg_PTLQ_svx, -1);
  rb_define_method(cgsl_matrix, "PTLQ_svx_T", rb_gsl_linalg_PTLQ_svx, -1);
  rb_define_method(cgsl_matrix_PTLQ, "svx_T", rb_gsl_linalg_PTLQ_svx, -1);

  rb_define_module_function(mgsl_linalg_PTLQ, "LQsolve_T", rb_gsl_linalg_PTLQ_LQsolve, 4);
  rb_define_module_function(mgsl_linalg_PTLQ, "update", rb_gsl_linalg_PTLQ_update, 5);

  rb_define_module_function(mgsl_linalg_PTLQ, "Lsolve_T", rb_gsl_linalg_PTLQ_Lsolve, -1);
  rb_define_method(cgsl_matrix_PTLQ, "Lsolve_T", rb_gsl_linalg_PTLQ_Lsolve, -1);
  rb_define_module_function(mgsl_linalg_PTLQ, "Lsvx_T", rb_gsl_linalg_PTLQ_Lsvx, -1);
  rb_define_method(cgsl_matrix_PTLQ, "Lsvx_T", rb_gsl_linalg_PTLQ_Lsvx, -1);

  mhessen = rb_define_module_under(mgsl_linalg, "Hessenberg");
  rb_define_module_function(mhessen, "decomp", rb_gsl_linalg_hessenberg_decomp, 1);
  rb_define_module_function(mgsl_linalg, "heesenberg_decomp", rb_gsl_linalg_hessenberg_decomp, 1);
  rb_define_module_function(mhessen, "unpack", rb_gsl_linalg_hessenberg_unpack, 2);
  rb_define_module_function(mgsl_linalg, "hessenberg_unpack", rb_gsl_linalg_hessenberg_unpack, 2);
  rb_define_module_function(mhessen, "unpack_accum", rb_gsl_linalg_hessenberg_unpack_accum, -1);
  rb_define_module_function(mgsl_linalg, "hessenberg_unpack_accum", rb_gsl_linalg_hessenberg_unpack_accum, -1);
  rb_define_module_function(mhessen, "set_zero", rb_gsl_linalg_hessenberg_set_zero, 1);
  rb_define_module_function(mgsl_linalg, "hessenberg_set_zero", rb_gsl_linalg_hessenberg_set_zero, 1);

  rb_define_module_function(mgsl_linalg, "hesstri_decomp", rb_gsl_linalg_hesstri_decomp, -1);
  rb_define_module_function(mgsl_linalg, "hesstri_decomp!", rb_gsl_linalg_hesstri_decomp_bang, -1);

  rb_define_module_function(mgsl_linalg, "balance_matrix", rb_gsl_linalg_balance_matrix, -1);
  rb_define_module_function(mgsl_linalg, "balance_matrix!", rb_gsl_linalg_balance_matrix2, -1);
  rb_define_module_function(mgsl_linalg, "balance", rb_gsl_linalg_balance_matrix, -1);
  rb_define_module_function(mgsl_linalg, "balance!", rb_gsl_linalg_balance_matrix2, -1);
}