[Dplr-commits] r848 - in pkg/dplR: . src

[noreply at r-forge.r-project.org]

Author: mvkorpel
Date: 2014-05-02 15:37:25 +0200 (Fri, 02 May 2014)
New Revision: 848

Modified:
   pkg/dplR/ChangeLog
   pkg/dplR/src/dplR.h
   pkg/dplR/src/redfit.c
Log:
redfit: Avoid using the long double type in some cases, hopefully
speeding up slow computation times on some systems
(r-patched-solaris-sparc test results on CRAN...)


Modified: pkg/dplR/ChangeLog
===================================================================
--- pkg/dplR/ChangeLog	2014-05-02 13:32:13 UTC (rev 847)
+++ pkg/dplR/ChangeLog	2014-05-02 13:37:25 UTC (rev 848)
@@ -6,6 +6,13 @@
 - Bug fix: make.plot=TRUE threw an error when input data.frame had leading
   or trailing all-NA rows
 
+File: redfit.c
+--------------
+
+- Avoid using the long double type in some cases. Hopefully this
+  will speed up otherwise unbearable computation times on some
+  systems.
+
 * CHANGES IN dplR VERSION 1.6.0
 
 File: TODO

Modified: pkg/dplR/src/dplR.h
===================================================================
--- pkg/dplR/src/dplR.h	2014-05-02 13:32:13 UTC (rev 847)
+++ pkg/dplR/src/dplR.h	2014-05-02 13:37:25 UTC (rev 848)
@@ -4,6 +4,7 @@
 #include <R.h>  /* to include Rconfig.h */
 #include <Rversion.h>
 #include <Rinternals.h>
+#include <float.h>
 size_t dplRlength(SEXP x);
           
 #ifdef ENABLE_NLS
@@ -18,4 +19,16 @@
 #define DPLR_RGEQ3
 #endif
 
+/*
+  dplr_ldouble is a 64 or 80 bit floating point type
+*/
+#if LDBL_MANT_DIG > 64
+typedef double dplr_ldouble;
+/* 64 bits */
+#else
+#define DPLR_LONG
+typedef long double dplr_ldouble;
+/* 64 or 80 bits */
 #endif
+
+#endif

Modified: pkg/dplR/src/redfit.c
===================================================================
--- pkg/dplR/src/redfit.c	2014-05-02 13:32:13 UTC (rev 847)
+++ pkg/dplR/src/redfit.c	2014-05-02 13:37:25 UTC (rev 848)
@@ -31,7 +31,7 @@
 void ftfix(const double *xx, const double *tsamp, const size_t nxx,
 	   const double *freq, const size_t nfreq, const double si,
 	   const size_t lfreq, const double tzero, const double *tcos,
-	   const double *tsin, const double *wtau, const long double sumbysqrt,
+	   const double *tsin, const double *wtau, const dplr_ldouble sumbysqrt,
 	   double *ftrx, double *ftix);
 SEXP makear1(SEXP t, SEXP np, SEXP tau);
 
@@ -165,7 +165,7 @@
 	    SEXP segskip, SEXP lmfit) {
     SEXP gxx, twk, xwk, ftrx, ftix, tmp, cbindcall, lengthfun;
     double dnseg, segskip_val, scal, np_val;
-    long double sumx, sqrt_nseg;
+    dplr_ldouble sumx, sqrt_nseg;
     size_t i, j, nseg_val, nfreq_val, n50_val, segstart, ncopy;
     size_t sincos_skip, wtau_skip;
     size_t wwidx = 0;
@@ -217,7 +217,11 @@
     xwk_data = REAL(xwk);
     ftrx_data = REAL(ftrx);
     ftix_data = REAL(ftix);
+#ifdef DPLR_LONG
     sqrt_nseg = sqrtl((long double) dnseg);
+#else
+    sqrt_nseg = sqrt(dnseg);
+#endif
     wtau_skip = nfreq_val - 1;
     sincos_skip = wtau_skip * nseg_val;
     for (i = 0; i < nfreq_val; i++) {
@@ -233,7 +237,7 @@
 	/* detrend data */
 	rmtrend(twk, xwk, lengthfun, lmfit);
         /* apply window to data */
-	sumx = 0.0L;
+	sumx = 0.0;
 	for (j = 0; j < nseg_val; j++) {
 	    xwk_data[j] *= ww_data[wwidx++];
 	    sumx += xwk_data[j];
@@ -284,14 +288,14 @@
 void ftfix(const double *xx, const double *tsamp, const size_t nxx,
 	   const double *freq, const size_t nfreq, const double si,
 	   const size_t lfreq, const double tzero, const double *tcos,
-	   const double *tsin, const double *wtau, const long double sumbysqrt,
+	   const double *tsin, const double *wtau, const dplr_ldouble sumbysqrt,
 	   double *ftrx, double *ftix) {
     const double_t tol1 = 1.0e-4;
     const double tol2 = 1.0e-8;
     const double_t const1 = M_SQRT1_2;
     double_t const2;
     double const3, ftrd, ftid, phase, wtnew, tmpsin, tmpcos, wrun;
-    long double cross, sumr, sumi, scos2, ssin2;
+    dplr_ldouble cross, sumr, sumi, scos2, ssin2;
     double complex work;
     size_t i, ii, iput;
     size_t idx = 0;
@@ -306,11 +310,11 @@
 	wrun = M_2PI * freq[ii]; /* omega = 2 * pi * freq */
     	wtnew = wtau[ii - 1];
     	/* summations over the sample */
-    	cross = 0.0L;
-    	scos2 = 0.0L;
-    	ssin2 = 0.0L;
-    	sumr = 0.0L;
-    	sumi = 0.0L;
+    	cross = 0.0;
+    	scos2 = 0.0;
+    	ssin2 = 0.0;
+    	sumr = 0.0;
+    	sumi = 0.0;
     	for (i = 0; i < nxx; i++) {
     	    tmpsin = tsin[idx];
     	    tmpcos = tcos[idx];