[GRASS-SVN] r49653 - grass/trunk/lib/gmath

[svn_grass at osgeo.org]

Author: hamish
Date: 2011-12-11 17:51:42 -0800 (Sun, 11 Dec 2011)
New Revision: 49653

Modified:
   grass/trunk/lib/gmath/lu.c
Log:
OpenMP pragma useful for v.surf.rst, doxygenize

Modified: grass/trunk/lib/gmath/lu.c
===================================================================
--- grass/trunk/lib/gmath/lu.c	2011-12-12 01:37:38 UTC (rev 49652)
+++ grass/trunk/lib/gmath/lu.c	2011-12-12 01:51:42 UTC (rev 49653)
@@ -5,26 +5,43 @@
 
 #define TINY 1.0e-20;
 
-
+/*!
+ * \brief LU decomposition
+ *
+ * \param a double **
+ * \param n int
+ * \param indx int *
+ * \param d double *
+ *
+ * \retrun 0 on singular matrix, 1 on success
+ */
 int G_ludcmp(double **a, int n, int *indx, double *d)
 {
     int i, imax = 0, j, k;
     double big, dum, sum, temp;
     double *vv;
+    int is_singular = FALSE;
 
     vv = G_alloc_vector(n);
     *d = 1.0;
+/* this pragma works, but doesn't really help speed things up */
+/* #pragma omp parallel for private(i, j, big, temp) shared(n, a, vv, is_singular) */
     for (i = 0; i < n; i++) {
 	big = 0.0;
 	for (j = 0; j < n; j++)
 	    if ((temp = fabs(a[i][j])) > big)
 		big = temp;
-	if (big == 0.0) {
-	    *d = 0.0;
-	    return 0;		/* Singular matrix  */
-	}
+
+	if (big == 0.0)
+	    is_singular = TRUE;
+
 	vv[i] = 1.0 / big;
     }
+    if (is_singular) {
+	*d = 0.0;
+	return 0;	/* Singular matrix  */
+    }
+
     for (j = 0; j < n; j++) {
 	for (i = 0; i < j; i++) {
 	    sum = a[i][j];
@@ -32,7 +49,10 @@
 		sum -= a[i][k] * a[k][j];
 	    a[i][j] = sum;
 	}
+
 	big = 0.0;
+/* not very efficent, but this pragma helps speed things up a bit */
+#pragma omp parallel for private(i, k, sum, dum) shared(j, n, a, vv, big, imax)
 	for (i = j; i < n; i++) {
 	    sum = a[i][j];
 	    for (k = 0; k < j; k++)
@@ -68,6 +88,17 @@
 
 #undef TINY
 
+
+/*!
+ * \brief LU backward substitution
+ *
+ * \param a double **
+ * \param n int
+ * \param indx int *
+ * \param b double []
+ *
+ * \retrun void
+ */
 void G_lubksb(double **a, int n, int *indx, double b[])
 {
     int i, ii, ip, j;