Use c-ctype functions in lexer.
[pspp-builds.git] / lib / linreg / sweep.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2005 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 "sweep.h"
44
45 /*
46   The matrix A will be overwritten. In ordinary uses of the sweep
47   operator, A will be the matrix
48
49    __       __
50   |X'X    X'Y|
51   |          |
52   |Y'X    Y'Y|
53    --        --
54
55    X refers to the design matrix and Y to the vector of dependent
56    observations. reg_sweep sweeps on the diagonal elements of
57    X'X.
58
59    The matrix A is assumed to be symmetric, so the sweep operation is
60    performed only for the upper triangle of A.
61  */
62
63 int
64 reg_sweep (gsl_matrix * A)
65 {
66   double sweep_element;
67   double tmp;
68   int i;
69   int j;
70   int k;
71   gsl_matrix *B;
72
73   if (A != NULL)
74     {
75       if (A->size1 == A->size2)
76         {
77           B = gsl_matrix_alloc (A->size1, A->size2);
78           for (k = 0; k < (A->size1 - 1); k++)
79             {
80               sweep_element = gsl_matrix_get (A, k, k);
81               if (fabs (sweep_element) > GSL_DBL_MIN)
82                 {
83                   tmp = -1.0 / sweep_element;
84                   gsl_matrix_set (B, k, k, tmp);
85                   /*
86                      Rows before current row k.
87                    */
88                   for (i = 0; i < k; i++)
89                     {
90                       for (j = i; j < A->size2; j++)
91                         {
92                           /*
93                              Use only the upper triangle of A.
94                            */
95                           if (j < k)
96                             {
97                               tmp = gsl_matrix_get (A, i, j) -
98                                 gsl_matrix_get (A, i, k)
99                                 * gsl_matrix_get (A, j, k) / sweep_element;
100                               gsl_matrix_set (B, i, j, tmp);
101                             }
102                           else if (j > k)
103                             {
104                               tmp = gsl_matrix_get (A, i, j) -
105                                 gsl_matrix_get (A, i, k)
106                                 * gsl_matrix_get (A, k, j) / sweep_element;
107                               gsl_matrix_set (B, i, j, tmp);
108                             }
109                           else
110                             {
111                               tmp = gsl_matrix_get (A, i, k) / sweep_element;
112                               gsl_matrix_set (B, i, j, tmp);
113                             }
114                         }
115                     }
116                   /*
117                      Current row k.
118                    */
119                   for (j = k + 1; j < A->size1; j++)
120                     {
121                       tmp = gsl_matrix_get (A, k, j) / sweep_element;
122                       gsl_matrix_set (B, k, j, tmp);
123                     }
124                   /*
125                      Rows after the current row k.
126                    */
127                   for (i = k + 1; i < A->size1; i++)
128                     {
129                       for (j = i; j < A->size2; j++)
130                         {
131                           tmp = gsl_matrix_get (A, i, j) -
132                             gsl_matrix_get (A, k, i)
133                             * gsl_matrix_get (A, k, j) / sweep_element;
134                           gsl_matrix_set (B, i, j, tmp);
135                         }
136                     }
137                 }
138               for (i = 0; i < A->size1; i++)
139                 for (j = i; j < A->size2; j++)
140                   {
141                     gsl_matrix_set (A, i, j, gsl_matrix_get (B, i, j));
142                   }
143             }
144           gsl_matrix_free (B);
145           return GSL_SUCCESS;
146         }
147       return GSL_ENOTSQR;
148     }
149   return GSL_EFAULT;
150 }