From acb57ec8dd3011534813bdbee84563f1a830f499 Mon Sep 17 00:00:00 2001
From: Maxim Koltsov <kolmax94@gmail.com>
Date: Thu, 8 Nov 2018 11:07:33 +0300
Subject: Add generalized eigenvalues via dggev and zggev

These lapack functions generalize dgeev and zgeev. Interface for them
was added similarly to eig* functions already present in hmatrix.
---
 packages/base/src/Internal/C/lapack-aux.c | 87 +++++++++++++++++++++++++++++++
 1 file changed, 87 insertions(+)

(limited to 'packages/base/src/Internal/C/lapack-aux.c')

diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 5018a45..dddbefb 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -415,6 +415,93 @@ int eig_l_R(ODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) {
     OK
 }
 
+//////////////////// generalized real eigensystem ////////////
+
+int dggev_(char *jobvl, char *jobvr, integer *n,
+    doublereal *a, integer *lda, doublereal *b, integer *ldb,
+    doublereal *alphar, doublereal *alphai, doublereal *beta,
+    doublereal *vl, integer *ldvl, doublereal *vr, integer *ldvr,
+    doublereal *work,
+	integer *lwork, integer *info);
+
+int eig_l_G(ODMAT(a), ODMAT(b), CVEC(alpha), DVEC(beta), ODMAT(vl), ODMAT(vr)) {
+    integer n = ar;
+    REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE);
+    REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE);
+    char jobvl = vlp==NULL?'N':'V';
+    REQUIRES(vrp==NULL || (vrr==n && vrc==n), BAD_SIZE);
+    char jobvr = vrp==NULL?'N':'V';
+    DEBUGMSG("eig_l_G");
+    integer lwork = -1;
+    integer res;
+    // ask for optimal lwork
+    double ans;
+    dggev_  (&jobvl,&jobvr,
+             &n,
+             ap,&n,bp,&n,
+             (double*)alphap, (double*)alphap+n, betap,
+             vlp, &n, vrp, &n,
+             &ans, &lwork,
+             &res);
+    lwork = ceil(ans);
+    double * work = (double*)malloc(lwork*sizeof(double));
+    CHECK(!work,MEM);
+    dggev_  (&jobvl,&jobvr,
+             &n,
+             ap,&n,bp,&n,
+             (double*)alphap, (double*)alphap+n, betap,
+             vlp, &n, vrp, &n,
+             work, &lwork,
+             &res);
+    CHECK(res,res);
+    free(work);
+    OK
+}
+
+//////////////////// generalized complex eigensystem ////////////
+
+int zggev_(char *jobvl, char *jobvr, integer *n,
+    doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
+    doublecomplex *alphar, doublecomplex *beta,
+    doublecomplex *vl, integer *ldvl, doublecomplex *vr, integer *ldvr,
+    doublecomplex *work, integer *lwork,
+    doublereal *rwork, integer *info);
+
+int eig_l_GC(OCMAT(a), OCMAT(b), CVEC(alpha), CVEC(beta), OCMAT(vl), OCMAT(vr)) {
+    integer n = ar;
+    REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE);
+    REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE);
+    char jobvl = vlp==NULL?'N':'V';
+    REQUIRES(vrp==NULL || (vrr==n && vrc==n), BAD_SIZE);
+    char jobvr = vrp==NULL?'N':'V';
+    DEBUGMSG("eig_l_GC");
+    double *rwork = (double*) malloc(8*n*sizeof(double));
+    CHECK(!rwork,MEM);
+    integer lwork = -1;
+    integer res;
+    // ask for optimal lwork
+    doublecomplex ans;
+    zggev_  (&jobvl,&jobvr,
+             &n,
+             ap,&n,bp,&n,
+             alphap, betap,
+             vlp, &n, vrp, &n,
+             &ans, &lwork,
+             rwork, &res);
+    lwork = ceil(ans.r);
+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    CHECK(!work,MEM);
+    zggev_  (&jobvl,&jobvr,
+             &n,
+             ap,&n,bp,&n,
+             alphap, betap,
+             vlp, &n, vrp, &n,
+             work, &lwork,
+             rwork, &res);
+    CHECK(res,res);
+    free(work);
+    OK
+}
 
 //////////////////// symmetric real eigensystem ////////////
 
-- 
cgit v1.2.3