sweep.c: swap rows/columns instead of using indirection for last_col

[pspp-builds.git] / lib / linreg / sweep.c

diff --git a/lib/linreg/sweep.c b/lib/linreg/sweep.c
index 000a5e3e0f740ed6b26957ccb2ee5703649db1e6..b8b57e89e599d301fae49474e10de0157326cd83 100644 (file)
--- a/lib/linreg/sweep.c
+++ b/lib/linreg/sweep.c
@@ -44,6 +44,7 @@
 
 #include "sweep.h"
 
+#include <assert.h>
 /*
   The matrix A will be overwritten. In ordinary uses of the sweep
   operator, A will be the matrix
@@ -60,10 +61,13 @@
 
    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;
@@ -76,6 +80,10 @@ reg_sweep (gsl_matrix * A)
     {
       if (A->size1 == A->size2)
        {
+          assert (last_col < A->size1);
+         gsl_matrix_swap_rows (A, A->size1 - 1, last_col);
+         gsl_matrix_swap_columns (A, A->size1 - 1 , last_col);
+         
          B = gsl_matrix_alloc (A->size1, A->size2);
          for (k = 0; k < (A->size1 - 1); k++)
            {
@@ -144,6 +152,10 @@ reg_sweep (gsl_matrix * A)
                  }
            }
          gsl_matrix_free (B);
+
+         gsl_matrix_swap_columns (A, A->size1 - 1 , last_col);
+         gsl_matrix_swap_rows (A, A->size1 - 1, last_col);
+
          return GSL_SUCCESS;
        }
       return GSL_ENOTSQR;