|
Packit |
67cb25 |
/* Cholesky Decomposition
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2000 Thomas Walter
|
|
Packit |
67cb25 |
* Copyright (C) 2000, 2001, 2002, 2003, 2005, 2007 Brian Gough, Gerard Jungman
|
|
Packit |
67cb25 |
* Copyright (C) 2016 Patrick Alken
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This is free software; you can redistribute it and/or modify it
|
|
Packit |
67cb25 |
* under the terms of the GNU General Public License as published by the
|
|
Packit |
67cb25 |
* Free Software Foundation; either version 3, or (at your option) any
|
|
Packit |
67cb25 |
* later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This source is distributed in the hope that it will be useful, but WITHOUT
|
|
Packit |
67cb25 |
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
Packit |
67cb25 |
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
Packit |
67cb25 |
* for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* 03 May 2000: Modified for GSL by Brian Gough
|
|
Packit |
67cb25 |
* 29 Jul 2005: Additions by Gerard Jungman
|
|
Packit |
67cb25 |
* 04 Mar 2016: Change Cholesky algorithm to gaxpy version by Patrick Alken
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* Cholesky decomposition of a symmetrix positive definite matrix.
|
|
Packit |
67cb25 |
* This is useful to solve the matrix arising in
|
|
Packit |
67cb25 |
* periodic cubic splines
|
|
Packit |
67cb25 |
* approximating splines
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This algorithm does:
|
|
Packit |
67cb25 |
* A = L * L'
|
|
Packit |
67cb25 |
* with
|
|
Packit |
67cb25 |
* L := lower left triangle matrix
|
|
Packit |
67cb25 |
* L' := the transposed form of L.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_linalg.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double cholesky_norm1(const gsl_matrix * LLT, gsl_vector * work);
|
|
Packit |
67cb25 |
static int cholesky_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
In GSL 2.2, we decided to modify the behavior of the Cholesky decomposition
|
|
Packit |
67cb25 |
to store the Cholesky factor in the lower triangle, and store the original
|
|
Packit |
67cb25 |
matrix in the upper triangle. Previous versions stored the Cholesky factor in
|
|
Packit |
67cb25 |
both places. The routine gsl_linalg_cholesky_decomp1 was added for the new
|
|
Packit |
67cb25 |
behavior, and gsl_linalg_cholesky_decomp is maintained for backward compatibility.
|
|
Packit |
67cb25 |
It will be removed in a future release.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_decomp (gsl_matrix * A)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_linalg_cholesky_decomp1(A);
|
|
Packit |
67cb25 |
if (status == GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_transpose_tricpy('L', 0, A, A);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_decomp1()
|
|
Packit |
67cb25 |
Perform Cholesky decomposition of a symmetric positive
|
|
Packit |
67cb25 |
definite matrix using lower triangle
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - (input) symmetric, positive definite matrix
|
|
Packit |
67cb25 |
(output) lower triangle contains Cholesky factor
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success/error
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes:
|
|
Packit |
67cb25 |
1) Based on algorithm 4.2.1 (Gaxpy Cholesky) of Golub and
|
|
Packit |
67cb25 |
Van Loan, Matrix Computations (4th ed).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
2) original matrix is saved in upper triangle on output
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_decomp1 (gsl_matrix * A)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = A->size1;
|
|
Packit |
67cb25 |
const size_t N = A->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (M != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* save original matrix in upper triangle for later rcond calculation */
|
|
Packit |
67cb25 |
gsl_matrix_transpose_tricpy('L', 0, A, A);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ajj;
|
|
Packit |
67cb25 |
gsl_vector_view v = gsl_matrix_subcolumn(A, j, j, N - j); /* A(j:n,j) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (j > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_view w = gsl_matrix_subrow(A, j, 0, j); /* A(j,1:j-1)^T */
|
|
Packit |
67cb25 |
gsl_matrix_view m = gsl_matrix_submatrix(A, j, 0, N - j, j); /* A(j:n,1:j-1) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dgemv(CblasNoTrans, -1.0, &m.matrix, &w.vector, 1.0, &v.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ajj = gsl_matrix_get(A, j, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (ajj <= 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("matrix is not positive definite", GSL_EDOM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ajj = sqrt(ajj);
|
|
Packit |
67cb25 |
gsl_vector_scale(&v.vector, 1.0 / ajj);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_solve (const gsl_matrix * LLT,
|
|
Packit |
67cb25 |
const gsl_vector * b,
|
|
Packit |
67cb25 |
gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (LLT->size1 != LLT->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LLT->size1 != b->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LLT->size2 != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* copy x <- b */
|
|
Packit |
67cb25 |
gsl_vector_memcpy (x, b);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_linalg_cholesky_svx(LLT, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_svx (const gsl_matrix * LLT,
|
|
Packit |
67cb25 |
gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (LLT->size1 != LLT->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LLT->size2 != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* solve for c using forward-substitution, L c = b */
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* perform back-substitution, L^T x = c */
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LLT, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_solve_mat (const gsl_matrix * LLT,
|
|
Packit |
67cb25 |
const gsl_matrix * B,
|
|
Packit |
67cb25 |
gsl_matrix * X)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (LLT->size1 != LLT->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LLT->size1 != B->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match B size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LLT->size2 != X->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* copy X <- B */
|
|
Packit |
67cb25 |
gsl_matrix_memcpy (X, B);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_linalg_cholesky_svx_mat(LLT, X);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_svx_mat (const gsl_matrix * LLT,
|
|
Packit |
67cb25 |
gsl_matrix * X)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (LLT->size1 != LLT->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LLT->size2 != X->size1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* solve for C using forward-substitution, L C = B */
|
|
Packit |
67cb25 |
gsl_blas_dtrsm (CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, 1.0,
|
|
Packit |
67cb25 |
LLT, X);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* perform back-substitution, L^T X = C */
|
|
Packit |
67cb25 |
gsl_blas_dtrsm (CblasLeft, CblasLower, CblasTrans, CblasNonUnit, 1.0,
|
|
Packit |
67cb25 |
LLT, X);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_invert()
|
|
Packit |
67cb25 |
Compute the inverse of a symmetric positive definite matrix in
|
|
Packit |
67cb25 |
Cholesky form.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: LLT - matrix in cholesky form on input
|
|
Packit |
67cb25 |
A^{-1} = L^{-t} L^{-1} on output
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success or error
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_invert(gsl_matrix * LLT)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (LLT->size1 != LLT->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = LLT->size1;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
gsl_vector_view v1, v2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* invert the lower triangle of LLT */
|
|
Packit |
67cb25 |
gsl_linalg_tri_lower_invert(LLT);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* The lower triangle of LLT now contains L^{-1}. Now compute
|
|
Packit |
67cb25 |
* A^{-1} = L^{-T} L^{-1}
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double aii = gsl_matrix_get(LLT, i, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i < N - 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double tmp;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
v1 = gsl_matrix_subcolumn(LLT, i, i, N - i);
|
|
Packit |
67cb25 |
gsl_blas_ddot(&v1.vector, &v1.vector, &tmp);
|
|
Packit |
67cb25 |
gsl_matrix_set(LLT, i, i, tmp);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_view m = gsl_matrix_submatrix(LLT, i + 1, 0, N - i - 1, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
v1 = gsl_matrix_subcolumn(LLT, i, i + 1, N - i - 1);
|
|
Packit |
67cb25 |
v2 = gsl_matrix_subrow(LLT, i, 0, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_blas_dgemv(CblasTrans, 1.0, &m.matrix, &v1.vector, aii, &v2.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
v1 = gsl_matrix_row(LLT, N - 1);
|
|
Packit |
67cb25 |
gsl_blas_dscal(aii, &v1.vector);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* copy lower triangle to upper */
|
|
Packit |
67cb25 |
gsl_matrix_transpose_tricpy('L', 0, LLT, LLT);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* gsl_linalg_cholesky_invert() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size1;
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* initial Cholesky */
|
|
Packit |
67cb25 |
int stat_chol = gsl_linalg_cholesky_decomp1(A);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(stat_chol == GSL_SUCCESS)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* calculate D from diagonal part of initial Cholesky */
|
|
Packit |
67cb25 |
for(i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double C_ii = gsl_matrix_get(A, i, i);
|
|
Packit |
67cb25 |
gsl_vector_set(D, i, C_ii*C_ii);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* multiply initial Cholesky by 1/sqrt(D) on the right */
|
|
Packit |
67cb25 |
for(i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for(j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix_set(A, i, j, gsl_matrix_get(A, i, j) / sqrt(gsl_vector_get(D, j)));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Because the initial Cholesky contained both L and transpose(L),
|
|
Packit |
67cb25 |
the result of the multiplication is not symmetric anymore;
|
|
Packit |
67cb25 |
but the lower triangle _is_ correct. Therefore we reflect
|
|
Packit |
67cb25 |
it to the upper triangle and declare victory.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
for(i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
for(j = i + 1; j < N; ++j)
|
|
Packit |
67cb25 |
gsl_matrix_set(A, i, j, gsl_matrix_get(A, j, i));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return stat_chol;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_scale()
|
|
Packit |
67cb25 |
This function computes scale factors diag(S), such that
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
diag(S) A diag(S)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
has a condition number within a factor N of the matrix
|
|
Packit |
67cb25 |
with the smallest condition number over all possible
|
|
Packit |
67cb25 |
diagonal scalings. See Corollary 7.6 of:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
N. J. Higham, Accuracy and Stability of Numerical Algorithms (2nd Edition),
|
|
Packit |
67cb25 |
SIAM, 2002.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - symmetric positive definite matrix
|
|
Packit |
67cb25 |
S - (output) scale factors, S_i = 1 / sqrt(A_ii)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_scale(const gsl_matrix * A, gsl_vector * S)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = A->size1;
|
|
Packit |
67cb25 |
const size_t N = A->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (M != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("A is not a square matrix", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (N != S->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("S must have length N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute S_i = 1/sqrt(A_{ii}) */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aii = gsl_matrix_get(A, i, i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (Aii <= 0.0)
|
|
Packit |
67cb25 |
gsl_vector_set(S, i, 1.0); /* matrix not positive definite */
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
gsl_vector_set(S, i, 1.0 / sqrt(Aii));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_scale_apply()
|
|
Packit |
67cb25 |
This function applies scale transformation to A:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A <- diag(S) A diag(S)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - (input/output)
|
|
Packit |
67cb25 |
on input, symmetric positive definite matrix
|
|
Packit |
67cb25 |
on output, diag(S) * A * diag(S) in lower triangle
|
|
Packit |
67cb25 |
S - (input) scale factors
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_scale_apply(gsl_matrix * A, const gsl_vector * S)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = A->size1;
|
|
Packit |
67cb25 |
const size_t N = A->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (M != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("A is not a square matrix", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (N != S->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("S must have length N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute: A <- diag(S) A diag(S) using lower triangle */
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double sj = gsl_vector_get(S, j);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = j; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double si = gsl_vector_get(S, i);
|
|
Packit |
67cb25 |
double *Aij = gsl_matrix_ptr(A, i, j);
|
|
Packit |
67cb25 |
*Aij *= si * sj;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_decomp2(gsl_matrix * A, gsl_vector * S)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = A->size1;
|
|
Packit |
67cb25 |
const size_t N = A->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (M != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (N != S->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR("S must have length N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute scaling factors to reduce cond(A) */
|
|
Packit |
67cb25 |
status = gsl_linalg_cholesky_scale(A, S);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* apply scaling factors */
|
|
Packit |
67cb25 |
status = gsl_linalg_cholesky_scale_apply(A, S);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute Cholesky decomposition of diag(S) A diag(S) */
|
|
Packit |
67cb25 |
status = gsl_linalg_cholesky_decomp1(A);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_svx2 (const gsl_matrix * LLT,
|
|
Packit |
67cb25 |
const gsl_vector * S,
|
|
Packit |
67cb25 |
gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (LLT->size1 != LLT->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LLT->size2 != S->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match S", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LLT->size2 != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* b~ = diag(S) b */
|
|
Packit |
67cb25 |
gsl_vector_mul(x, S);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Solve for c using forward-substitution, L c = b~ */
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Perform back-substitution, L^T x~ = c */
|
|
Packit |
67cb25 |
gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LLT, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute original solution vector x = S x~ */
|
|
Packit |
67cb25 |
gsl_vector_mul(x, S);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_solve2 (const gsl_matrix * LLT,
|
|
Packit |
67cb25 |
const gsl_vector * S,
|
|
Packit |
67cb25 |
const gsl_vector * b,
|
|
Packit |
67cb25 |
gsl_vector * x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (LLT->size1 != LLT->size2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LLT->size1 != S->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match S size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LLT->size1 != b->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (LLT->size2 != x->size)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Copy x <- b */
|
|
Packit |
67cb25 |
gsl_vector_memcpy (x, b);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = gsl_linalg_cholesky_svx2(LLT, S, x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_cholesky_rcond (const gsl_matrix * LLT, double * rcond,
|
|
Packit |
67cb25 |
gsl_vector * work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t M = LLT->size1;
|
|
Packit |
67cb25 |
const size_t N = LLT->size2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (M != N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (work->size != 3 * N)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("work vector must have length 3*N", GSL_EBADLEN);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
double Anorm = cholesky_norm1(LLT, work); /* ||A||_1 */
|
|
Packit |
67cb25 |
double Ainvnorm; /* ||A^{-1}||_1 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*rcond = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* don't continue if matrix is singular */
|
|
Packit |
67cb25 |
if (Anorm == 0.0)
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* estimate ||A^{-1}||_1 */
|
|
Packit |
67cb25 |
status = gsl_linalg_invnorm1(N, cholesky_Ainv, (void *) LLT, &Ainvnorm, work);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (Ainvnorm != 0.0)
|
|
Packit |
67cb25 |
*rcond = (1.0 / Anorm) / Ainvnorm;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute 1-norm of original matrix, stored in upper triangle of LLT;
|
|
Packit |
67cb25 |
* diagonal entries have to be reconstructed */
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
cholesky_norm1(const gsl_matrix * LLT, gsl_vector * work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = LLT->size1;
|
|
Packit |
67cb25 |
double max = 0.0;
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double sum = 0.0;
|
|
Packit |
67cb25 |
gsl_vector_const_view lj = gsl_matrix_const_subrow(LLT, j, 0, j + 1);
|
|
Packit |
67cb25 |
double Ajj;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute diagonal (j,j) entry of A */
|
|
Packit |
67cb25 |
gsl_blas_ddot(&lj.vector, &lj.vector, &Ajj);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < j; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double *wi = gsl_vector_ptr(work, i);
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get(LLT, i, j);
|
|
Packit |
67cb25 |
double absAij = fabs(Aij);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
sum += absAij;
|
|
Packit |
67cb25 |
*wi += absAij;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_set(work, j, sum + fabs(Ajj));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double wi = gsl_vector_get(work, i);
|
|
Packit |
67cb25 |
max = GSL_MAX(max, wi);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return max;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* x := A^{-1} x = A^{-t} x, A = L L^T */
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
cholesky_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
gsl_matrix * A = (gsl_matrix * ) params;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
(void) TransA; /* unused parameter warning */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute L^{-1} x */
|
|
Packit |
67cb25 |
status = gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasNonUnit, A, x);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute L^{-t} x */
|
|
Packit |
67cb25 |
status = gsl_blas_dtrsv(CblasLower, CblasTrans, CblasNonUnit, A, x);
|
|
Packit |
67cb25 |
if (status)
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|