diff options
Diffstat (limited to 'gsl-1.9/multifit/covar.c')
| -rw-r--r-- | gsl-1.9/multifit/covar.c | 192 |
1 files changed, 192 insertions, 0 deletions
diff --git a/gsl-1.9/multifit/covar.c b/gsl-1.9/multifit/covar.c new file mode 100644 index 0000000..05683a7 --- /dev/null +++ b/gsl-1.9/multifit/covar.c @@ -0,0 +1,192 @@ +/* multifit/covar.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Brian Gough + * + * 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 + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_permutation.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_multifit_nlin.h> + +/* Compute the covariance matrix + + cov = inv (J^T J) + + by QRP^T decomposition of J +*/ + +int +gsl_multifit_covar (const gsl_matrix * J, double epsrel, gsl_matrix * covar) +{ + double tolr; + + size_t i, j, k; + size_t kmax = 0; + + gsl_matrix * r; + gsl_vector * tau; + gsl_vector * norm; + gsl_permutation * perm; + + size_t m = J->size1, n = J->size2 ; + + if (m < n) + { + GSL_ERROR ("Jacobian be rectangular M x N with M >= N", GSL_EBADLEN); + } + + if (covar->size1 != covar->size2 || covar->size1 != n) + { + GSL_ERROR ("covariance matrix must be square and match second dimension of jacobian", GSL_EBADLEN); + } + + r = gsl_matrix_alloc (m, n); + tau = gsl_vector_alloc (n); + perm = gsl_permutation_alloc (n) ; + norm = gsl_vector_alloc (n) ; + + { + int signum = 0; + gsl_matrix_memcpy (r, J); + gsl_linalg_QRPT_decomp (r, tau, perm, &signum, norm); + } + + + /* Form the inverse of R in the full upper triangle of R */ + + tolr = epsrel * fabs(gsl_matrix_get(r, 0, 0)); + + for (k = 0 ; k < n ; k++) + { + double rkk = gsl_matrix_get(r, k, k); + + if (fabs(rkk) <= tolr) + { + break; + } + + gsl_matrix_set(r, k, k, 1.0/rkk); + + for (j = 0; j < k ; j++) + { + double t = gsl_matrix_get(r, j, k) / rkk; + gsl_matrix_set (r, j, k, 0.0); + + for (i = 0; i <= j; i++) + { + double rik = gsl_matrix_get (r, i, k); + double rij = gsl_matrix_get (r, i, j); + + gsl_matrix_set (r, i, k, rik - t * rij); + } + } + kmax = k; + } + + /* Form the full upper triangle of the inverse of R^T R in the full + upper triangle of R */ + + for (k = 0; k <= kmax ; k++) + { + for (j = 0; j < k; j++) + { + double rjk = gsl_matrix_get (r, j, k); + + for (i = 0; i <= j ; i++) + { + double rij = gsl_matrix_get (r, i, j); + double rik = gsl_matrix_get (r, i, k); + + gsl_matrix_set (r, i, j, rij + rjk * rik); + } + } + + { + double t = gsl_matrix_get (r, k, k); + + for (i = 0; i <= k; i++) + { + double rik = gsl_matrix_get (r, i, k); + + gsl_matrix_set (r, i, k, t * rik); + }; + } + } + + /* Form the full lower triangle of the covariance matrix in the + strict lower triangle of R and in w */ + + for (j = 0 ; j < n ; j++) + { + size_t pj = gsl_permutation_get (perm, j); + + for (i = 0; i <= j; i++) + { + size_t pi = gsl_permutation_get (perm, i); + + double rij; + + if (j > kmax) + { + gsl_matrix_set (r, i, j, 0.0); + rij = 0.0 ; + } + else + { + rij = gsl_matrix_get (r, i, j); + } + + if (pi > pj) + { + gsl_matrix_set (r, pi, pj, rij); + } + else if (pi < pj) + { + gsl_matrix_set (r, pj, pi, rij); + } + + } + + { + double rjj = gsl_matrix_get (r, j, j); + gsl_matrix_set (covar, pj, pj, rjj); + } + } + + + /* symmetrize the covariance matrix */ + + for (j = 0 ; j < n ; j++) + { + for (i = 0; i < j ; i++) + { + double rji = gsl_matrix_get (r, j, i); + + gsl_matrix_set (covar, j, i, rji); + gsl_matrix_set (covar, i, j, rji); + } + } + + gsl_matrix_free (r); + gsl_permutation_free (perm); + gsl_vector_free (tau); + gsl_vector_free (norm); + + return GSL_SUCCESS; +} |