From: Jason H Stover <jhs@franz.gcsu.edu>
Date: Tue, 20 Jul 2010 22:32:35 +0000 (-0400)
Subject: Add argument specifying column containing dependent variable.
X-Git-Tag: v0.7.6~331
X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f338eb59ae54d23708eaf2ffa009239807794653;p=pspp-builds.git

Add argument specifying column containing dependent variable.
---

diff --git a/lib/linreg/sweep.c b/lib/linreg/sweep.c
index 000a5e3e..7d7dfc91 100644
--- a/lib/linreg/sweep.c
+++ b/lib/linreg/sweep.c
@@ -60,58 +60,78 @@
 
    The matrix A is assumed to be symmetric, so the sweep operation is
    performed only for the upper triangle of A.
+
+   LAST_COL is considered to be the final column in the augmented matrix,
+   that is, the column to the right of the '=' sign of the system.
  */
 
 int
-reg_sweep (gsl_matrix * A)
+reg_sweep (gsl_matrix * A, int last_col)
 {
   double sweep_element;
   double tmp;
   int i;
   int j;
   int k;
+  int row_i;
+  int row_k;
+  int row_j;
+  int *ordered_cols;
   gsl_matrix *B;
 
   if (A != NULL)
     {
       if (A->size1 == A->size2)
 	{
+	  ordered_cols = malloc (A->size1 * sizeof (*ordered_cols));
+	  for (i = 0; i < last_col; i++)
+	    {
+	      ordered_cols[i] = i;
+	    }
+	  for (i = last_col + 1; i < A->size1 - 1; i++)
+	    {
+	      ordered_cols[i - 1] = i;
+	    }
+	  ordered_cols[A->size1 - 1] = last_col;
 	  B = gsl_matrix_alloc (A->size1, A->size2);
 	  for (k = 0; k < (A->size1 - 1); k++)
 	    {
-	      sweep_element = gsl_matrix_get (A, k, k);
+	      row_k = ordered_cols[k];
+	      sweep_element = gsl_matrix_get (A, row_k, row_k);
 	      if (fabs (sweep_element) > GSL_DBL_MIN)
 		{
 		  tmp = -1.0 / sweep_element;
-		  gsl_matrix_set (B, k, k, tmp);
+		  gsl_matrix_set (B, row_k, row_k, tmp);
 		  /*
 		     Rows before current row k.
 		   */
 		  for (i = 0; i < k; i++)
 		    {
+		      row_i = ordered_cols[i];
 		      for (j = i; j < A->size2; j++)
 			{
+			  row_j = ordered_cols[j];
 			  /*
 			     Use only the upper triangle of A.
 			   */
-			  if (j < k)
+			  if (row_j < row_k)
 			    {
-			      tmp = gsl_matrix_get (A, i, j) -
-				gsl_matrix_get (A, i, k)
-				* gsl_matrix_get (A, j, k) / sweep_element;
-			      gsl_matrix_set (B, i, j, tmp);
+			      tmp = gsl_matrix_get (A, row_i, row_j) -
+				gsl_matrix_get (A, row_i, row_k)
+				* gsl_matrix_get (A, row_j, row_k) / sweep_element;
+			      gsl_matrix_set (B, row_i, row_j, tmp);
 			    }
-			  else if (j > k)
+			  else if (row_j > row_k)
 			    {
-			      tmp = gsl_matrix_get (A, i, j) -
-				gsl_matrix_get (A, i, k)
-				* gsl_matrix_get (A, k, j) / sweep_element;
-			      gsl_matrix_set (B, i, j, tmp);
+			      tmp = gsl_matrix_get (A, row_i, row_j) -
+				gsl_matrix_get (A, row_i, row_k)
+				* gsl_matrix_get (A, row_k, row_j) / sweep_element;
+			      gsl_matrix_set (B, row_i, row_j, tmp);
 			    }
 			  else
 			    {
-			      tmp = gsl_matrix_get (A, i, k) / sweep_element;
-			      gsl_matrix_set (B, i, j, tmp);
+			      tmp = gsl_matrix_get (A, row_i, row_k) / sweep_element;
+			      gsl_matrix_set (B, row_i, row_j, tmp);
 			    }
 			}
 		    }
@@ -120,31 +140,37 @@ reg_sweep (gsl_matrix * A)
 		   */
 		  for (j = k + 1; j < A->size1; j++)
 		    {
-		      tmp = gsl_matrix_get (A, k, j) / sweep_element;
-		      gsl_matrix_set (B, k, j, tmp);
+		      row_j = ordered_cols[j];
+		      tmp = gsl_matrix_get (A, row_k, row_j) / sweep_element;
+		      gsl_matrix_set (B, row_k, row_j, tmp);
 		    }
 		  /*
 		     Rows after the current row k.
 		   */
 		  for (i = k + 1; i < A->size1; i++)
 		    {
+		      row_i = ordered_cols[i];
 		      for (j = i; j < A->size2; j++)
 			{
-			  tmp = gsl_matrix_get (A, i, j) -
-			    gsl_matrix_get (A, k, i)
-			    * gsl_matrix_get (A, k, j) / sweep_element;
-			  gsl_matrix_set (B, i, j, tmp);
+			  row_j = ordered_cols[j];
+			  tmp = gsl_matrix_get (A, row_i, row_j) -
+			    gsl_matrix_get (A, row_k, row_i)
+			    * gsl_matrix_get (A, row_k, row_j) / sweep_element;
+			  gsl_matrix_set (B, row_i, row_j, tmp);
 			}
 		    }
 		}
 	      for (i = 0; i < A->size1; i++)
 		for (j = i; j < A->size2; j++)
 		  {
-		    gsl_matrix_set (A, i, j, gsl_matrix_get (B, i, j));
+		    row_i = ordered_cols[i];
+		    row_j = ordered_cols[j];
+		    gsl_matrix_set (A, row_i, row_j, gsl_matrix_get (B, row_i, row_j));
 		  }
 	    }
 	  gsl_matrix_free (B);
 	  return GSL_SUCCESS;
+	  free (ordered_cols);
 	}
       return GSL_ENOTSQR;
     }
diff --git a/lib/linreg/sweep.h b/lib/linreg/sweep.h
index 7bdd5c7c..a1567a16 100644
--- a/lib/linreg/sweep.h
+++ b/lib/linreg/sweep.h
@@ -65,6 +65,6 @@
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_math.h>
 
-int reg_sweep (gsl_matrix * A);
+int reg_sweep (gsl_matrix *, int);
 
 #endif
diff --git a/src/math/linreg.c b/src/math/linreg.c
index e74646b9..fe01ac39 100644
--- a/src/math/linreg.c
+++ b/src/math/linreg.c
@@ -74,6 +74,7 @@ linreg_alloc (const struct variable *depvar, const struct variable **indep_vars,
   c = xmalloc (sizeof (*c));
   c->depvar = depvar;
   c->indep_vars = xnmalloc (p, sizeof (*indep_vars));
+  c->dependent_column = p;
   for (i = 0; i < p; i++)
     {
       c->indep_vars[i] = indep_vars[i];
@@ -375,7 +376,7 @@ linreg_fit (const gsl_matrix *cov, linreg *l)
       gsl_matrix *params;
       params = gsl_matrix_calloc (cov->size1, cov->size2);
       gsl_matrix_memcpy (params, cov);
-      reg_sweep (params);
+      reg_sweep (params, l->dependent_column);
       post_sweep_computations (l, params);  
       gsl_matrix_free (params);
     }
@@ -460,3 +461,9 @@ linreg_get_depvar_mean (linreg *c)
 {
   return c->depvar_mean;
 }
+static
+linreg_set_dependent_column (linreg *c, int column)
+{
+  c->dependent_column = column;
+}
+
diff --git a/src/math/linreg.h b/src/math/linreg.h
index 349d5a90..31bf5651 100644
--- a/src/math/linreg.h
+++ b/src/math/linreg.h
@@ -142,6 +142,7 @@ struct linreg_struct
 
   struct variable *pred;
   struct variable *resid;
+  int dependent_column; /* Column containing the dependent variable. Defaults to last column. */
 };
 
 typedef struct linreg_struct linreg;