Numerical Linear Algebra for Applications in Statistics. JE Gentle.
Springer. 1998. ISBN 0-387-98542-5.
- */
+*/
#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
- __ __
+ __ __
|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;
+ if (A == NULL)
+ return GSL_EFAULT;
+
+ if (A->size1 != A->size2)
+ return GSL_ENOTSQR;
+
int i;
int j;
int k;
- int row_i;
- int row_k;
- int row_j;
- int *ordered_cols;
gsl_matrix *B;
- if (A != NULL)
+
+ 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 - 1; 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++)
+ gsl_matrix_set (B, k, k, -1.0 / sweep_element);
+ /*
+ Rows before current row k.
+ */
+ for (i = 0; i < k; i++)
{
- 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;
+ gsl_matrix_set (B, i, j, tmp);
}
- /*
- 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;
+ gsl_matrix_set (B, i, j, tmp);
}
- /*
- 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);
}
}
- 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);
- return GSL_SUCCESS;
- free (ordered_cols);
+ /*
+ 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);
+ }
+ }
}
- return GSL_ENOTSQR;
+ for (i = 0; i < A->size1; i++)
+ for (j = i; j < A->size2; j++)
+ {
+ gsl_matrix_set (A, i, j, gsl_matrix_get (B, i, j));
+ }
}
- 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;
}
projects / pspp-builds.git / blobdiff