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 6e114266c7ac0f64b5d2d89729a3fb7f568555ce..b8b57e89e599d301fae49474e10de0157326cd83 100644 (file)
--- a/lib/linreg/sweep.c
+++ b/lib/linreg/sweep.c
@@ -1,30 +1,25 @@
-/* lib/linreg/sweep.c
+/* PSPP - a program for statistical analysis.
+   Copyright (C) 2005, 2009 Free Software Foundation, Inc.
 
- Copyright (C) 2005 Free Software Foundation, Inc.
- Written by Jason H Stover.
+   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
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
 
- 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
- the Free Software Foundation; either version 2 of the License, or (at
- your option) any later version.
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
 
- This program is distributed in the hope that it will be useful, but
- WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program; if not, write to the Free Software
- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
- 02111-1307, USA.
- */
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>. */
 
 /*
   Find the least-squares estimate of b for the linear model:
 
   Y = Xb + Z
 
-  where Y is an n-by-1 column vector, X is an n-by-p matrix of 
+  where Y is an n-by-1 column vector, X is an n-by-p matrix of
   independent variables, b is a p-by-1 vector of regression coefficients,
   and Z is an n-by-1 normally-distributed random vector with independent
   identically distributed components with mean 0.
@@ -45,8 +40,11 @@
   Springer. 1998. ISBN 0-387-98542-5.
  */
 
-#include "pspp_linreg.h"
+#include <config.h>
+
+#include "sweep.h"
 
+#include <assert.h>
 /*
   The matrix A will be overwritten. In ordinary uses of the sweep
   operator, A will be the matrix
@@ -58,15 +56,18 @@
    --        --
 
    X refers to the design matrix and Y to the vector of dependent
-   observations. pspp_reg_sweep sweeps on the diagonal elements of 
+   observations. reg_sweep sweeps on the diagonal elements of
    X'X.
 
    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
-pspp_reg_sweep (gsl_matrix * A)
+reg_sweep (gsl_matrix * A, int last_col)
 {
   double sweep_element;
   double tmp;
@@ -79,6 +80,10 @@ pspp_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++)
            {
@@ -95,7 +100,7 @@ pspp_reg_sweep (gsl_matrix * A)
                      for (j = i; j < A->size2; j++)
                        {
                          /*
-                            Use only the upper triangle of A. 
+                            Use only the upper triangle of A.
                           */
                          if (j < k)
                            {
@@ -147,6 +152,10 @@ pspp_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;