X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=lib%2Flinreg%2Fsweep.c;h=b72d24a7edfebfba71201b28a34d13a374c08804;hb=4ab3b7df3996d059d8efd8e9b8c5bdbd1cae84a1;hp=bbab8e2aad4adcdf9f9cbe41ec5d863a5084837b;hpb=43b1296aafe7582e7dbe6c2b6a8b478d7d9b0fcf;p=pspp-builds.git

diff --git a/lib/linreg/sweep.c b/lib/linreg/sweep.c
index bbab8e2a..b72d24a7 100644
--- a/lib/linreg/sweep.c
+++ b/lib/linreg/sweep.c
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2005 Free Software Foundation, Inc.
+   Copyright (C) 2005, 2009 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -40,6 +40,8 @@
   Springer. 1998. ISBN 0-387-98542-5.
  */
 
+#include <config.h>
+
 #include "sweep.h"
 
 /*
@@ -58,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; 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);
 			    }
 			}
 		    }
@@ -118,30 +140,36 @@ 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);
+	  free (ordered_cols);
 	  return GSL_SUCCESS;
 	}
       return GSL_ENOTSQR;