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