X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=lib%2Flinreg%2Fsweep.c;h=62dda200701ceda43e4040e3128b1b8121b9e279;hb=f2d4cc6e7a4d5948a2c0cf70883347000a79a2b0;hp=b72d24a7edfebfba71201b28a34d13a374c08804;hpb=6be497eee03defa0af28afe9662cd1c0e6a7ab41;p=pspp diff --git a/lib/linreg/sweep.c b/lib/linreg/sweep.c index b72d24a7ed..62dda20070 100644 --- a/lib/linreg/sweep.c +++ b/lib/linreg/sweep.c @@ -1,5 +1,5 @@ /* PSPP - a program for statistical analysis. - Copyright (C) 2005, 2009 Free Software Foundation, Inc. + Copyright (C) 2005, 2009, 2011 Free Software Foundation, Inc. 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 @@ -38,141 +38,115 @@ Numerical Linear Algebra for Applications in Statistics. JE Gentle. Springer. 1998. ISBN 0-387-98542-5. - */ +*/ #include #include "sweep.h" +#include /* The matrix A will be overwritten. In ordinary uses of the sweep operator, A will be the matrix - __ __ + __ __ |X'X X'Y| | | |Y'X Y'Y| - -- -- + -- -- - X refers to the design matrix and Y to the vector of dependent - observations. reg_sweep sweeps on the diagonal elements of - X'X. + X refers to the design matrix and Y to the vector of dependent + 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. + 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. - */ + 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, 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 == NULL) + return GSL_EFAULT; + + if (A->size1 != A->size2) + return GSL_ENOTSQR; + + 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++) { - if (A->size1 == A->size2) + const double sweep_element = gsl_matrix_get (A, k, k); + if (fabs (sweep_element) > GSL_DBL_MIN) { - 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; i++) + gsl_matrix_set (B, k, k, -1.0 / sweep_element); + /* + Rows before current row k. + */ + for (i = 0; i < k; 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++) - { - row_k = ordered_cols[k]; - sweep_element = gsl_matrix_get (A, row_k, row_k); - if (fabs (sweep_element) > GSL_DBL_MIN) + for (j = i; j < A->size2; j++) { - tmp = -1.0 / sweep_element; - gsl_matrix_set (B, row_k, row_k, tmp); - /* - Rows before current row k. - */ - for (i = 0; i < k; i++) + /* Use only the upper triangle of A. */ + double tmp; + if (j < k) { - 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 (row_j < row_k) - { - 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 (row_j > row_k) - { - 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, row_i, row_k) / sweep_element; - gsl_matrix_set (B, row_i, row_j, tmp); - } - } + tmp = gsl_matrix_get (A, i, j) - + gsl_matrix_get (A, i, k) + * gsl_matrix_get (A, j, k) / sweep_element; } - /* - Current row k. - */ - for (j = k + 1; j < A->size1; j++) + else if (j > k) { - 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); + tmp = gsl_matrix_get (A, i, j) - + gsl_matrix_get (A, i, k) + * gsl_matrix_get (A, k, j) / sweep_element; } - /* - Rows after the current row k. - */ - for (i = k + 1; i < A->size1; i++) + else { - row_i = ordered_cols[i]; - for (j = i; j < A->size2; j++) - { - 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); - } + tmp = gsl_matrix_get (A, i, k) / sweep_element; } + gsl_matrix_set (B, i, j, tmp); + } + } + /* + Current row k. + */ + for (j = k + 1; j < A->size1; j++) + { + double tmp = gsl_matrix_get (A, k, j) / sweep_element; + gsl_matrix_set (B, k, j, tmp); + } + /* + Rows after the current row k. + */ + for (i = k + 1; i < A->size1; i++) + { + for (j = i; j < A->size2; j++) + { + double 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); } - for (i = 0; i < A->size1; i++) - for (j = i; j < A->size2; 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); - free (ordered_cols); - return GSL_SUCCESS; } - return GSL_ENOTSQR; + gsl_matrix_memcpy (A, B); } - return GSL_EFAULT; + 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; }