#include <ruby.h>
#include <math.h>
#include <stdio.h>
#include <narray.h>
#include "ssl2.h"

extern void eig1_(float*,int*,int*,int*,float*,float*,float*,float*,int*);
extern void ceig2_(scomplex*,int*,int*,int*,scomplex*,scomplex*,float*,int*,int*);
extern void seig1_(float*,int*,float*,float*,int*,int*,float*,int*);
extern void seig2_(float*,int*,int*,float*,float*,int*,float*,int*);
extern void heig2_(float*,int*,int*,int*,float*,float*,float*,float*,int*);
extern void bseg_(float*,int*,int*,int*,int*,float*,float*,float*,int*,float*,int*);
extern void bsegj_(float*,int*,int*,int*,float*,int*,float*,float*,int*,int*,float*,int*);
extern void teig1_(float*,float*,int*,float*,float*,int*,int*,int*);
extern void teig2_(float*,float*,int*,int*,float*,float*,int*,float*,int*);
extern void gseg2_(float*,float*,int*,int*,float*,float*,float*,float*,int*,float*,int*);
extern void gbseg_(float*,float*,int*,int*,int*,float*,float*,int*,float*,float*,int*,int*,float*,int*);

extern void hsqr_(float*,int*,int*,float*,float*,int*,int*);
extern void chsqr_(scomplex*,int*,int*,scomplex*,int*,int*);
extern void trql_(float*,float*,int*,float*,int*,int*);
extern void bsct1_(float*,float*,int*,int*,float*,float*,float*,int*);

extern void hvec_(float*,int*,int*,float*,float*,int*,int*,float*,int*,float*,int*);
extern void chvec_(scomplex*,int*,int*,scomplex*,int*,int*,scomplex*,scomplex*,int*);
extern void bsvec_(float*,int*,int*,int*,float*,float*,int*,float*,int*);

extern void blnc_(float*,int*,int*,float*,int*);
extern void cblnc_(scomplex*,int*,int*,float*,int*);
extern void hes1_(float*,int*,int*,float*,int*);
extern void ches2_(scomplex*,int*,int*,int*,int*);
extern void trid1_(float*,int*,float*,float*,int*);
extern void tridh_(float*,int*,int*,float*,float*,float*,int*);
extern void btrid_(float*,int*,int*,float*,float*,int*);
extern void hbk1_(float*,int*,int*,int*,int*,float*,float*,float*,int*);
extern void chbk2_(scomplex*,int*,int*,int*,int*,scomplex*,int*,float*,int*);
extern void trbk_(float*,int*,int*,int*,float*,int*);
extern void trbkh_(float*,float*,int*,int*,int*,float*,float*,int*);
extern void nrml_(float*,int*,int*,int*,int*,int*,int*);
extern void cnrml_(scomplex*,int*,int*,int*,int*,int*,int*);
extern void gschl_(float*,float*,int*,float*,int*);
extern void gsbk_(float*,int*,int*,int*,float*,int*);


static VALUE
rb_eig1(argc, argv, self)
     int argc; VALUE *argv, self;
{
  float *ca;
  int n;
  int cmode;
  float *cer, *cei;
  float *cev;
  float *vw;
  int icon;

  VALUE a, mode;
  VALUE e,ev;
  scomplex *cec;
  int shape[2];

  int i,nn;

  rb_scan_args(argc, argv, "11", &a, &mode);
  if ( mode==Qnil )
    cmode = 0;
  else {
    switch ( TYPE(mode) ) {
    case T_FIXNUM:
      cmode = NUM2INT(mode);
      break;
    case T_TRUE:
      cmode = 1;
      break;
    case T_FALSE:
      cmode = 0;
      break;
    default:
      rb_raise(rb_eTypeError, "mode is invalid type");
      break;
    }
  }

  ca = ssl2_getcary(a,&nn);
  n = sqrt(nn);
  if ( nn != n*n ) rb_raise(rb_eRuntimeError, "a.length != n*n");
  cer = ALLOC_N(float,n);
  cei = ALLOC_N(float,n);
  cev = ALLOC_N(float,n*n);
  vw = ALLOC_N(float,n);

  eig1_(ca,&n,&n,&cmode,cer,cei,cev,vw,&icon);
  free(ca);
  free(vw);

  cec = ALLOC_N(scomplex,2*n);
  for (i=0;i<n;i++){
    cec[i].r = cer[i];
    cec[i].i = cei[i];
  }
  free(cer);
  free(cei);
  shape[0] = n;
  e = ssl2_getaryc(cec,1,shape);
  free(cec);
  shape[1] = n;
  ev = ssl2_getary(cev,2,shape);
  free(cev);

  ssl2_error(icon);
  return rb_ary_new3(2,e,ev);
}

static VALUE
rb_ceig2(argc, argv, self)
     int argc; VALUE *argv, self;
{
  scomplex *cza;
  int n;
  int cmode;
  scomplex *cze,*czev;
  float *vw;
  int *ivw;
  int icon;

  VALUE za, mode;
  VALUE ze,zev;
  int shape[2];

  int nn;

  rb_scan_args(argc, argv, "11", &za, &mode);
  if ( mode==Qnil )
    cmode = 0;
  else {
    switch ( TYPE(mode) ) {
    case T_FIXNUM:
      cmode = NUM2INT(mode);
      break;
    case T_TRUE:
      cmode = 1;
      break;
    case T_FALSE:
      cmode = 0;
      break;
    default:
      rb_raise(rb_eTypeError, "mode is invalid type");
      break;
    }
  }

  cza = ssl2_getcaryc(za,&nn);
  n = sqrt(nn);
  if ( nn != n*n ) rb_raise(rb_eRuntimeError, "za.length != n*n");

  cze = ALLOC_N(scomplex,n);
  czev = ALLOC_N(scomplex,n*n);
  vw = ALLOC_N(float,n);
  ivw = ALLOC_N(int,n);
  ceig2_(cza,&n,&n,&cmode,cze,czev,vw,ivw,&icon);
  free(cza);
  free(vw);
  free(ivw);

  shape[0] = n;
  ze = ssl2_getaryc(cze,1,shape);
  free(cze);
  shape[1] = n;
  zev = ssl2_getaryc(czev,2,shape);
  free(czev);

  ssl2_error(icon);
  return rb_ary_new3(2,ze,zev);
}

static VALUE
rb_seig1(self, a)
     VALUE self, a;
{
  float *ca;
  int n;
  float *ce,*cev;
  int m;
  float *vw;
  int icon;

  VALUE e,ev;
  int shape[2];

  int nn;

  ca = ssl2_getcary(a,&nn);
  n = (sqrt(8*nn+1)-1)/2;
  if ( nn != n*(n+1)/2 ) rb_raise(rb_eRuntimeError, "a.length != n(n+1)/2");

  ce = ALLOC_N(float,n); 
  cev = ALLOC_N(float,n*n);
  vw = ALLOC_N(float,2*n);
  seig1_(ca,&n,ce,cev,&n,&m,vw,&icon);
  free(ca);
  free(vw);

  shape[0] = m;
  e = ssl2_getary(ce,1,shape);
  free(ce);
  shape[0] = n;
  shape[1] = m;
  ev = ssl2_getary(cev,2,shape);
  free(cev);

  ssl2_error(icon);
  return rb_ary_new3(3,e,ev,INT2NUM(m));
}

static VALUE
rb_seig2(self, a, m)
     VALUE self, a, m;
{
  float *ca;
  int n;
  int cm;
  float *ce,*cev;
  float *vw;
  int icon;

  VALUE e,ev;
  int shape[2];

  int nn;

  cm = NUM2INT(m);
  ca = ssl2_getcary(a,&nn);
  n = (sqrt(8*nn+1)-1)/2;
  if ( nn != n*(n+1)/2 ) rb_raise(rb_eRuntimeError, "a.length != n(n+1)/2");
  if ( abs(cm)<1 )
    rb_raise(rb_eRuntimeError, "abs(m) must be >= 1");
  else if ( abs(cm)>n )
    rb_raise(rb_eRuntimeError, "abs(m) must be <= n");

  ce = ALLOC_N(float,abs(cm)); 
  cev = ALLOC_N(float,n*abs(cm));
  vw = ALLOC_N(float,7*n);
  seig2_(ca,&n,&cm,ce,cev,&n,vw,&icon);
  free(ca);
  free(vw);

  shape[0] = abs(cm);
  e = ssl2_getary(ce,1,shape);
  free(ce);
  shape[0] = n;
  shape[1] = abs(cm);
  ev = ssl2_getary(cev,2,shape);
  free(cev);

  ssl2_error(icon);
  return rb_ary_new3(2,e,ev);
}

static VALUE
rb_heig2(self, a, m)
     VALUE self, a, m;
{
  float *ca;
  int n;
  int cm;
  float *ce;
  float *cevr, *cevi;
  float *vw;
  int icon;

  VALUE e;
  VALUE evc;
  scomplex *cevc;
  int shape[2];

  int nn,i;

  cm = NUM2INT(m);
  ca = ssl2_getcary(a,&nn);
  n = sqrt(nn);
  if ( nn != n*n ) rb_raise(rb_eRuntimeError, "a.length != n*n");
  if ( abs(cm)<1 ) rb_raise(rb_eRuntimeError, "abs(m) must be >= 1");
  if ( abs(cm)>n ) rb_raise(rb_eRuntimeError, "abs(m) must be <= n");

  ce = ALLOC_N(float,abs(cm)); 
  cevr = ALLOC_N(float,n*abs(cm));
  cevi = ALLOC_N(float,n*abs(cm));
  vw = ALLOC_N(float,9*n);
  heig2_(ca,&n,&n,&cm,ce,cevr,cevi,vw,&icon);
  free(ca);
  free(vw);

  shape[0] = abs(cm);
  e = ssl2_getary(ce,1,shape);
  free(ce);
  cevc = ALLOC_N(scomplex,n*abs(cm));
  for (i=0;i<n*abs(cm);i++) {
    cevc[i].r = cevr[i];
    cevc[i].i = cevi[i];
  }
  free(cevr);
  free(cevi);
  shape[0] = n;
  shape[1] = abs(cm);
  evc = ssl2_getaryc(cevc,2,shape);
  free(cevc);

  ssl2_error(icon);
  return rb_ary_new3(2,e,evc);
}

static VALUE
rb_bseg(self, a, nh, m, nv, epst)
     VALUE self, a, nh, m, nv, epst;
{
  float *ca;
  int n;
  int cnh;
  int cm;
  int cnv;
  float cepst;
  float *ce,*cev;
  float *vw;
  int icon;

  VALUE e,ev;
  int shape[2];

  int nn;

  cnh = NUM2INT(nh);
  if ( cnh<0 ) rb_raise(rb_eRuntimeError, "nh must be >= 0");
  ca = ssl2_getcary(a,&nn);
  n = nn/(cnh+1)+cnh/2;
  if (nn!=n*(cnh+1)-cnh*(cnh+1)/2) rb_raise(rb_eRuntimeError, "a.length != n(nh+1)-nh(nh+1)/2");
  if ( cnh>=n ) rb_raise(rb_eRuntimeError, "nh >= n");
  cm = NUM2INT(m);
  if (abs(cm)<1) rb_raise(rb_eRuntimeError, "abs(m) must be >= 1");
  if (abs(cm)>n) rb_raise(rb_eRuntimeError, "abs(m) must be <= n");
  cnv = NUM2INT(nv);
  if (cnv<0) rb_raise(rb_eRuntimeError, "nv must be >= 0");
  if (cnv>m) rb_raise(rb_eRuntimeError, "nv must be <= m");
  cepst = NUM2FLT(epst);

  ce = ALLOC_N(float,abs(cm)); 
  cev = ALLOC_N(float,n*cnv);
  vw = ALLOC_N(float,MAX(3*n+2*cm,2*n*(cnh+1)));
  bseg_(ca,&n,&cnh,&cm,&cnv,&cepst,ce,cev,&n,vw,&icon);
  free(ca);
  free(vw);

  shape[0] = abs(cm);
  e = ssl2_getary(ce,1,shape);
  free(ce);
  shape[0] = n;
  shape[1] = cnv;
  ev = ssl2_getary(cev,2,shape);
  free(cev);

  ssl2_error(icon);
  return rb_ary_new3(2,e,ev);
}

static VALUE
rb_bsegj(self, a, nh, m, epst, lm, ev)
     VALUE self, a, nh, m, epst, lm, ev;
{
  float *ca;
  int n;
  int cnh;
  int cm;
  float cepst;
  int clm;
  float *ce;
  float *cevi, *cev;
  int it;
  float *vw;
  int icon;

  VALUE e;
  int shape[2];

  int i,nn;

  cnh = NUM2INT(nh);
  if ( cnh<0 ) rb_raise(rb_eRuntimeError, "nh must be >= 0");
  ca = ssl2_getcary(a,&nn);
  n = (nn+cnh*(cnh+1)/2)/(cnh+1);
  if ( nn!=n*(cnh+1)-cnh*(cnh+1)/2 ) rb_raise(rb_eRuntimeError, "a.length != n(nh+1)-nh(nh+1)/2");
  if ( cnh>=n ) rb_raise(rb_eRuntimeError, "nh>=n");
  cm = NUM2INT(m);
  if ( abs(cm)<1 ) rb_raise(rb_eRuntimeError, "abs(m) must be >= 1");
  if ( abs(cm)>n ) rb_raise(rb_eRuntimeError, "abs(m) must be <= n");
  cepst = NUM2FLT(epst);
  clm = NUM2INT(lm);
  if ( clm<1 ) rb_raise(rb_eRuntimeError, "lm must be > 0");
  cevi = ssl2_getcary(ev,&nn);
  if ( nn!=n*cm ) rb_raise(rb_eRuntimeError, "ve.length != n*cm");

  cev = ALLOC_N(float,n*(abs(m)+2));
  for (i=0;i<n*abs(m);i++)
    cev[i] = cevi[i];
  free(cevi);
  ce = ALLOC_N(float,abs(cm));
  vw = ALLOC_N(float,MAX(n,2*abs(cm))+abs(cm)*(3*abs(cm)+1)/2);
  bsegj_(ca,&n,&cnh,&cm,&cepst,&clm,ce,cev,&n,&it,vw,&icon);
  free(ca);
  free(vw);

  shape[0] = abs(cm);
  e = ssl2_getary(ce,1,shape);
  free(ce);
  shape[0] = n;
  shape[1] = abs(cm)+2;
  ev = ssl2_getary(cev,2,shape);
  free(cev);

  ssl2_error(icon);
  return rb_ary_new3(3,e,ev,INT2NUM(it));
}

static VALUE
rb_teig1(self, d, sd)
     VALUE self, d, sd;
{
  float *cd, *csd;
  int n;
  float *ce, *cev;
  int m;
  int icon;

  VALUE e, ev;
  int shape[2];

  int nn;

  cd = ssl2_getcary(d,&nn);
  n = nn;
  if ( n<1 ) rb_raise(rb_eRuntimeError, "d.length must be >= !");
  csd = ssl2_getcary(sd,&nn);
  if ( nn!=n ) rb_raise(rb_eRuntimeError, "sd.length != d.length");

  ce = ALLOC_N(float,n); 
  cev = ALLOC_N(float,n*n); 
  teig1_(cd, csd, &n, ce, cev, &n, &m, &icon);
  free(cd);
  free(csd);

  shape[0] = n;
  e = ssl2_getary(ce,1,shape);
  free(ce);
  shape[0] = n;
  shape[1] = n;
  ev = ssl2_getary(cev,2,shape);
  free(cev);

  ssl2_error(icon);
  return rb_ary_new3(3,e,ev,INT2NUM(m));
}

static VALUE
rb_teig2(self, d, sd, m)
     VALUE self, d, sd, m;
{
  float *cd, *csd;
  int n;
  int cm;
  float *ce, *cev;
  float *vw;
  int icon;

  VALUE e, ev;
  int shape[2];

  int nn;

  cm = NUM2INT(m);
  if ( abs(cm)<1 ) rb_raise(rb_eRuntimeError, "abs(m) must be >= 1");
  cd = ssl2_getcary(d,&n);
  if ( n<1 ) rb_raise(rb_eRuntimeError, "d.length must be >= !");
  csd = ssl2_getcary(sd,&nn);
  if ( nn!=n ) rb_raise(rb_eRuntimeError, "sd.length != d.length");
  if ( abs(cm)>n ) rb_raise(rb_eRuntimeError, "abs(m)>n");

  ce = ALLOC_N(float,abs(cm)); 
  cev = ALLOC_N(float,n*abs(cm)); 
  vw = ALLOC_N(float,5*n);
  teig2_(cd, csd, &n, &cm, ce, cev, &n, vw, &icon);
  free(cd);
  free(csd);
  free(vw);

  shape[0] = n;
  e = ssl2_getary(ce,1,shape);
  free(ce);
  shape[1] = n;
  ev = ssl2_getary(cev,2,shape);
  free(cev);

  ssl2_error(icon);
  return rb_ary_new3(2,e,ev);
}

static VALUE
rb_gseg2(self, a, b, m, epsz, epst)
     VALUE self, a, b, m, epsz, epst;
{
  float *ca, *cb;
  int n;
  int cm;
  float cepsz, cepst;
  float *ce, *cev;
  float *vw;
  int icon;

  VALUE e, ev;
  int shape[2];

  int nn;

  ca = ssl2_getcary(a,&nn);
  n = (sqrt(8*nn+1)-1)/2;
  if ( nn!=n*(n+1)/2 ) rb_raise(rb_eRuntimeError, "a.length must be = n*(n+1)/2");
  cb = ssl2_getcary(b,&nn);
  if ( nn!=n*(n+1)/2 ) rb_raise(rb_eRuntimeError, "b.length must be = n*(n+1)/2");
  cm = NUM2INT(m);
  if ( abs(cm)>n ) rb_raise(rb_eRuntimeError, "abs(m) must be <= n");
  cepsz = NUM2FLT(epsz);
  cepst = NUM2FLT(epst);

  ce = ALLOC_N(float,abs(cm)); 
  cev = ALLOC_N(float,n*abs(cm)); 
  vw = ALLOC_N(float,7*n);
  gseg2_(ca, cb, &n, &cm, &cepsz, &cepst, ce, cev, &n, vw, &icon);
  free(ca);
  free(cb);
  free(vw);

  shape[0] = abs(cm);
  e = ssl2_getary(ce,1,shape);
  free(ce);
  shape[0] = n;
  shape[1] = abs(cm);
  ev = ssl2_getary(cev,2,shape);
  free(cev);

  ssl2_error(icon);
  return rb_ary_new3(2,e,ev);
}

static VALUE
rb_gbseg(self, a, b, nh, m, epsz, epst, lm, ev)
     VALUE self, a, b, nh, m, epsz, epst, lm, ev;
{
  float *ca, *cb;
  int n;
  int cnh;
  int cm;
  float cepsz, cepst;
  int clm;
  float *ce;
  float *cevi, *cev;
  int it;
  float *vw;
  int icon;

  VALUE e;
  int shape[2];

  int i,nn;

  cnh = NUM2INT(nh);
  if ( cnh<0 ) rb_raise(rb_eRuntimeError, "nh must be >= 0");
  ca = ssl2_getcary(a,&nn);
  n = nn/(cnh+1)+cnh/2;
  if (nn!=n*(cnh+1)-cnh*(cnh+1)/2) rb_raise(rb_eRuntimeError, "a.length must be = n(nh+1)-nh(nh+1)/2");
  if ( cnh>=n ) rb_raise(rb_eRuntimeError, "nh>=n");
  cb = ssl2_getcary(b,&nn);
  if (nn!=n*(cnh+1)-cnh*(cnh+1)/2) rb_raise(rb_eRuntimeError, "b.length must be = n(nh+1)-nh(nh+1)/2");
  cm = NUM2INT(m);
  if ( abs(cm)>n ) rb_raise(rb_eRuntimeError, "abs(m) must be <= n");
  clm = NUM2INT(lm);
  if ( clm==0 ) rb_raise(rb_eRuntimeError, "lm must be >= 0");
  cepsz = NUM2FLT(epsz);
  cepst = NUM2FLT(epst);
  cevi = ssl2_getcary(ev,&nn);
  if ( nn!=n*cm ) rb_raise(rb_eRuntimeError, "ev.length must be = n*m");

  cev = ALLOC_N(float,n*(abs(cm)+2));
  for (i=0;i<n*abs(cm);i++)
    cev[i] = cevi[i];
  free(cevi);
  ce = ALLOC_N(float,abs(cm)); 
  vw = ALLOC_N(float,2*n+abs(cm)*(3*abs(cm)+1)/2);
  gbseg_(ca, cb, &n, &cnh, &cm, &cepsz, &cepst, &clm, ce, cev, &n, &it, vw, &icon);
  free(ca);
  free(cb);
  free(vw);

  shape[0] = abs(cm);
  e = ssl2_getary(ce,1,shape);
  free(ce);
  shape[0] = abs(n);
  shape[1] = abs(cm)+2;
  ev = ssl2_getary(cev,2,shape);
  free(cev);

  ssl2_error(icon);
  return rb_ary_new3(3,e,ev,INT2NUM(it));
}

static VALUE
rb_hsqr(self,a)
     VALUE self,a;
{
  float *ca;
  int n;
  float *cer,*cei;
  int m;
  int icon;

  VALUE ec;
  scomplex *cec;
  int *shape,rank;

  int i;

  ca = ssl2_getcary_withshape(a,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"a.rank!=2");
  n = shape[0];
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"shape[0]!=shape[1]");
  if (n==1) rb_raise(rb_eRuntimeError,"n==1");

  cer = ALLOC_N(float,n);
  cei = ALLOC_N(float,n);
  hsqr_(ca,&n,&n,cer,cei,&m,&icon);
  free(ca);

  cec = ALLOC_N(scomplex,m);
  for (i=0;i<m;i++) {
    cec[i].r = cer[i];
    cec[i].i = cei[i];
  }
  free(cer);
  free(cei);

  shape[0] = m;
  ec = ssl2_getaryc(cec,1,shape);
  free(cec);
  free(shape);

  ssl2_error(icon);
  return ec;
}

static VALUE
rb_chsqr(self,za)
     VALUE self,za;
{
  scomplex *cza;
  int n;
  scomplex *cze;
  int m;
  int icon;

  VALUE ze;
  int *shape,rank;

  cza = ssl2_getcaryc_withshape(za,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"za.rank!=2");
  n = shape[0];
  if (shape[1]) rb_raise(rb_eRuntimeError,"za.shape[0]!=za.shape[1]");
  if (n<1) rb_raise(rb_eRuntimeError,"za.shape[0]<1");

  cze = ALLOC_N(scomplex,n);
  chsqr_(cza,&n,&n,cze,&m,&icon);
  free(cza);

  shape[0] = m;
  ze = ssl2_getaryc(cze,1,shape);
  free(cze);
  free(shape);

  ssl2_error(icon);
  return ze;
}

static VALUE
rb_trql(self,d,sd)
     VALUE self,d,sd;
{
  float *cd,*csd;
  int n;
  float *ce;
  int m;
  int icon;

  VALUE e;
  int shape[1];

  int nn;

  cd = ssl2_getcary(d,&n);
  if (n<=1) rb_raise(rb_eRuntimeError,"n<=1");
  csd = ssl2_getcary(sd,&nn);
  if (nn!=n) rb_raise(rb_eRuntimeError,"sd.length != d.length");

  ce = ALLOC_N(float,n);
  trql_(cd,csd,&n,ce,&m,&icon);
  free(cd);
  free(csd);

  shape[0] = m;
  e = ssl2_getary(ce,1,shape);
  free(ce);

  ssl2_error(icon);
  return e;
}

static VALUE
rb_bsct1(self,d,sd,m,epst)
     VALUE self,d,sd,m,epst;
{
  float *cd,*csd;
  int n;
  int cm;
  float cepst;
  float *ce;
  float *vw;
  int icon;

  VALUE e;
  int shape[1];

  int nn;

  cd = ssl2_getcary(d,&n);
  if (n==1) rb_raise(rb_eRuntimeError,"n==1");
  csd = ssl2_getcary(sd,&nn);
  if (nn!=n) rb_raise(rb_eRuntimeError,"sd.length != d.length");
  cm = NUM2INT(m);
  if (abs(cm)>n) rb_raise(rb_eRuntimeError,"abs(m)>n");
  if (cm==0) rb_raise(rb_eRuntimeError,"m==0");
  cepst = NUM2FLT(epst);

  ce = ALLOC_N(float,abs(cm));
  vw = ALLOC_N(float,n+2*cm);
  bsct1_(cd,csd,&n,&cm,&cepst,ce,vw,&icon);
  free(cd);
  free(csd);
  free(vw);

  shape[0] = m;
  e = ssl2_getary(ce,1,shape);
  free(ce);

  ssl2_error(icon);
  return e;
}

static VALUE
rb_hvec(self,a,e,ind)
     VALUE self,a,e,ind;
{
  float *ca;
  int n;
  scomplex *ce;
  float *cer,*cei;
  int *cind;
  int m;
  float *cev;
  float *aw;
  int icon;

  VALUE ev;
  int *shape,rank;

  int i,nn;

  ca = ssl2_getcary_withshape(a,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"a.rank != 2");
  n = shape[0];
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"shape[0] != shape[1]");
  if (n==1) rb_raise(rb_eRuntimeError,"n==1");
  ce = ssl2_getcaryc(e,&m);
  if (n<m) rb_raise(rb_eRuntimeError,"n<m");
  cer = ALLOC_N(float,m);
  cei = ALLOC_N(float,m);
  for (i=0;i<m;i++) {
    cer[i] = ce[i].r;
    cei[i] = ce[i].i;
  }
  free(ce);
  cind = ssl2_getcaryi(ind,&nn);
  if (nn!=m) rb_raise(rb_eRuntimeError,"ind.length != e.length");

  cev = ALLOC_N(float,n*m);
  aw = ALLOC_N(float,n*(n+4));
  hvec_(ca,&n,&n,cer,cei,cind,&m,cev,&m,aw,&icon);
  free(ca);
  free(cer);
  free(cei);
  free(cind);
  free(aw);

  shape[0] = m;
  ev = ssl2_getary(cev,1,shape);
  free(cev);
  free(shape);

  ssl2_error(icon);
  return ev;
}

static VALUE
rb_chvec(self,za,ze,ind)
     VALUE self,za,ze,ind;
{
  scomplex *cza;
  int n;
  scomplex *cze;
  int *cind;
  int m;
  scomplex *czev;
  scomplex *zaw;
  int icon;

  VALUE zev;
  int *shape,rank;

  int nn;

  cza = ssl2_getcaryc_withshape(za,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"za.rank!=2");
  n = shape[0];
  if (n<1) rb_raise(rb_eRuntimeError,"za.shape[0]<1");
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"za.shape[0]!=za.shape[1]");
  cze = ssl2_getcaryc(ze,&m);
  if (m<1) rb_raise(rb_eRuntimeError,"ze.length<1");
  if (m>n) rb_raise(rb_eRuntimeError,"ze.length>n");
  cind = ssl2_getcaryi(ind,&nn);
  if (nn!=m) rb_raise(rb_eRuntimeError,"ind.length != ze.length");

  czev = ALLOC_N(scomplex,n*m);
  zaw = ALLOC_N(scomplex,n*(n+1));
  chvec_(cza,&n,&n,cze,cind,&m,czev,zaw,&icon);
  free(cza);
  free(cze);
  free(cind);

  shape[1] = m;
  zev = ssl2_getaryc(cze,2,shape);
  free(cze);
  free(shape);

  ssl2_error(icon);
  return zev;
}

static VALUE
rb_bsvec(self,a,nh,e)
     VALUE self,a,nh,e;
{
  float *ca;
  int n;
  int cnh;
  int nv;
  float *ce;
  float *cev;
  float *vw;
  int icon;

  VALUE ev;
  int shape[2];

  int nn;

  cnh = NUM2INT(nh);
  if (cnh<0) rb_raise(rb_eRuntimeError,"nh<0");
  ca = ssl2_getcary(a,&nn);
  n = (nn+cnh*(cnh+1)/2)/(cnh+1);
  if (n*(cnh+1)-cnh*(cnh+1)/2!=nn) rb_raise(rb_eRuntimeError,"a.length != n(nh+1)-nh(nh+1)/2");
  ce = ssl2_getcary(e,&nv);
  if (nv>n) rb_raise(rb_eRuntimeError,"e.length > n");

  cev = ALLOC_N(float,n*nv);
  vw = ALLOC_N(float,2*n*(cnh+1));
  bsvec_(ca,&n,&cnh,&nv,ce,cev,&n,vw,&icon);
  free(ca);
  free(ce);
  free(vw);

  shape[0] = n;
  shape[1] = nv;
  ev = ssl2_getary(cev,2,shape);
  free(cev);

  ssl2_error(icon);
  return ev;
}

static VALUE
rb_blnc(self,a)
     VALUE self,a;
{
  float *ca;
  int n;
  float *cdv;
  int icon;

  VALUE dv;
  int *shape,rank;

  ca = ssl2_getcary_withshape(a,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"a.rank != 2");
  n = shape[0];
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"shape[0] != shape[1]");

  cdv = ALLOC_N(float,n);
  blnc_(ca,&n,&n,cdv,&icon);

  a = ssl2_getary(ca,2,shape);
  free(ca);
  shape[0] = n;
  dv = ssl2_getary(cdv,1,shape);
  free(cdv);
  free(shape);

  ssl2_error(icon);
  return rb_ary_new3(2,a,dv);
}

static VALUE
rb_cblnc(self,za)
     VALUE self,za;
{
  scomplex *cza;
  int n;
  float *cdv;
  int icon;

  VALUE dv;
  int *shape,rank;

  cza = ssl2_getcaryc_withshape(za,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"za.rank != 2");
  n = shape[0];
  if (n<=1) rb_raise(rb_eRuntimeError,"za.shape[0] <= 1");
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"za.shape[0]!=za.shape[1]");

  cdv = ALLOC_N(float,n);
  cblnc_(cza,&n,&n,cdv,&icon);

  za = ssl2_getaryc(cza,2,shape);
  free(cza);
  dv = ssl2_getary(cdv,1,shape);
  free(cdv);
  free(shape);

  ssl2_error(icon);
  return rb_ary_new3(2,za,dv);
}

static VALUE
rb_hes1(self,a)
     VALUE self,a;
{
  float *ca;
  int n;
  float *cpv;
  int icon;

  VALUE pv;
  int *shape,rank;

  ca = ssl2_getcary_withshape(a,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"a.rank != 2");
  n = shape[0];
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"shape[0] != shape[1]");
  if (n<=2) rb_raise(rb_eRuntimeError,"n<=2");

  cpv = ALLOC_N(float,n);
  hes1_(ca,&n,&n,cpv,&icon);

  a = ssl2_getary(ca,2,shape);
  free(ca);
  shape[0] = n;
  pv = ssl2_getary(cpv,1,shape);
  free(cpv);
  free(shape);

  ssl2_error(icon);
  return rb_ary_new3(2,a,pv);
}

static VALUE
rb_ches2(self,za)
     VALUE self,za;
{
  scomplex *cza;
  int n;
  int *cip;
  int icon;

  VALUE ip;
  int *shape,rank;

  cza = ssl2_getcaryc_withshape(za,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"za.rank!=2");
  n = shape[0];
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"za.shape[0]!=za.shape[1]");

  cip = ALLOC_N(int,n);
  ches2_(cza,&n,&n,cip,&icon);

  za = ssl2_getaryc(cza,2,shape);
  free(cza);
  ip = ssl2_getaryi(cip,1,shape);
  free(cip);
  free(shape);

  ssl2_error(icon);
  return rb_ary_new3(za,ip);
}

static VALUE
rb_trid1(self,a)
     VALUE self,a;
{
  float *ca;
  int n;
  float *cd,*csd;
  int icon;

  VALUE d,sd;
  int shape[1];

  int nn;

  ca = ssl2_getcary(a,&nn);
  n = (sqrt(8*nn+1)-1)/2;
  if (nn!=n*(n+1)/2) rb_raise(rb_eRuntimeError,"a.length != n(n+1)/2");

  cd = ALLOC_N(float,n);
  csd = ALLOC_N(float,n);
  trid1_(ca,&n,cd,csd,&icon);

  shape[0] = n*(n+1)/2;
  a = ssl2_getary(ca,1,shape);
  free(ca);
  shape[0] = n;
  d = ssl2_getary(cd,1,shape);
  free(cd);
  sd =ssl2_getary(csd,1,shape);
  free(csd);

  ssl2_error(icon);
  return rb_ary_new3(3,a,d,sd);
}

static VALUE
rb_tridh(self,a)
     VALUE self,a;
{
  float *ca;
  int n;
  float *cd,*csd,*cpv;
  int icon;

  VALUE d,sd,pv;
  int *shape,rank;

  ca = ssl2_getcary_withshape(a,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"a.rank!=2");
  n = shape[0];
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"a.shape[0] != a.shape[1]");
  if (n==1) rb_raise(rb_eRuntimeError,"n==1");

  cd = ALLOC_N(float,n);
  csd = ALLOC_N(float,n);
  cpv = ALLOC_N(float,2*n);
  tridh_(ca,&n,&n,cd,csd,cpv,&icon);

  a = ssl2_getary(ca,2,shape);
  free(ca);
  shape[0] = n;
  d = ssl2_getary(cd,1,shape);
  free(cd);
  sd = ssl2_getary(csd,1,shape);
  free(csd);
  shape[0] = 2*n;
  pv = ssl2_getary(cpv,1,shape);
  free(cpv);
  free(shape);

  ssl2_error(icon);
  return rb_ary_new3(4,a,d,sd,pv);
}

static VALUE
rb_btrid(self,a,nh)
     VALUE self,a,nh;
{
  float *ca;
  int n;
  int cnh;
  float *cd,*csd;
  int icon;

  VALUE d,sd;
  int shape[1];

  int nn;

  cnh = NUM2INT(nh);
  if(cnh<=1) rb_raise(rb_eRuntimeError,"nh<=1");
  ca = ssl2_getcary(a,&nn);
  n = (nn+cnh*(cnh+1)/2)/(cnh+1);
  if (nn!=n*(cnh+1)-cnh*(cnh+1)/2) rb_raise(rb_eRuntimeError,"a.lenght != n(nh+1)-nh(nh+1)/2");
  if (cnh>=n) rb_raise(rb_eRuntimeError,"nh>=n");

  cd = ALLOC_N(float,n);
  csd = ALLOC_N(float,n);
  btrid_(ca,&n,&cnh,cd,csd,&icon);
  free(ca);

  shape[0] = n;
  d = ssl2_getary(cd,1,shape);
  free(cd);
  sd = ssl2_getary(csd,1,shape);
  free(csd);

  ssl2_error(icon);
  return rb_ary_new3(2,d,sd);
}

static VALUE
rb_hbk1(self,ev,ind,p,pv,dv)
     VALUE self,ev,ind,p,pv,dv;
{
  float *cev;
  int n;
  int *cind;
  int m;
  float *cp,*cpv,*cdv;
  int icon;

  int *shape,rank;

  int nn;

  cev = ssl2_getcary_withshape(ev,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"ev.rank != 2");
  n = shape[0];
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"ev.shape[0] != ev.shape[1]");
  if (n==1) rb_raise(rb_eRuntimeError,"n==1");
  free(shape);
  cind = ssl2_getcaryi(ind,&m);
  if (m<1) rb_raise(rb_eRuntimeError,"m<1");
  if (n<m) rb_raise(rb_eRuntimeError,"ind.length > n");
  cp = ssl2_getcary_withshape(p,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"p.rank != 2");
  if (shape[0]!=n) rb_raise(rb_eRuntimeError,"p.shape[0] != n");
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"p.shape[1] != n");
  cpv = ssl2_getcary(pv,&nn);
  if (nn!=n) rb_raise(rb_eRuntimeError,"pv.length != n");
  cdv = ssl2_getcary(dv,&nn);
  if (nn!=n) rb_raise(rb_eRuntimeError,"dv.length != n");

  hbk1_(cev,&n,&n,cind,&m,cp,cpv,cdv,&icon);
  free(cind);
  free(cp);
  free(cpv);
  free(cdv);

  ev = ssl2_getary(cev,2,shape);
  free(cev);
  free(shape);

  ssl2_error(icon);
  return ev;
}

static VALUE
rb_chbk2(self,zev,ind,zp,ip,dv)
     VALUE self,zev,ind,zp,ip,dv;
{
  scomplex *czev;
  int n;
  int *cind;
  int m;
  scomplex *czp;
  int *cip;
  float *cdv;
  int icon;

  int *shape,rank;

  int nn;

  czev = ssl2_getcaryc_withshape(zev,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"zev.rank!=2");
  n = shape[0];
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"zev.shape[0]!=zev.shape[1]");
  free(shape);
  cind = ssl2_getcaryi(ind,&m);
  if (m<1) rb_raise(rb_eRuntimeError,"ind.lenght<1");
  if (m>n) rb_raise(rb_eRuntimeError,"ind.lenght>n");
  czp = ssl2_getcaryc_withshape(zp,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"zp.rank!=2");
  if (shape[0]!=n) rb_raise(rb_eRuntimeError,"zp.shape[0]!=n");
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"zp.shape[1]!=n");
  cip = ssl2_getcaryi(ip,&nn);
  if (nn!=n) rb_raise(rb_eRuntimeError,"ip.length!=n");
  cdv = ssl2_getcary(dv,&nn);
  if (nn!=n) rb_raise(rb_eRuntimeError,"dv.length!=n");

  chbk2_(czev,&n,&n,cind,&m,czp,cip,cdv,&icon);
  free(cind);
  free(czp);
  free(cip);
  free(cdv);

  zev = ssl2_getaryc(czev,2,shape);
  free(czev);
  free(shape);

  ssl2_error(icon);
  return zev;
}

static VALUE
rb_trbk(self,ev,p)
     VALUE self,ev,p;
{
  float *cev;
  int n;
  int m;
  float *cp;
  int icon;

  int *shape,rank;

  int nn;

  cev = ssl2_getcary_withshape(ev,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"ev.rank != 2");
  n = shape[0];
  m = shape[1];
  if (n==1) rb_raise(rb_eRuntimeError,"n==1");
  if (n<m) rb_raise(rb_eRuntimeError,"n<m");
  cp = ssl2_getcary(p,&nn);
  if (nn!=n*(n+1)/2) rb_raise(rb_eRuntimeError,"p.length != n(n+1)/2");

  trbk_(cev,&n,&n,&m,cp,&icon);
  free(cp);

  ev = ssl2_getary(cev,2,shape);
  free(cev);
  free(shape);

  ssl2_error(icon);
  return ev;
}

static VALUE
rb_trbkh(self,evr,p,pv)
     VALUE self,evr,p,pv;
{
  float *cevr,*cevi;
  int n,m;
  float *cp,*cpv;
  int icon;

  scomplex *cev;
  VALUE ev;
  int *shape,rank;

  int i,nn;

  cevr = ssl2_getcary_withshape(evr,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"evr.rank != 2");
  n = shape[0];
  m = shape[1];
  if (n==1) rb_raise(rb_eRuntimeError,"evr.shape[0] == 1");
  if (m==0) rb_raise(rb_eRuntimeError,"evr.shape[1] == 0");
  if (n<m) rb_raise(rb_eRuntimeError,"n<m");
  cp =  ssl2_getcary_withshape(evr,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"p.rank != 2");
  if (shape[0]!=n) rb_raise(rb_eRuntimeError,"p.shape[0] != n");
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"p.shape[1] != m");
  cpv = ssl2_getcary(pv,&nn);
  if (nn!=2*n) rb_raise(rb_eRuntimeError,"pv.length != 2n");

  cevi = ALLOC_N(float,n*m);
  trbkh_(cevr,cevi,&n,&n,&m,cp,cpv,&icon);
  free(cp);
  free(cpv);

  cev = ALLOC_N(scomplex,n*m);
  for (i=0;i<n*m;i++) {
    cev[i].r = cevr[i];
    cev[i].i = cevi[i];
  }
  free(cevr);
  free(cevi);
  ev = ssl2_getaryc(cev,2,shape);
  free(cev);
  free(shape);

  ssl2_error(icon);
  return ev;
}

static VALUE
rb_nrml(self,ev,ind,mode)
     VALUE self,ev,ind,mode;
{
  float *cev;
  int n,m;
  int *cind;
  int cmode;  
  int icon;

  int *shape,rank;

  int nn;

  cev = ssl2_getcary_withshape(ev,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"ev.rank != ");
  n = shape[0];
  m = shape[1];
  if (m<1) rb_raise(rb_eRuntimeError,"m<1");
  if (n<m) rb_raise(rb_eRuntimeError,"n<m");
  cind = ssl2_getcaryi(ind,&nn);
  if (nn!=m) rb_raise(rb_eRuntimeError,"ind.lenght != m");
  cmode = NUM2INT(mode);
  if (cmode!=1&&cmode!=2) rb_raise(rb_eRuntimeError,"mode must be 1 or 2");

  nrml_(cev,&n,&n,cind,&m,&cmode,&icon);
  free(cind);

  ev = ssl2_getary(cev,2,shape);
  free(cev);
  free(shape);

  ssl2_error(icon);
  return ev;
}

static VALUE
rb_cnrml(self,zev,ind,mode)
     VALUE self,zev,ind,mode;
{
  scomplex *czev;
  int n;
  int *cind;
  int m;
  int cmode;
  int icon;

  int *shape,rank;

  czev = ssl2_getcaryc_withshape(zev,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"zev.rank!=2");
  n = shape[0];
  if (shape[1]!=n) rb_raise(rb_eRuntimeError,"zev.shape[0]!=zev.shape[1]");
  cind = ssl2_getcaryi(ind,&m);
  if (m>n) rb_raise(rb_eRuntimeError,"ind.length>n");
  if (m<1) rb_raise(rb_eRuntimeError,"ind.length<1");
  cmode = NUM2INT(mode);
  if (cmode!=1&&cmode!=2) rb_raise(rb_eRuntimeError,"mode must be 1 or 2");

  cnrml_(czev,&n,&n,cind,&m,&cmode,&icon);
  free(cind);

  zev = ssl2_getaryc(czev,2,shape);
  free(czev);
  free(shape);

  ssl2_error(icon);
  return zev;
}

static VALUE
rb_gschl(self,a,b,epsz)
     VALUE self,a,b,epsz;
{
  float *ca,*cb;
  int n;
  float cepsz;
  int icon;

  int shape[1];

  int nn;

  ca = ssl2_getcary(a,&nn);
  n = (sqrt(8*nn+1)-1)/2;
  if (nn!=n*(n+1)/2) rb_raise(rb_eRuntimeError,"a.length != n(n+1)/2");
  cb = ssl2_getcary(b,&nn);
  if (nn!=n*(n+1)/2) rb_raise(rb_eRuntimeError,"b.length != n(n+1)/2");
  cepsz = NUM2FLT(epsz);

  gschl_(ca,cb,&n,&cepsz,&icon);

  shape[0] = n*(n+1)/2;
  a = ssl2_getary(ca,1,shape);
  free(ca);
  b = ssl2_getary(cb,1,shape);
  free(cb);

  ssl2_error(icon);
  return rb_ary_new3(2,a,b);
}

static VALUE
rb_gsbk(self,ev,b)
     VALUE self,ev,b;
{
  float *cev;
  int n,m;
  float *cb;
  int icon;

  int *shape,rank;

  int nn;

  cev = ssl2_getcary_withshape(ev,&shape,&rank);
  if (rank!=2) rb_raise(rb_eRuntimeError,"ev.rank != 2");
  n = shape[0];
  m = shape[1];
  if (m==0) rb_raise(rb_eRuntimeError,"ev.shape[1]==0");
  if (m>n) rb_raise(rb_eRuntimeError,"m<n");
  cb = ssl2_getcary(b,&nn);
  if (nn!=n*(n+1)/2) rb_raise(rb_eRuntimeError,"b.length != n(n+1)/2");

  gsbk_(cev,&n,&n,&m,cb,&icon);
  free(cb);

  ev = ssl2_getary(cev,2,shape);
  free(cev);
  feee(shape);

  ssl2_error(icon);
  return ev;
}


void init_eigenvalue_eigenvector(mSSL2)
     VALUE mSSL2;
{
  rb_define_module_function(mSSL2, "eig1", rb_eig1, -1);
  rb_define_module_function(mSSL2, "ceig2", rb_ceig2, -1);
  rb_define_module_function(mSSL2, "seig1", rb_seig1, 1);
  rb_define_module_function(mSSL2, "seig2", rb_seig2, 2);
  rb_define_module_function(mSSL2, "heig2", rb_heig2, 2);
  rb_define_module_function(mSSL2, "bseg", rb_bseg, 5);
  rb_define_module_function(mSSL2, "bsegj", rb_bsegj, 6);
  rb_define_module_function(mSSL2, "teig1", rb_teig1, 2);
  rb_define_module_function(mSSL2, "teig2", rb_teig2, 3);
  rb_define_module_function(mSSL2, "gseg2", rb_gseg2, 5);
  rb_define_module_function(mSSL2, "gbseg", rb_gbseg, 8);

  rb_define_module_function(mSSL2, "hsqr", rb_hsqr, 1);
  rb_define_module_function(mSSL2, "chsqr", rb_chsqr, 1);
  rb_define_module_function(mSSL2, "trql", rb_trql, 2);
  rb_define_module_function(mSSL2, "bsct1", rb_bsct1, 4);

  rb_define_module_function(mSSL2, "hvec", rb_hvec, 3);
  rb_define_module_function(mSSL2, "chvec", rb_chvec, 3);
  rb_define_module_function(mSSL2, "bsvec", rb_bsvec, 3);

  rb_define_module_function(mSSL2, "blnc", rb_blnc, 1);
  rb_define_module_function(mSSL2, "hes1", rb_hes1, 1);
  rb_define_module_function(mSSL2, "ches2", rb_ches2, 1);
  rb_define_module_function(mSSL2, "trid1", rb_trid1, 1);
  rb_define_module_function(mSSL2, "tridh", rb_tridh, 1);
  rb_define_module_function(mSSL2, "btrid", rb_btrid, 2);
  rb_define_module_function(mSSL2, "hbk1", rb_hbk1, 5);
  rb_define_module_function(mSSL2, "chbk2", rb_chbk2, 5);
  rb_define_module_function(mSSL2, "trbk", rb_trbk, 2);
  rb_define_module_function(mSSL2, "trbkh", rb_trbkh, 3);
  rb_define_module_function(mSSL2, "nrml", rb_nrml, 3);
  rb_define_module_function(mSSL2, "cnrml", rb_cnrml, 3);
  rb_define_module_function(mSSL2, "gschl", rb_gschl, 3);
  rb_define_module_function(mSSL2, "gsbk", rb_gsbk, 2);
}