Add non-sym eigen value computation.
diff --git a/TensorMath.lua b/TensorMath.lua
index 50d3158..74bc89c 100644
--- a/TensorMath.lua
+++ b/TensorMath.lua
@@ -986,6 +986,18 @@
{name='charoption', values={'N', 'V'}, default='N'},
{name='charoption', values={'U', 'L'}, default='U'}}
)
+ interface:wrap("reig",
+ cname("geev"),
+ {{name=Tensor, returned=true},
+ {name=Tensor, returned=true},
+ {name=Tensor},
+ {name='charoption', values={'N', 'V'}, default='N'}},
+ cname("geev"),
+ {{name=Tensor, default=true, returned=true, invisible=true},
+ {name=Tensor, default=true, returned=true, invisible=true},
+ {name=Tensor},
+ {name='charoption', values={'N', 'V'}, default='N'}}
+ )
interface:wrap("svd",
cname("gesvd"),
diff --git a/lib/TH/generic/THLapack.c b/lib/TH/generic/THLapack.c
index 1816edb..44ac046 100644
--- a/lib/TH/generic/THLapack.c
+++ b/lib/TH/generic/THLapack.c
@@ -48,6 +48,21 @@
#endif
}
+void THLapack_(geev)(char jobvl, char jobvr, int n, real *a, int lda, real *wr, real *wi, real* vl, int ldvl, real *vr, int ldvr, real *work, int lwork, int *info)
+{
+#ifdef USE_LAPACK
+#if defined(TH_REAL_IS_DOUBLE)
+ extern void dgeev_(char *jobvl, char *jobvr, int *n, double *a, int *lda, double *wr, double *wi, double* vl, int *ldvl, double *vr, int *ldvr, double *work, int *lwork, int *info);
+ dgeev_(&jobvl, &jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, info);
+#else
+ extern void sgeev_(char *jobvl, char *jobvr, int *n, float *a, int *lda, float *wr, float *wi, float* vl, int *ldvl, float *vr, int *ldvr, float *work, int *lwork, int *info);
+ sgeev_(&jobvl, &jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, info);
+#endif
+#else
+ THError("geev : Lapack library not found in compile time\n");
+#endif
+}
+
void THLapack_(gesvd)(char jobu, char jobvt, int m, int n, real *a, int lda, real *s, real *u, int ldu, real *vt, int ldvt, real *work, int lwork, int *info)
{
#ifdef USE_LAPACK
diff --git a/lib/TH/generic/THLapack.h b/lib/TH/generic/THLapack.h
index 44eb093..3c2b4e8 100644
--- a/lib/TH/generic/THLapack.h
+++ b/lib/TH/generic/THLapack.h
@@ -10,6 +10,8 @@
void THLapack_(gels)(char trans, int m, int n, int nrhs, real *a, int lda, real *b, int ldb, real *work, int lwork, int *info);
/* Eigenvals */
void THLapack_(syev)(char jobz, char uplo, int n, real *a, int lda, real *w, real *work, int lwork, int *info);
+/* Non-sym eigenvals */
+void THLapack_(geev)(char jobvl, char jobvr, int n, real *a, int lda, real *wr, real *wi, real* vl, int ldvl, real *vr, int ldvr, real *work, int lwork, int *info);
/* svd */
void THLapack_(gesvd)(char jobu, char jobvt, int m, int n, real *a, int lda, real *s, real *u, int ldu, real *vt, int ldvt, real *work, int lwork, int *info);
/* LU decomposition */
diff --git a/lib/TH/generic/THTensorLapack.c b/lib/TH/generic/THTensorLapack.c
index 1ba025a..c8e9751 100644
--- a/lib/TH/generic/THTensorLapack.c
+++ b/lib/TH/generic/THTensorLapack.c
@@ -190,6 +190,72 @@
THTensor_(free)(work);
}
+TH_API void THTensor_(geev)(THTensor *re_, THTensor *rv_, THTensor *a_, const char *jobvr)
+{
+ int n, lda, lwork, info, ldvr;
+ THTensor *work, *wi, *wr, *a;
+ real wkopt;
+ real *rv_data;
+ long i;
+
+ THArgCheck(a_->nDimension == 2, 3, "A should be 2 dimensional");
+ THArgCheck(a_->size[0] == a_->size[1], 3,"A should be square");
+
+ /* we want to definitely clone */
+ a = THTensor_(new)();
+ THTensor_(lapackClone)(a,a_,1);
+
+ n = a->size[0];
+ lda = n;
+
+ wi = THTensor_(new)();
+ wr = THTensor_(new)();
+ THTensor_(resize2d)(re_,n,2);
+ THTensor_(resize1d)(wi,n);
+ THTensor_(resize1d)(wr,n);
+
+ rv_data = NULL;
+ ldvr = 1;
+ if (*jobvr == 'V')
+ {
+ THTensor_(resize2d)(rv_,n,n);
+ rv_data = THTensor_(data)(rv_);
+ ldvr = n;
+ }
+ // get optimal workspace size
+ THLapack_(geev)('N', jobvr[0], n, THTensor_(data)(a), lda, THTensor_(data)(wr), THTensor_(data)(wi),
+ NULL, 1, rv_data, ldvr, &wkopt, -1, &info);
+
+ lwork = (int)wkopt;
+ work = THTensor_(newWithSize1d)(lwork);
+
+ THLapack_(geev)('N', jobvr[0], n, THTensor_(data)(a), lda, THTensor_(data)(wr), THTensor_(data)(wi),
+ NULL, 1, rv_data, ldvr, THTensor_(data)(work), lwork, &info);
+
+ if (info > 0)
+ {
+ THError(" Lapack geev : Failed to converge. %d off-diagonal elements of an didn't converge to zero",info);
+ }
+ else if (info < 0)
+ {
+ THError("Lapack geev : Argument %d : illegal value", -info);
+ }
+
+ real *re_data = THTensor_(data)(re_);
+ real *wi_data = THTensor_(data)(wi);
+ real *wr_data = THTensor_(data)(wr);
+ for (i=0; i<n; i++)
+ {
+ re_data[2*i] = wr_data[i];
+ re_data[2*i+1] = wi_data[i];
+ }
+ THTensor_(transpose)(rv_,NULL,0,1);
+ THTensor_(free)(a);
+ THTensor_(free)(wi);
+ THTensor_(free)(wr);
+ THTensor_(free)(work);
+}
+
TH_API void THTensor_(syev)(THTensor *re_, THTensor *rv_, THTensor *a, const char *jobz, const char *uplo)
{
int n, lda, lwork, info;
diff --git a/lib/TH/generic/THTensorLapack.h b/lib/TH/generic/THTensorLapack.h
index de7715e..1b17907 100644
--- a/lib/TH/generic/THTensorLapack.h
+++ b/lib/TH/generic/THTensorLapack.h
@@ -5,6 +5,7 @@
TH_API void THTensor_(gesv)(THTensor *rb_, THTensor *ra_, THTensor *b_, THTensor *a_);
TH_API void THTensor_(gels)(THTensor *rb_, THTensor *ra_, THTensor *b_, THTensor *a_);
TH_API void THTensor_(syev)(THTensor *re_, THTensor *rv_, THTensor *a_, const char *jobz, const char *uplo);
+TH_API void THTensor_(geev)(THTensor *re_, THTensor *rv_, THTensor *a_, const char *jobvr);
TH_API void THTensor_(gesvd)(THTensor *ru_, THTensor *rs_, THTensor *rv_, THTensor *a, const char *jobu);
TH_API void THTensor_(gesvd2)(THTensor *ru_, THTensor *rs_, THTensor *rv_, THTensor *ra_, THTensor *a, const char *jobu);
TH_API void THTensor_(getri)(THTensor *ra_, THTensor *a);
