Antonio Galea: reduced computational time to 30-40% by skipping unnecessary computations

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@36018 15284696-431f-4ddb-bdfa-cd5b030d7da7

vector/v.vol.rst/user.h

@@ -14,6 +14,7 @@ double crs();
 double crsd();
 double crsdr2();
 double crsdd();
+void crs_full();
 int COGRR1();
 int POINT();
 int LINEQS();

vector/v.vol.rst/user3.c

@@ -235,6 +235,7 @@ COGRR1(double x_or, double y_or, double z_or, int n_rows, int n_cols,
     double wm, dx, dy, dz, dxx, dyy, dxy, dxz, dyz, dzz, h, bmgd1,
 	bmgd2, etar, zcon, r, ww, wz, r2, hcell, zzcell2,
 	etarcell, rcell, wwcell, zzcell;
+    double x_crs,x_crsd,x_crsdd,x_crsdr2;
     int n1, k1, k2, k, i1, l, l1, n4, n5, m, i;
     int NGST, LSIZE, ngstc, nszc, ngstr, nszr, ngstl, nszl;
     int POINT();
@@ -244,6 +245,8 @@ COGRR1(double x_or, double y_or, double z_or, int n_rows, int n_cols,
     int bmask = 1;
     static FCELL *cell = NULL;
 
+    int cond1 = (gradient != NULL) || (aspect1 != NULL) || (aspect2 != NULL);
+    int cond2 = (ncurv != NULL) || (gcurv != NULL) || (mcurv != NULL);
 
 #define CEULER .57721566
     /*
@@ -450,25 +453,31 @@ COGRR1(double x_or, double y_or, double z_or, int n_rows, int n_cols,
 			    r = sqrt(r2);
 			    etar = (fi * r) / 2.;
 
-			    h = h + b[m] * crs(etar);
-			    /* partial derivatives */
-			    bmgd1 = b[m] * crsd(etar, fi);
+                            crs_full(
+                              etar,fi,
+                              &x_crs,
+                              cond1?&x_crsd:NULL,
+                              cond2?&x_crsdr2:NULL,
+                              cond2?&x_crsdd:NULL
+                            );
+                            h = h + b[m] * x_crs;
+                            if(cond1)
+                            {
+                                   bmgd1 = b[m] * x_crsd;
 			    dx = dx + bmgd1 * xx;
 			    dy = dy + bmgd1 * w[m];
-			    /* dz = dz + bmgd1 * sqrt(wz2[m]); */
 			    dz = dz + bmgd1 * wz1[m];
-			    bmgd2 = b[m] * crsdd(etar, fi);
-			    bmgd1 = b[m] * crsdr2(etar, fi);
-			    dxx = dxx + bmgd2 * xx2 + bmgd1 * xx2;
+                            }
+                            if(cond2)
+                            {
+                                   bmgd2 = b[m] * x_crsdd;
+                                   bmgd1 = b[m] * x_crsdr2;
 			    dyy = dyy + bmgd2 * w2[m] + bmgd1 * w2[m];
 			    dzz = dzz + bmgd2 * wz2[m] + bmgd1 * wz2[m];
 			    dxy = dxy + bmgd2 * xx * w[m] + bmgd1 * xx * w[m];
-			    dxz =
-				dxz + bmgd2 * xx * wz1[m] +
-				bmgd1 * xx * wz1[m];
-			    dyz =
-				dyz + bmgd2 * w[m] * wz1[m] +
-				bmgd1 * w[m] * wz1[m];
+                                   dxz = dxz + bmgd2 * xx * wz1[m] + bmgd1 * xx * wz1[m];
+                                   dyz = dyz + bmgd2 * w[m] * wz1[m] + bmgd1 * w[m] * wz1[m];
+                            }                            
 			}
 			ww = h + wmin;
 			if ((cellinp != NULL) && (cellout != NULL) &&

vector/v.vol.rst/user4.c

@@ -145,7 +145,7 @@ double erfr(double rf2)
     static double a[5] = { 0.254829592, -0.284496736, 1.421413741,
 	-1.453152027, 1.061405429
     };
-    double p = 0.3275911;
+    static double p = 0.3275911;
     double erf, t;
 
 
@@ -167,21 +167,21 @@ double crs(double x)
  */
 {
 
-    static double a[10];
-    double tp = 1.1283791671 / 2.;
+    static double a[10] = {
+      1.,
+      -1. / 3.,
+      1. / 10.,
+      -1. / 42.,
+      1. / (24. * 9.),
+      -1. / (120. * 11.),
+      1. / (720. * 13.),
+      -1. / (5040. * 15.),
+      1. / (40320. * 17.),
+      -1. / (362880. * 19.),
+    };
+    static double tp = 1.1283791671 / 2.;
     double xx, res;
 
-    a[0] = 1.;
-    a[1] = -1. / 3.;
-    a[2] = 1. / 10.;
-    a[3] = -1. / 42.;
-    a[4] = 1. / (24. * 9.);
-    a[5] = -1. / (120. * 11.);
-    a[6] = 1. / (720. * 13.);
-    a[7] = -1. / (5040. * 15.);
-    a[8] = 1. / (40320. * 17.);
-    a[9] = -1. / (362880. * 19.);
-
     if (x < XA) {
 	xx = x * x;
 	res =
@@ -204,37 +204,24 @@ double crs(double x)
 }
 
 
+void crs_full(double x, double fi, double *crs, double *crsd, double *crsdr2, double *crsdd){
+    static double a[10] = { 1., -1. / 3., 1. / 10., -1. / 42., 1. / (24. * 9.), -1. / (120. * 11.),
+      1. / (720. * 13.), -1. / (5040. * 15.), 1. / (40320. * 17.), -1. / (362880. * 19.) };
+    static double b[10] = { 0, -2. / 3., 4. / 10., -6. / 42., 8. / (24. * 9.), -10. / (120. * 11.),
+      12. / (720. * 13.), -14. / (5040. * 15.), 16. / (40320. * 17.), -18. / (362880. * 19.) };
+    static double c[10] = { 0, 0, 8. / 10., -24. / 42., 48. / (24. * 9.), -80. / (120. * 11.),
+      120. / (720. * 13.), -168. / (5040. * 15.), 16. * 14. / (40320. * 17.), -18. * 16. / (362880. * 19.) };
+    static double tp = 1.1283791671 / 2.;
+    double xx, r, r2, fi2, fi4, fi8, tmp1, tmp2;
 
-
-double crsd(double x, double fi)
-/*
-   function for first derivatives 
- */
-{
-
-    static double a[10];
-    double tp = 1.1283791671 / 2.;
-    double xx, r, r2, fi2, fi4, res;
-
-    a[1] = -2. / 3.;
-    a[2] = 4. / 10.;
-    a[3] = -6. / 42.;
-    a[4] = 8. / (24. * 9.);
-    a[5] = -10. / (120. * 11.);
-    a[6] = 12. / (720. * 13.);
-    a[7] = -14. / (5040. * 15.);
-    a[8] = 16. / (40320. * 17.);
-    a[9] = -18. / (362880. * 19.);
-
-    r = 2. * x / fi;
-    r2 = r * r;
     fi2 = fi / 2.;
     fi4 = fi2 * fi2;
+    fi8 = fi4 * fi4;
 
     if (x < XA) {
 	xx = x * x;
-	res =
-	    tp * fi4 * (a[1] +
+	*crs =
+	    tp * xx * (a[1] +
 			xx * (a[2] +
 			      xx * (a[3] +
 				    xx * (a[4] +
@@ -243,111 +230,52 @@ double crsd(double x, double fi)
 						      xx * (a[7] +
 							    xx * (a[8] +
 								  xx *
-								  a
-								  [9]))))))));
-    }
-    else {
-	res = (2. * tp * exp(-x * x) - (erf(x) / x)) / (2. * r2);
-    }
-    return (res);
-}
-
-double crsdr2(double x, double fi)
-/*
-   function for first derivatives
- */
-{
-
-    static double a[10];
-    double tp = 1.1283791671 / 2.;
-    double xx, r, r2, fi2, fi4, res;
-
-    a[1] = -2. / 3.;
-    a[2] = 4. / 10.;
-    a[3] = -6. / 42.;
-    a[4] = 8. / (24. * 9.);
-    a[5] = -10. / (120. * 11.);
-    a[6] = 12. / (720. * 13.);
-    a[7] = -14. / (5040. * 15.);
-    a[8] = 16. / (40320. * 17.);
-    a[9] = -18. / (362880. * 19.);
-
-    r = 2. * x / fi;
-    r2 = r * r;
-    fi2 = fi / 2.;
-    fi4 = fi2 * fi2;
-
-    if (x < XA) {
-	xx = x * x;
-	res =
-	    tp * fi4 * (a[1] +
-			xx * (a[2] +
-			      xx * (a[3] +
-				    xx * (a[4] +
-					  xx * (a[5] +
-						xx * (a[6] +
-						      xx * (a[7] +
-							    xx * (a[8] +
+                                                                      a[9]))))))));
+	if(crsd!=NULL){
+          *crsd =
+              tp * fi4 * (b[1] +
+                          xx * (b[2] +
+                                xx * (b[3] +
+                                      xx * (b[4] +
+                                            xx * (b[5] +
+                                                  xx * (b[6] +
+                                                        xx * (b[7] +
+                                                              xx * (b[8] +
 								  xx *
-								  a
-								  [9]))))))));
+                                                                    b[9]))))))));
     }
-    else {
-	res = (2. * tp * exp(-x * x) - (erf(x) / x)) / (2. * r2 * r2);
+	if(crsdr2!=NULL){
+	  *crsdr2 = *crsd;
     }
-    return (res);
-}
-
-double crsdd(double x, double fi)
-/*
-   function for second derivatives 
- */
-{
-
-    static double a[10];
-    double tp = 1.1283791671 / 2.;
-    double xx, r, r2, fi8, res;
-
-    a[2] = 8. / 10.;
-    a[3] = -24. / 42.;
-    a[4] = 48. / (24. * 9.);
-    a[5] = -80. / (120. * 11.);
-    a[6] = 120. / (720. * 13.);
-    a[7] = -168. / (5040. * 15.);
-    a[8] = 16. * 14. / (40320. * 17.);
-    a[9] = -18. * 16. / (362880. * 19.);
-
+	if(crsdd!=NULL){
+          *crsdd =
+              tp * fi8 * (c[2] +
+                          xx * (c[3] +
+                                xx * (c[4] +
+                                      xx * (c[5] +
+                                            xx * (c[6] +
+                                                  xx * (c[7] +
+                                                        xx * (c[8] +
+                                                              xx * c[9])))))));
+        }
+    } else {
+        tmp1 = erf(x) / x;
+	*crs = tmp1 / 2. - tp;
+        if((crsd!=NULL)||(crsdr2!=NULL)||(crsdd!=NULL)){
     r = 2. * x / fi;
     r2 = r * r;
-    fi8 = fi * fi * fi * fi / 16.;
-
-    if (x < XA) {
-	xx = x * x;
-	res =
-	    tp * fi8 * (a[2] +
-			xx * (a[3] +
-			      xx * (a[4] +
-				    xx * (a[5] +
-					  xx * (a[6] +
-						xx * (a[7] +
-						      xx * (a[8] +
-							    xx * a[9])))))));
+          tmp2 = exp(-x * x);
     }
-    else {
+        if(crsd!=NULL){ *crsd   = (2. * tp * tmp2 - tmp1) / (2. * r2); }
+	if(crsdr2!=NULL){ *crsdr2 = *crsd / r2; }
 	/* Jarov vzorec */
-	res =
-	    (erf(x) / (x * r2) -
-	     exp(-x * x) * (2 / r2 + fi * fi / 2) * tp) / (r2 * r2);;
+	if(crsdd!=NULL){ *crsdd  = (tmp1 / r2 - tmp2 * (2 / r2 + fi * fi / 2) * tp) / (r2 * r2); }
     }
-    return (res);
+    return;
 }
 
 
 
-
-
-
-
 /*********solution of system of lin. equations*********/
 
 int LINEQS(int DIM1, int N1, int N2, int *NERROR, double *DETERM)