Add argument specifying column containing dependent variable.

[pspp] / lib / linreg / sweep.c

diff --git a/lib/linreg/sweep.c b/lib/linreg/sweep.c
index 000a5e3e0f740ed6b26957ccb2ee5703649db1e6..7d7dfc917fb3d846bdb29d7cc34564a8ea06b1d5 100644 (file)
--- 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;
     }