|
Packit |
67cb25 |
/* linalg/condest.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2016 Patrick Alken
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <string.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.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 |
/*
|
|
Packit |
67cb25 |
* This module contains routines for estimating the condition number
|
|
Packit |
67cb25 |
* of matrices in the 1-norm. The algorithm is based on the paper,
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* [1] N. J. Higham, "FORTRAN codes for estimating the one-norm of
|
|
Packit |
67cb25 |
* a real or complex matrix, with applications to condition estimation",
|
|
Packit |
67cb25 |
* ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double condest_tri_norm1(CBLAS_UPLO_t Uplo, const gsl_matrix * A);
|
|
Packit |
67cb25 |
static int condest_tri_rcond(CBLAS_UPLO_t Uplo, const gsl_matrix * A,
|
|
Packit |
67cb25 |
double * rcond, gsl_vector * work);
|
|
Packit |
67cb25 |
static int condest_same_sign(const gsl_vector * x, const gsl_vector * y);
|
|
Packit |
67cb25 |
static int condest_invtriu(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);
|
|
Packit |
67cb25 |
static int condest_invtril(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_tri_upper_rcond()
|
|
Packit |
67cb25 |
Estimate reciprocal condition number of upper triangular matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - upper triangular matrix, N-by-N
|
|
Packit |
67cb25 |
rcond - (output) reciprocal condition number estimate
|
|
Packit |
67cb25 |
work - workspace, length 3*N
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success/error
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_tri_upper_rcond(const gsl_matrix * A, double * rcond, gsl_vector * work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status = condest_tri_rcond(CblasUpper, A, rcond, work);
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_tri_lower_rcond()
|
|
Packit |
67cb25 |
Estimate reciprocal condition number of lower triangular matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - lower triangular matrix, N-by-N
|
|
Packit |
67cb25 |
rcond - (output) reciprocal condition number estimate
|
|
Packit |
67cb25 |
work - workspace, length 3*N
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success/error
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_tri_lower_rcond(const gsl_matrix * A, double * rcond, gsl_vector * work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status = condest_tri_rcond(CblasLower, A, rcond, work);
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_linalg_invnorm1()
|
|
Packit |
67cb25 |
Estimate the 1-norm of ||A^{-1}||, where A is a square
|
|
Packit |
67cb25 |
N-by-N matrix
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: N - size of matrix
|
|
Packit |
67cb25 |
Ainvx - pointer to function which calculates:
|
|
Packit |
67cb25 |
x := A^{-1} x or x := A^{-t} x
|
|
Packit |
67cb25 |
params - parameters to pass to Ainvx
|
|
Packit |
67cb25 |
Ainvnorm - (output) estimate of ||A^{-1}||_1
|
|
Packit |
67cb25 |
work - workspace, length 3*N
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_linalg_invnorm1(const size_t N,
|
|
Packit |
67cb25 |
int (* Ainvx)(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params),
|
|
Packit |
67cb25 |
void * params, double * Ainvnorm, gsl_vector * work)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
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 |
const size_t maxit = 5;
|
|
Packit |
67cb25 |
gsl_vector_view x = gsl_vector_subvector(work, 0, N);
|
|
Packit |
67cb25 |
gsl_vector_view v = gsl_vector_subvector(work, N, N);
|
|
Packit |
67cb25 |
gsl_vector_view xi = gsl_vector_subvector(work, 2*N, N);
|
|
Packit |
67cb25 |
double gamma, gamma_old, temp;
|
|
Packit |
67cb25 |
size_t i, k;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
gsl_vector_set(&x.vector, i, 1.0 / (double) N);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute v = A^{-1} x */
|
|
Packit |
67cb25 |
gsl_vector_memcpy(&v.vector, &x.vector);
|
|
Packit |
67cb25 |
(*Ainvx)(CblasNoTrans, &v.vector, params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* gamma = ||v||_1 */
|
|
Packit |
67cb25 |
gamma = gsl_blas_dasum(&v.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* xi = sign(v) */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double vi = gsl_vector_get(&v.vector, i);
|
|
Packit |
67cb25 |
gsl_vector_set(&xi.vector, i, GSL_SIGN(vi));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* x = A^{-t} xi */
|
|
Packit |
67cb25 |
gsl_vector_memcpy(&x.vector, &xi.vector);
|
|
Packit |
67cb25 |
(*Ainvx)(CblasTrans, &x.vector, params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (k = 0; k < maxit; ++k)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t j = (size_t) gsl_blas_idamax(&x.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* v := A^{-1} e_j */
|
|
Packit |
67cb25 |
gsl_vector_set_zero(&v.vector);
|
|
Packit |
67cb25 |
gsl_vector_set(&v.vector, j, 1.0);
|
|
Packit |
67cb25 |
(*Ainvx)(CblasNoTrans, &v.vector, params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gamma_old = gamma;
|
|
Packit |
67cb25 |
gamma = gsl_blas_dasum(&v.vector);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check for repeated sign vector (algorithm has converged) */
|
|
Packit |
67cb25 |
if (condest_same_sign(&v.vector, &xi.vector) || (gamma < gamma_old))
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* xi = sign(v) */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double vi = gsl_vector_get(&v.vector, i);
|
|
Packit |
67cb25 |
gsl_vector_set(&xi.vector, i, GSL_SIGN(vi));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* x = A^{-t} sign(v) */
|
|
Packit |
67cb25 |
gsl_vector_memcpy(&x.vector, &xi.vector);
|
|
Packit |
67cb25 |
(*Ainvx)(CblasTrans, &x.vector, params);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
temp = 1.0; /* (-1)^i */
|
|
Packit |
67cb25 |
for (i = 0; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double term = 1.0 + (double) i / (N - 1.0);
|
|
Packit |
67cb25 |
gsl_vector_set(&x.vector, i, temp * term);
|
|
Packit |
67cb25 |
temp = -temp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* x := A^{-1} x */
|
|
Packit |
67cb25 |
(*Ainvx)(CblasNoTrans, &x.vector, params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
temp = 2.0 * gsl_blas_dasum(&x.vector) / (3.0 * N);
|
|
Packit |
67cb25 |
if (temp > gamma)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_memcpy(&v.vector, &x.vector);
|
|
Packit |
67cb25 |
gamma = temp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*Ainvnorm = gamma;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
condest_tri_rcond(CBLAS_UPLO_t Uplo, const gsl_matrix * A, double * rcond, gsl_vector * work)
|
|
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 ("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 = condest_tri_norm1(Uplo, A); /* ||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 |
if (Uplo == CblasUpper)
|
|
Packit |
67cb25 |
status = gsl_linalg_invnorm1(N, condest_invtriu, (void *) A, &Ainvnorm, work);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
status = gsl_linalg_invnorm1(N, condest_invtril, (void *) A, &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 |
/* calculate 1 norm of triangular matrix */
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
condest_tri_norm1(CBLAS_UPLO_t Uplo, const gsl_matrix * A)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t N = A->size2;
|
|
Packit |
67cb25 |
double max = 0.0;
|
|
Packit |
67cb25 |
size_t i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (Uplo == CblasUpper)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double sum = 0.0;
|
|
Packit |
67cb25 |
for (i = 0; i <= j; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get(A, i, j);
|
|
Packit |
67cb25 |
sum += fabs(Aij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
max = GSL_MAX(max, sum);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
for (j = 0; j < N; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double sum = 0.0;
|
|
Packit |
67cb25 |
for (i = j; i < N; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double Aij = gsl_matrix_get(A, i, j);
|
|
Packit |
67cb25 |
sum += fabs(Aij);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
max = GSL_MAX(max, sum);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return max;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* return 1 if sign(x) = sign(y), 0 otherwise */
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
condest_same_sign(const gsl_vector * x, const gsl_vector * y)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const size_t n = x->size;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < n; ++i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double xi = gsl_vector_get(x, i);
|
|
Packit |
67cb25 |
double yi = gsl_vector_get(y, i);
|
|
Packit |
67cb25 |
if (GSL_SIGN(xi) != GSL_SIGN(yi))
|
|
Packit |
67cb25 |
return 0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* x := A^{-1} x, A upper triangular */
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
condest_invtriu(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * A = (gsl_matrix *) params;
|
|
Packit |
67cb25 |
return gsl_blas_dtrsv(CblasUpper, TransA, CblasNonUnit, A, x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* x := A^{-1} x, A lower triangular */
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
condest_invtril(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_matrix * A = (gsl_matrix *) params;
|
|
Packit |
67cb25 |
return gsl_blas_dtrsv(CblasLower, TransA, CblasNonUnit, A, x);
|
|
Packit |
67cb25 |
}
|