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.time-shift
parent
acc2317953
commit
b81fcbeeb1
|
@ -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);
|
||||
}
|
||||
|
|
Loading…
Reference in New Issue