From b929089586ee6da5ab92e6890b68876a598112bc Mon Sep 17 00:00:00 2001
From: Fidler <matthew.fidler@gmail.com>
Date: Tue, 19 Jun 2018 11:30:46 -0500
Subject: [PATCH] Split out thread safe into blas0.  See Issue #15

---
 src/blas.h         | 18 +++++++++---------
 src/correction.c   |  4 ++--
 src/daxpy.c        |  2 +-
 src/ddot.c         |  2 +-
 src/dgefa.c        |  8 ++++----
 src/dgesl.c        | 10 +++++-----
 src/dscal.c        |  2 +-
 src/fnorm.c        |  2 +-
 src/idamax.c       |  2 +-
 src/lsoda.c        | 20 ++++++++++----------
 src/methodswitch.c |  4 ++--
 src/orderswitch.c  |  2 +-
 src/prja.c         |  4 ++--
 src/stoda.c        |  6 +++---
 src/vmnorm.c       |  2 +-
 15 files changed, 44 insertions(+), 44 deletions(-)

diff --git a/src/blas.h b/src/blas.h
index e6426f2..a405b86 100644
--- a/src/blas.h
+++ b/src/blas.h
@@ -1,14 +1,14 @@
 
-double ddot(int n, double dx[], int incx, double dy[], int incy);
-void dgesl(double **a, int n, int * ipvt, double b[], int job);
-void dgefa(double ** a, int n, int * ipvt, int * info);
-void daxpy(int n, double da, double dx[], int incx, double dy[], int incy);
-int idamax(int n, double dx[], int incx);
-void dscal(int n, double da, double dx[], int incx);
-double vmnorm(int n, double *v, double *w);
-double fnorm(int n, double **a, double *w);
+double ddot0(int n, double dx[], int incx, double dy[], int incy);
+void dgesl0(double **a, int n, int * ipvt, double b[], int job);
+void dgefa0(double ** a, int n, int * ipvt, int * info);
+void daxpy0(int n, double da, double dx[], int incx, double dy[], int incy);
+int idamax0(int n, double dx[], int incx);
+void dscal0(int n, double da, double dx[], int incx);
+double vmnorm0(int n, double *v, double *w);
+double fnorm0(int n, double **a, double *w);
 #if 0
-static double vmnorm(int n, double *v, double *w)
+static double vmnorm0(int n, double *v, double *w)
 
 /*
    This function routine computes the weighted max-norm
diff --git a/src/correction.c b/src/correction.c
index b95bf41..bd1ad14 100644
--- a/src/correction.c
+++ b/src/correction.c
@@ -61,7 +61,7 @@ int correction(struct lsoda_context_t * ctx, double *y, double pnorm, double *de
 				_C(savf)[i] = _C(h) * _C(savf)[i] - _C(yh)[2][i];
 				y[i] = _C(savf)[i] - _C(acor)[i];
 			}
-			*del = vmnorm(neq, y, _C(ewt));
+			*del = vmnorm0(neq, y, _C(ewt));
 			for (i = 1; i <= neq; i++) {
 				y[i] = _C(yh)[1][i] + _C(el)[1] * _C(savf)[i];
 				_C(acor)[i] = _C(savf)[i];
@@ -77,7 +77,7 @@ int correction(struct lsoda_context_t * ctx, double *y, double pnorm, double *de
 			for (i = 1; i <= neq; i++)
 				y[i] = _C(h) * _C(savf)[i] - (_C(yh)[2][i] + _C(acor)[i]);
 			solsy(ctx, y);
-			*del = vmnorm(neq, y, _C(ewt));
+			*del = vmnorm0(neq, y, _C(ewt));
 			for (i = 1; i <= neq; i++) {
 				_C(acor)[i] += y[i];
 				y[i] = _C(yh)[1][i] + _C(el)[1] * _C(acor)[i];
diff --git a/src/daxpy.c b/src/daxpy.c
index 5609f82..8a6e226 100644
--- a/src/daxpy.c
+++ b/src/daxpy.c
@@ -11,7 +11,7 @@ To: whitbeck@sanjuan.wrc.unr.edu
 */
 
 void 
-daxpy(n, da, dx, incx, dy, incy)
+daxpy0(n, da, dx, incx, dy, incy)
 	double          da, *dx, *dy;
 	int             n, incx, incy;
 
diff --git a/src/ddot.c b/src/ddot.c
index 85c88eb..7b88a11 100644
--- a/src/ddot.c
+++ b/src/ddot.c
@@ -3,7 +3,7 @@
  **********/
 
 double 
-ddot(n, dx, incx, dy, incy)
+ddot0(n, dx, incx, dy, incy)
 	double         *dx, *dy;
 	int             n, incx, incy;
 
diff --git a/src/dgefa.c b/src/dgefa.c
index 72a965b..dc8e472 100644
--- a/src/dgefa.c
+++ b/src/dgefa.c
@@ -5,7 +5,7 @@
  ***********/
 
 void 
-dgefa(a, n, ipvt, info)
+dgefa0(a, n, ipvt, info)
 	double        **a;
 	int             n, *ipvt, *info;
 
@@ -60,7 +60,7 @@ dgefa(a, n, ipvt, info)
    Find j = pivot index.  Note that a[k]+k-1 is the address of
    the 0-th element of the row vector whose 1st element is a[k][k].
 */
-		j = idamax(n - k + 1, a[k] + k - 1, 1) + k - 1;
+		j = idamax0(n - k + 1, a[k] + k - 1, 1) + k - 1;
 		ipvt[k] = j;
 /*
    Zero pivot implies this row already triangularized.
@@ -81,7 +81,7 @@ dgefa(a, n, ipvt, info)
    Compute multipliers.
 */
 		t = -1. / a[k][k];
-		dscal(n - k, t, a[k] + k, 1);
+		dscal0(n - k, t, a[k] + k, 1);
 /*
    Column elimination with row indexing.
 */
@@ -91,7 +91,7 @@ dgefa(a, n, ipvt, info)
 				a[i][j] = a[i][k];
 				a[i][k] = t;
 			}
-			daxpy(n - k, t, a[k] + k, 1, a[i] + k, 1);
+			daxpy0(n - k, t, a[k] + k, 1, a[i] + k, 1);
 		}
 	}			/* end k-loop  */
 
diff --git a/src/dgesl.c b/src/dgesl.c
index df77c91..6bc7c12 100644
--- a/src/dgesl.c
+++ b/src/dgesl.c
@@ -4,7 +4,7 @@
  ***********/
 
 void 
-dgesl(a, n, ipvt, b, job)
+dgesl0(a, n, ipvt, b, job)
 	double        **a, *b;
 	int             n, *ipvt, job;
 
@@ -58,14 +58,14 @@ dgesl(a, n, ipvt, b, job)
    First solve L * y = b.
 */
 		for (k = 1; k <= n; k++) {
-			t = ddot(k - 1, a[k], 1, b, 1);
+			t = ddot0(k - 1, a[k], 1, b, 1);
 			b[k] = (b[k] - t) / a[k][k];
 		}
 /*
    Now solve U * x = y.
 */
 		for (k = n - 1; k >= 1; k--) {
-			b[k] = b[k] + ddot(n - k, a[k] + k, 1, b + k, 1);
+			b[k] = b[k] + ddot0(n - k, a[k] + k, 1, b + k, 1);
 			j = ipvt[k];
 			if (j != k) {
 				t = b[j];
@@ -87,7 +87,7 @@ dgesl(a, n, ipvt, b, job)
 			b[j] = b[k];
 			b[k] = t;
 		}
-		daxpy(n - k, t, a[k] + k, 1, b + k, 1);
+		daxpy0(n - k, t, a[k] + k, 1, b + k, 1);
 	}
 /*
    Now solve Transpose(L) * x = y.
@@ -95,7 +95,7 @@ dgesl(a, n, ipvt, b, job)
 	for (k = n; k >= 1; k--) {
 		b[k] = b[k] / a[k][k];
 		t = -b[k];
-		daxpy(k - 1, t, a[k], 1, b, 1);
+		daxpy0(k - 1, t, a[k], 1, b, 1);
 	}
 
 }
diff --git a/src/dscal.c b/src/dscal.c
index 6bde224..442e28f 100644
--- a/src/dscal.c
+++ b/src/dscal.c
@@ -4,7 +4,7 @@
  ***********/
 
 void 
-dscal(n, da, dx, incx)
+dscal0(n, da, dx, incx)
 	double          da, *dx;
 	int             n, incx;
 
diff --git a/src/fnorm.c b/src/fnorm.c
index 1ab2909..a21ec0c 100644
--- a/src/fnorm.c
+++ b/src/fnorm.c
@@ -1,6 +1,6 @@
 #include <math.h>
 
-double fnorm(int n, double **a, double *w)
+double fnorm0(int n, double **a, double *w)
 
 /*
    This subroutine computes the norm of a full n by n matrix,
diff --git a/src/idamax.c b/src/idamax.c
index 60e903a..7be7d15 100644
--- a/src/idamax.c
+++ b/src/idamax.c
@@ -5,7 +5,7 @@
 #include <math.h>
 
 int 
-idamax(n, dx, incx)
+idamax0(n, dx, incx)
 	double         *dx;
 	int             n, incx;
 
diff --git a/src/lsoda.c b/src/lsoda.c
index 58f4e7d..3d869f2 100644
--- a/src/lsoda.c
+++ b/src/lsoda.c
@@ -118,7 +118,7 @@ tam@wri.com
 /* Terminate lsoda due to various error conditions. */
 #define softfailure(code, fmt,...) \
 { \
-	int             i; \
+	int i=0; \
 	int neq = ctx->neq; \
  \
 	ERROR(fmt, __VA_ARGS__); \
@@ -131,7 +131,7 @@ tam@wri.com
 
 #define softfailure0(code, fmt,...)              \
   {						 \
-    int             i;				 \
+    int i=0;				 \
     int neq = ctx->neq;				 \
 						 \
     ERROR0(fmt);				 \
@@ -151,7 +151,7 @@ tam@wri.com
 
 #define successreturn() \
 { \
-	int             i; \
+	int i=0; \
 	int neq = ctx->neq; \
  \
 	for (i = 1; i <= neq; i++) \
@@ -202,7 +202,7 @@ static int check_opt(struct lsoda_context_t * ctx, struct lsoda_opt_t * opt) {
 	 */
 	if (ctx->state == 1 || ctx->state == 3) {
 		/* c convention starts from 0. converted fortran code expects 1 */
-		int i;
+		int i=0;
 		for (i = 1; i <= ctx->neq; i++) {
 			double rtoli = rtol[i];
 			double atoli = atol[i];
@@ -274,7 +274,7 @@ static int alloc_mem(struct lsoda_context_t * ctx) {
 	int nyh = ctx->neq;
 	int lenyh = 1 + max(ctx->opt->mxordn, ctx->opt->mxords);
 	long offset = 0;
-	int i;
+	int i=0;
 	long yhoff = offset;
 	/* _C(yh) */
 	offset += (1 + lenyh) * sizeof(double *);
@@ -487,7 +487,7 @@ void lsoda_reset(struct lsoda_context_t * ctx) {
 void lsoda_free(struct lsoda_context_t * ctx) {
 	free(ctx->common->memory);
 	if(ctx->error) {
-		fprintf(stderr, "unhandled error message: %s\n", ctx->error);
+		REprintf("unhandled error message: %s\n", ctx->error);
 		free(ctx->error);
 	}
 	free(ctx->common);
@@ -495,7 +495,7 @@ void lsoda_free(struct lsoda_context_t * ctx) {
 
 #define ewset(ycur)  \
 { \
-	int             i; \
+	int i=0; \
  \
 	for (i = 1; i <= neq; i++) \
 		_C(ewt)[i] = rtol[i] * fabs((ycur)[i]) + atol[i]; \
@@ -521,7 +521,7 @@ int lsoda(struct lsoda_context_t * ctx, double *y, double *t, double tout) {
          * in C y[] starts from 0, but the converted fortran code starts from 1 */
 		y--;
 
-		int             i=0, ihit;
+		int             i=0, ihit=0;
 		const int neq = ctx->neq;
 		double          big, h0, hmx, rh, tcrit, tdist, tnext, tol,
 						tolsf, tp, size, sum, w0;
@@ -656,7 +656,7 @@ int lsoda(struct lsoda_context_t * ctx, double *y, double *t, double tout) {
 				}
 				tol = fmax(tol, 100. * ETA);
 				tol = fmin(tol, 0.001);
-				sum = vmnorm(neq, _C(yh)[2], _C(ewt));
+				sum = vmnorm0(neq, _C(yh)[2], _C(ewt));
 				sum = 1. / (tol * w0 * w0) + tol * sum * sum;
 				h0 = 1. / sqrt(sum);
 				h0 = fmin(h0, tdist);
@@ -754,7 +754,7 @@ int lsoda(struct lsoda_context_t * ctx, double *y, double *t, double tout) {
 					}
 				}
 			}
-			tolsf = ETA * vmnorm(neq, _C(yh)[1], _C(ewt));
+			tolsf = ETA * vmnorm0(neq, _C(yh)[1], _C(ewt));
 			if (tolsf > 0.01) {
 				tolsf = tolsf * 200.;
 				if (_C(nst) == 0) {
diff --git a/src/methodswitch.c b/src/methodswitch.c
index 185a13c..aad1e46 100644
--- a/src/methodswitch.c
+++ b/src/methodswitch.c
@@ -70,7 +70,7 @@ void methodswitch(struct lsoda_context_t * ctx, double dsm, double pnorm, double
 				lm2 = mxords + 1;
 				exm2 = 1. / (double) lm2;
 				lm2p1 = lm2 + 1;
-				dm2 = vmnorm(neq, _C(yh)[lm2p1], _C(ewt)) / cm2[mxords];
+				dm2 = vmnorm0(neq, _C(yh)[lm2p1], _C(ewt)) / cm2[mxords];
 				rh2 = 1. / (1.2 * pow(dm2, exm2) + 0.0000012);
 			} else {
 				dm2 = dsm * (cm1[_C(nq)] / cm2[_C(nq)]);
@@ -107,7 +107,7 @@ void methodswitch(struct lsoda_context_t * ctx, double dsm, double pnorm, double
 		lm1 = mxordn + 1;
 		exm1 = 1. / (double) lm1;
 		lm1p1 = lm1 + 1;
-		dm1 = vmnorm(neq, _C(yh)[lm1p1], _C(ewt)) / cm1[mxordn];
+		dm1 = vmnorm0(neq, _C(yh)[lm1p1], _C(ewt)) / cm1[mxordn];
 		rh1 = 1. / (1.2 * pow(dm1, exm1) + 0.0000012);
 	} else {
 		dm1 = dsm * (cm2[_C(nq)] / cm1[_C(nq)]);
diff --git a/src/orderswitch.c b/src/orderswitch.c
index ee89051..4a940a2 100644
--- a/src/orderswitch.c
+++ b/src/orderswitch.c
@@ -31,7 +31,7 @@ int orderswitch(struct lsoda_context_t * ctx, double rhup, double dsm, double *r
 
 	rhdn = 0.;
 	if (_C(nq) != 1) {
-		ddn = vmnorm(neq, _C(yh)[(_C(nq) + 1)], _C(ewt)) / _C(tesco)[_C(nq)][1];
+		ddn = vmnorm0(neq, _C(yh)[(_C(nq) + 1)], _C(ewt)) / _C(tesco)[_C(nq)][1];
 		exdn = 1. / (double) _C(nq);
 		rhdn = 1. / (1.3 * pow(ddn, exdn) + 0.0000013);
 	}
diff --git a/src/prja.c b/src/prja.c
index 101d0dc..80339dc 100644
--- a/src/prja.c
+++ b/src/prja.c
@@ -49,7 +49,7 @@ int prja(struct lsoda_context_t * ctx, double *y)
 /*
    Compute norm of Jacobian.
 */
-		_C(pdnorm) = fnorm(neq, _C(wm), _C(ewt)) / fabs(hl0);
+		_C(pdnorm) = fnorm0(neq, _C(wm), _C(ewt)) / fabs(hl0);
 /*
    Add identity matrix.
 */
@@ -58,7 +58,7 @@ int prja(struct lsoda_context_t * ctx, double *y)
 /*
    Do LU decomposition on P.
 */
-		dgefa(_C(wm), neq, _C(ipvt), &ier);
+		dgefa0(_C(wm), neq, _C(ipvt), &ier);
 		if (ier != 0)
 		return 0;
 	}
diff --git a/src/stoda.c b/src/stoda.c
index 2e2e70d..7bbf181 100644
--- a/src/stoda.c
+++ b/src/stoda.c
@@ -177,7 +177,7 @@ int stoda(struct lsoda_context_t * ctx, double *y, int jstart)
 					for (i = 1; i <= neq; i++)
 						_C(yh)[i1][i] += _C(yh)[i1 + 1][i];
 				}
-			pnorm = vmnorm(neq, _C(yh)[1], _C(ewt));
+			pnorm = vmnorm0(neq, _C(yh)[1], _C(ewt));
 
 			int corflag = correction(ctx, y, pnorm, &del, &delp, told, &m);
 			if (corflag == 0)
@@ -200,7 +200,7 @@ int stoda(struct lsoda_context_t * ctx, double *y, int jstart)
 		if (m == 0)
 			dsm = del / _C(tesco)[_C(nq)][2];
 		if (m > 0)
-			dsm = vmnorm(neq, _C(acor), _C(ewt)) / _C(tesco)[_C(nq)][2];
+			dsm = vmnorm0(neq, _C(acor), _C(ewt)) / _C(tesco)[_C(nq)][2];
 		if (dsm <= 1.) {
 /*
    After a successful step, update the _C(yh) array.
@@ -245,7 +245,7 @@ int stoda(struct lsoda_context_t * ctx, double *y, int jstart)
 				if ((_C(nq) + 1) != maxord + 1) {
 					for (i = 1; i <= neq; i++)
 						_C(savf)[i] = _C(acor)[i] - _C(yh)[maxord + 1][i];
-					dup = vmnorm(neq, _C(savf), _C(ewt)) / _C(tesco)[_C(nq)][3];
+					dup = vmnorm0(neq, _C(savf), _C(ewt)) / _C(tesco)[_C(nq)][3];
 					exup = 1. / (double) ((_C(nq) + 1) + 1);
 					rhup = 1. / (1.4 * pow(dup, exup) + 0.0000014);
 				}
diff --git a/src/vmnorm.c b/src/vmnorm.c
index 9e5b5ae..8faf2db 100644
--- a/src/vmnorm.c
+++ b/src/vmnorm.c
@@ -1,6 +1,6 @@
 #include <math.h>
 
-double vmnorm(int n, double *v, double *w)
+double vmnorm0(int n, double *v, double *w)
 
 /*
    This function routine computes the weighted max-norm