Skip to content

Commit

Permalink
patch source code for compile errors (Jagan/Bolker)
Browse files Browse the repository at this point in the history
  • Loading branch information
bbolker committed Mar 12, 2023
1 parent b4edd5d commit a72d4af
Show file tree
Hide file tree
Showing 3 changed files with 26 additions and 14 deletions.
8 changes: 4 additions & 4 deletions src/chebpol.c
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ static void chebcoef(double *x, int *dims, const int rank, double *F, int dct, i
src = dest;
dest = p;

F77_CALL(dgemm)("TRANS","TRANS",&N,&stride,&N,&alpha,mat[i],&N,src,&stride,&beta,dest,&N);
F77_CALL(dgemm)("T", "T", &N, &stride, &N, &alpha, mat[i], &N, src, &stride, &beta, dest, &N FCONE FCONE);
// Fix the first element in each vector, nah do it at the end
// for(int k = 0; k < siz; k+= N) dest[k] *= 0.5;
}
Expand Down Expand Up @@ -193,7 +193,7 @@ static double FH(double *fv, double *x, double **knots, int *dims, const int ran

// Special case:
for(int i = 0; i < N; i++) {
if(fabs(xx - kn[i]) < 10.0*DOUBLE_EPS) return FH(&fv[i*siz], x, knots, dims, newrank, weights);
if(fabs(xx - kn[i]) < 10.0*DBL_EPSILON) return FH(&fv[i*siz], x, knots, dims, newrank, weights);
}

#if 1
Expand Down Expand Up @@ -755,7 +755,7 @@ static double findsimplex(double *x, double *knots, int *dtri, double *lumats, i
int *ipiv = ipivs + simplex*N;
for(int d = 0; d < dim; d++) vec[d] = x[d];
vec[dim] = 1.0;
F77_CALL(dgetrs)("N", &N, &one, lumat, &N, ipiv, vec, &N, &info);
F77_CALL(dgetrs)("N", &N, &one, lumat, &N, ipiv, vec, &N, &info FCONE);
for(int d = 0; d < N; d++) if(vec[d] < -1e-10) {bad = 1; break;}
if(bad) continue;

Expand Down Expand Up @@ -794,7 +794,7 @@ static double findsimplex(double *x, double *knots, int *dtri, double *lumats, i
const int *tri = dtri + nearest * (dim+1);
for(int d = 0; d < dim; d++) vec[d] = x[d];
vec[dim] = 1.0;
F77_CALL(dgetrs)("N", &N, &one, lumat, &N, ipiv, vec, &N, &info);
F77_CALL(dgetrs)("N", &N, &one, lumat, &N, ipiv, vec, &N, &info FCONE);
double sum = 0;
for(int d = 0; d <= dim; d++) sum += vec[d]*val[tri[d]-1];
return sum;
Expand Down
20 changes: 16 additions & 4 deletions src/chebpol.h
Original file line number Diff line number Diff line change
@@ -1,14 +1,26 @@
#ifndef _CHEBPOL
#define _CHEBPOL

#include <float.h>
#include <math.h>
#include <Rmath.h>
#include <R.h>
#include <Rdefines.h>
#include <R_ext/Rdynload.h>

#ifndef USE_FC_LEN_T
# define USE_FC_LEN_T
#endif
#include <Rconfig.h>
#include <R_ext/BLAS.h>
#include <R_ext/Lapack.h>
#ifndef FCONE
# define FCONE
#endif

#include <R.h>
#include <Rdefines.h>
#include <Rmath.h>
#include <R_ext/Applic.h>
#include <R_ext/Rdynload.h>
#include <R_ext/Visibility.h>

#include "config.h"

SEXP R_makerbf(SEXP, SEXP, SEXP, SEXP);
Expand Down
12 changes: 6 additions & 6 deletions src/stalker.c
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
#include <gsl/gsl_roots.h>
#endif

#define MYEPS (1e3*sqrt(DOUBLE_EPS))

This comment has been minimized.

Copy link
@jaganmn

jaganmn Mar 12, 2023

Thanks - I must have forgotten to grep and missed these.

#define MYEPS (1e3*sqrt(DBL_EPSILON))
typedef struct {double dmin, dplus, vmin, vplus; int pos;} rblock;
typedef struct {double kq, vmin, vplus;} sblock;
typedef struct {double *det, *pmin, *pplus;int uniform;} precomp;
Expand Down Expand Up @@ -135,15 +135,15 @@ static R_INLINE double stalkhyp(double x, double vmin, double vplus,double dmin,
// double iD = 1.0/(dmin*pow(dplus,r) + pow(dmin,r)*dplus);
// double b = (vplus*pow(dmin,r) - vmin*pow(dplus,r))*iD;
// double c = (vplus*dmin + vmin*dplus)*iD;
if(fabs(x - dplus) < 1e3*DOUBLE_EPS*dplus) return vplus;
if(fabs(x + dmin) < 1e3*DOUBLE_EPS*dmin) return vmin;
if(fabs(x - dplus) < 1e3*DBL_EPSILON*dplus) return vplus;
if(fabs(x + dmin) < 1e3*DBL_EPSILON*dmin) return vmin;
// hyperbolic function
// a+bx-a/(cx+1)
if(fabs(vplus*dmin + vmin*dplus) < 1e3*DOUBLE_EPS*dmin*dplus) return x*(vplus-vmin)/(dplus+dmin); //linear
if(fabs(vplus*dmin + vmin*dplus) < 1e3*DBL_EPSILON*dmin*dplus) return x*(vplus-vmin)/(dplus+dmin); //linear
if(sign(vplus*vmin) < 0) {
//monotonic, set b=0
double a,c;
if(fabs(vmin-vplus) < 1e3*DOUBLE_EPS) return 0;
if(fabs(vmin-vplus) < 1e3*DBL_EPSILON) return 0;
#if 1
c = (vmin*dplus + vplus*dmin)/(dmin*dplus*(vmin-vplus));
a = vplus + vplus/(c*dplus);
Expand All @@ -170,7 +170,7 @@ static R_INLINE double stalkhyp(double x, double vmin, double vplus,double dmin,
// non-monotonic
double D = dplus*dplus*vmin - dmin*dmin*vplus;

if(fabs(D) < 1e3*DOUBLE_EPS) {
if(fabs(D) < 1e3*DBL_EPSILON) {
double a;
// if c==0, it's a parabola y = a*x^2
// we have (0,0), (-dmin,vmin), and (dplus,vplus)
Expand Down

0 comments on commit a72d4af

Please sign in to comment.