From b81fcbeeb1418765fdca0edc2b1923a24efb646c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Charybdis=29?=
 <stephane.adjemian@univ-lemans.fr>
Date: Wed, 12 Jun 2013 09:46:09 +0200
Subject: [PATCH] Added an option for the threshold level of the 0/0
 generalized eigenvalue test.

The option is passed in the fourth input argument of mjdgges. If the
mex is called with less than four arguments, then the threshold level
takes its previous (hardcoded) default value: 1e-6.
---
 mex/sources/mjdgges/mjdgges.c | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/mex/sources/mjdgges/mjdgges.c b/mex/sources/mjdgges/mjdgges.c
index 145439bb9..8d15061fc 100644
--- a/mex/sources/mjdgges/mjdgges.c
+++ b/mex/sources/mjdgges/mjdgges.c
@@ -32,7 +32,7 @@ my_criteria(const double *alphar, const double *alphai, const double *beta)
 }
 
 void
-mjdgges(double *a, double *b, double *z, double *n, double *sdim, double *eval_r, double *eval_i, double *info)
+mjdgges(double *a, double *b, double *z, double *n, double *sdim, double *eval_r, double *eval_i, double zhreshold, double *info)
 {
   lapack_int i_n, i_info, i_sdim, lwork;
   double *alphar, *alphai, *beta, *work, *par, *pai, *pb, *per, *pei;
@@ -60,7 +60,7 @@ mjdgges(double *a, double *b, double *z, double *n, double *sdim, double *eval_r
   pei = eval_i;
   for (per = eval_r; per <= &eval_r[i_n-1]; ++per)
     {
-      if ((fabs(*par) > 1e-6) || (fabs(*pb) > 1e-6)) 
+      if ((fabs(*par) > zhreshold) || (fabs(*pb) > zhreshold)) 
 	*per = *par / *pb;
       else
 	{
@@ -136,6 +136,18 @@ mexFunction(int nlhs, mxArray *plhs[],
       criterium = 1+1e-6;
     }
 
+  /* set criterium for 0/0 generalized eigenvalues */
+  double zhreshold;
+  if (nrhs == 4 && mxGetM(prhs[3]) > 0)
+    {
+      zhreshold = *mxGetPr(prhs[3]);
+    }
+  else
+    {
+      zhreshold = 1e-6;
+    }
+
+
   /* keep a and b intact */
   memcpy(s, a, sizeof(double)*n1*n1);
   memcpy(t, b, sizeof(double)*n1*n1);
@@ -143,7 +155,7 @@ mexFunction(int nlhs, mxArray *plhs[],
   n = n1;
 
   /* Do the actual computations in a subroutine */
-  mjdgges(s, t, z, &n, sdim, eval_r, eval_i, info);
+  mjdgges(s, t, z, &n, sdim, eval_r, eval_i, zhreshold, info);
 
   plhs[0] = mxCreateDoubleScalar(0);
 }