Remove duplicated functions
[pspp] / lib / linreg / sweep.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2005, 2009 Free Software Foundation, Inc.
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation, either version 3 of the License, or
7    (at your option) any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License for more details.
13
14    You should have received a copy of the GNU General Public License
15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
16
17 /*
18   Find the least-squares estimate of b for the linear model:
19
20   Y = Xb + Z
21
22   where Y is an n-by-1 column vector, X is an n-by-p matrix of
23   independent variables, b is a p-by-1 vector of regression coefficients,
24   and Z is an n-by-1 normally-distributed random vector with independent
25   identically distributed components with mean 0.
26
27   This estimate is found via the sweep operator, which is a modification
28   of Gauss-Jordan pivoting.
29
30
31   References:
32
33   Matrix Computations, third edition. GH Golub and CF Van Loan.
34   The Johns Hopkins University Press. 1996. ISBN 0-8018-5414-8.
35
36   Numerical Analysis for Statisticians. K Lange. Springer. 1999.
37   ISBN 0-387-94979-8.
38
39   Numerical Linear Algebra for Applications in Statistics. JE Gentle.
40   Springer. 1998. ISBN 0-387-98542-5.
41  */
42
43 #include <config.h>
44
45 #include "sweep.h"
46
47 /*
48   The matrix A will be overwritten. In ordinary uses of the sweep
49   operator, A will be the matrix
50
51    __       __
52   |X'X    X'Y|
53   |          |
54   |Y'X    Y'Y|
55    --        --
56
57    X refers to the design matrix and Y to the vector of dependent
58    observations. reg_sweep sweeps on the diagonal elements of
59    X'X.
60
61    The matrix A is assumed to be symmetric, so the sweep operation is
62    performed only for the upper triangle of A.
63
64    LAST_COL is considered to be the final column in the augmented matrix,
65    that is, the column to the right of the '=' sign of the system.
66  */
67
68 int
69 reg_sweep (gsl_matrix * A, int last_col)
70 {
71   double sweep_element;
72   double tmp;
73   int i;
74   int j;
75   int k;
76   int row_i;
77   int row_k;
78   int row_j;
79   int *ordered_cols;
80   gsl_matrix *B;
81
82   if (A != NULL)
83     {
84       if (A->size1 == A->size2)
85         {
86           ordered_cols = malloc (A->size1 * sizeof (*ordered_cols));
87           for (i = 0; i < last_col; i++)
88             {
89               ordered_cols[i] = i;
90             }
91           for (i = last_col + 1; i < A->size1; i++)
92             {
93               ordered_cols[i - 1] = i;
94             }
95           ordered_cols[A->size1 - 1] = last_col;
96           B = gsl_matrix_alloc (A->size1, A->size2);
97           for (k = 0; k < (A->size1 - 1); k++)
98             {
99               row_k = ordered_cols[k];
100               sweep_element = gsl_matrix_get (A, row_k, row_k);
101               if (fabs (sweep_element) > GSL_DBL_MIN)
102                 {
103                   tmp = -1.0 / sweep_element;
104                   gsl_matrix_set (B, row_k, row_k, tmp);
105                   /*
106                      Rows before current row k.
107                    */
108                   for (i = 0; i < k; i++)
109                     {
110                       row_i = ordered_cols[i];
111                       for (j = i; j < A->size2; j++)
112                         {
113                           row_j = ordered_cols[j];
114                           /*
115                              Use only the upper triangle of A.
116                            */
117                           if (row_j < row_k)
118                             {
119                               tmp = gsl_matrix_get (A, row_i, row_j) -
120                                 gsl_matrix_get (A, row_i, row_k)
121                                 * gsl_matrix_get (A, row_j, row_k) / sweep_element;
122                               gsl_matrix_set (B, row_i, row_j, tmp);
123                             }
124                           else if (row_j > row_k)
125                             {
126                               tmp = gsl_matrix_get (A, row_i, row_j) -
127                                 gsl_matrix_get (A, row_i, row_k)
128                                 * gsl_matrix_get (A, row_k, row_j) / sweep_element;
129                               gsl_matrix_set (B, row_i, row_j, tmp);
130                             }
131                           else
132                             {
133                               tmp = gsl_matrix_get (A, row_i, row_k) / sweep_element;
134                               gsl_matrix_set (B, row_i, row_j, tmp);
135                             }
136                         }
137                     }
138                   /*
139                      Current row k.
140                    */
141                   for (j = k + 1; j < A->size1; j++)
142                     {
143                       row_j = ordered_cols[j];
144                       tmp = gsl_matrix_get (A, row_k, row_j) / sweep_element;
145                       gsl_matrix_set (B, row_k, row_j, tmp);
146                     }
147                   /*
148                      Rows after the current row k.
149                    */
150                   for (i = k + 1; i < A->size1; i++)
151                     {
152                       row_i = ordered_cols[i];
153                       for (j = i; j < A->size2; j++)
154                         {
155                           row_j = ordered_cols[j];
156                           tmp = gsl_matrix_get (A, row_i, row_j) -
157                             gsl_matrix_get (A, row_k, row_i)
158                             * gsl_matrix_get (A, row_k, row_j) / sweep_element;
159                           gsl_matrix_set (B, row_i, row_j, tmp);
160                         }
161                     }
162                 }
163               for (i = 0; i < A->size1; i++)
164                 for (j = i; j < A->size2; j++)
165                   {
166                     row_i = ordered_cols[i];
167                     row_j = ordered_cols[j];
168                     gsl_matrix_set (A, row_i, row_j, gsl_matrix_get (B, row_i, row_j));
169                   }
170             }
171           gsl_matrix_free (B);
172           free (ordered_cols);
173           return GSL_SUCCESS;
174         }
175       return GSL_ENOTSQR;
176     }
177   return GSL_EFAULT;
178 }