/* 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
int
reg_sweep (gsl_matrix * A, int last_col)
{
- if (A == NULL)
- return GSL_EFAULT;
-
- if (A->size1 != A->size2)
- return GSL_ENOTSQR;
-
- double sweep_element;
- double tmp;
int i;
int j;
int k;
gsl_matrix *B;
+ 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);
B = gsl_matrix_alloc (A->size1, A->size2);
for (k = 0; k < (A->size1 - 1); k++)
{
- sweep_element = gsl_matrix_get (A, k, k);
+ const double sweep_element = gsl_matrix_get (A, k, k);
if (fabs (sweep_element) > GSL_DBL_MIN)
{
- tmp = -1.0 / sweep_element;
- gsl_matrix_set (B, k, k, tmp);
+ gsl_matrix_set (B, k, k, -1.0 / sweep_element);
/*
Rows before current row k.
*/
{
for (j = i; j < A->size2; j++)
{
- /*
- Use only the upper triangle of A.
- */
+ /* Use only the upper triangle of A. */
+ double tmp;
if (j < k)
{
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);
}
else if (j > k)
{
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);
}
else
{
tmp = gsl_matrix_get (A, i, k) / sweep_element;
- gsl_matrix_set (B, i, j, tmp);
}
+ gsl_matrix_set (B, i, j, tmp);
}
}
/*
*/
for (j = k + 1; j < A->size1; j++)
{
- tmp = gsl_matrix_get (A, k, j) / sweep_element;
+ double tmp = gsl_matrix_get (A, k, j) / sweep_element;
gsl_matrix_set (B, k, j, tmp);
}
/*
{
for (j = i; j < A->size2; j++)
{
- tmp = gsl_matrix_get (A, i, 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++)
- {
- gsl_matrix_set (A, i, j, gsl_matrix_get (B, i, j));
- }
+ gsl_matrix_memcpy (A, B);
}
gsl_matrix_free (B);