Please find attached the Cholesky invert for complex matrices and the testing
functions.
I used gcc to compile and test it.
Regards,
Huan
-----Original Message-----
From: Brian Gough [mailto:[email protected]]
Sent: Monday, May 24, 2010 1:56 PM
To: Huan Wu
Cc: Patrick Alken
Subject: Re: gsl_linalg_complex_cholesky_invert
At Fri, 14 May 2010 11:18:37 -0400,
Huan Wu wrote:
> I noticed that the Cholesky invert function for complex Hermitian matrix
> is not available yet from release 1.14. I'm working on a project that
> needs it so I did the extension of the available function for real
> symmetric matrix (gsl_linalg_cholesky_invert, added to GSL in release
> 1.12), to the one for complex Hermitian matrix. I'd like to contribute
> the new one (gsl_linalg_complex_cholesky_invert) to the GSL but not sure
> if anybody else has already planned to do that (forgive my ignorance to
> the GSL community if I missed any should-to-know information as a first
> time user).
> Thanks,
>
> Huan
>
> Huan Wu
>
> Research In Motion Limited
>
Hello
Thanks for your email. Sorry for the delay in replying, I was on
holiday last week. I'd like to include your function if you can also
make a test function, similar to test_cholesky_invert in test.c
(please send any patches to [email protected]). I am pretty sure no-one
else is working on this.
--
Brian Gough
---------------------------------------------------------------------
This transmission (including any attachments) may contain confidential
information, privileged material (including material protected by the
solicitor-client or other applicable privileges), or constitute non-public
information. Any use of this information by anyone other than the intended
recipient is prohibited. If you have received this transmission in error,
please immediately reply to the sender and delete this information from your
system. Use, dissemination, distribution, or reproduction of this transmission
by unintended recipients is not authorized and may be unlawful.
/******************************************************************************
gsl_linalg_complex_cholesky_invert()
Compute the inverse of an Hermitian positive definite matrix in
Cholesky form.
Inputs: LLT - matrix in cholesky form on input
A^{-1} = L^{-H} L^{-1} on output
Return: success or error
******************************************************************************/
int
gsl_linalg_complex_cholesky_invert(gsl_matrix_complex * LLT)
{
if (LLT->size1 != LLT->size2)
{
GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
}
else
{
size_t N = LLT->size1;
size_t i, j;
gsl_vector_complex_view v1;
/* invert the lower triangle of LLT */
for (i = 0; i < N; ++i)
{
double ajj;
gsl_complex z;
j = N - i - 1;
ajj = 1.0 / GSL_REAL(gsl_matrix_complex_get(LLT, j, j));
GSL_SET_COMPLEX(&z, ajj, 0.0);
gsl_matrix_complex_set(LLT, j, j, z);
ajj = -GSL_REAL(gsl_matrix_complex_get(LLT, j, j));
if (j < N - 1)
{
gsl_matrix_complex_view m;
m = gsl_matrix_complex_submatrix(LLT, j + 1, j + 1,
N - j - 1, N - j - 1);
v1 = gsl_matrix_complex_subcolumn(LLT, j, j + 1, N - j - 1);
gsl_blas_ztrmv(CblasLower, CblasNoTrans, CblasNonUnit,
&m.matrix, &v1.vector);
gsl_blas_zdscal(ajj, &v1.vector);
}
} /* for (i = 0; i < N; ++i) */
/*
* The lower triangle of LLT now contains L^{-1}. Now compute
* A^{-1} = L^{-H} L^{-1}
*
* The (ij) element of A^{-1} is column i of conj(L^{-1}) dotted into
* column j of L^{-1}
*/
for (i = 0; i < N; ++i)
{
gsl_complex sum;
for (j = i + 1; j < N; ++j)
{
gsl_vector_complex_view v2;
v1 = gsl_matrix_complex_subcolumn(LLT, i, j, N - j);
v2 = gsl_matrix_complex_subcolumn(LLT, j, j, N - j);
/* compute Ainv[i,j] = sum_k{conj(Linv[k,i]) * Linv[k,j]} */
gsl_blas_zdotc(&v1.vector, &v2.vector, &sum);
/* store in upper triangle */
gsl_matrix_complex_set(LLT, i, j, sum);
}
/* now compute the diagonal element */
v1 = gsl_matrix_complex_subcolumn(LLT, i, i, N - i);
gsl_blas_zdotc(&v1.vector, &v1.vector, &sum);
gsl_matrix_complex_set(LLT, i, i, sum);
}
/* copy the Hermitian upper triangle to the lower triangle */
for (j = 1; j < N; j++)
{
for (i = 0; i < j; i++)
{
gsl_complex z = gsl_matrix_complex_get(LLT, i, j);
gsl_matrix_complex_set(LLT, j, i, gsl_complex_conjugate(z));
}
}
return GSL_SUCCESS;
}
} /* gsl_linalg_complex_cholesky_invert() */
/*----------------
Testing Functions
----------------*/
int
test_choleskyc_invert_dim(const gsl_matrix_complex * m, double eps)
{
int s = 0;
unsigned long i, j, N = m->size1;
gsl_complex af, bt;
gsl_matrix_complex * v = gsl_matrix_complex_alloc(N, N);
gsl_matrix_complex * c = gsl_matrix_complex_alloc(N, N);
gsl_matrix_complex_memcpy(v, m);
s += gsl_linalg_complex_cholesky_decomp(v);
s += gsl_linalg_complex_cholesky_invert(v);
GSL_SET_COMPLEX(&af, 1.0, 0.0);
GSL_SET_COMPLEX(&bt, 0.0, 0.0);
gsl_blas_zhemm(CblasLeft, CblasUpper, af, m, v, bt, c);
/* c should be the identity matrix */
for (i = 0; i < N; ++i)
{
for (j = 0; j < N; ++j)
{
int foo;
gsl_complex cij = gsl_matrix_complex_get(c, i, j);
double expected, actual;
/* check real part */
if (i == j)
expected = 1.0;
else
expected = 0.0;
actual = GSL_REAL(cij);
foo = check(actual, expected, eps);
if (foo)
printf("REAL(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
N, N, i,j, actual, expected);
s += foo;
/* check imaginary part */
expected = 0.0;
actual = GSL_IMAG(cij);
foo = check(actual, expected, eps);
if (foo)
printf("IMAG(%3lu,%3lu)[%lu,%lu]: %22.18g %22.18g\n",
N, N, i,j, actual, expected);
s += foo;
}
}
gsl_matrix_complex_free(v);
gsl_matrix_complex_free(c);
return s;
}
int
test_choleskyc_invert(void)
{
int f;
int s = 0;
double dat2[] = {
92.303, 0.000, 10.858, 1.798,
10.858, -1.798, 89.027, 0.000
};
double dat3[] = {
59.75,0, 49.25,172.25, 66.75,-162.75,
49.25,-172.25, 555.5,0, -429,-333.5,
66.75,162.75, -429,333.5, 536.5,0
};
double dat4[] = {
102.108, 0.000, 14.721, 1.343, -17.480, 15.591, 3.308, -2.936,
14.721, -1.343, 101.970, 0.000, 11.671, -6.776, -5.009, -2.665,
-17.480, -15.591, 11.671, 6.776, 105.071, 0.000, 3.396, 6.276,
3.308, 2.936, -5.009, 2.665, 3.396, -6.276, 107.128, 0.000
};
gsl_matrix_complex_view rv2 = gsl_matrix_complex_view_array(dat2, 2, 2);
f = test_choleskyc_invert_dim(&rv2.matrix, 2 * 8.0 * GSL_DBL_EPSILON);
gsl_test(f, " choleskyc_invert 2x2 Hermitian");
s += f;
gsl_matrix_complex_view rv3 = gsl_matrix_complex_view_array(dat3, 3, 3);
f = test_choleskyc_invert_dim(&rv3.matrix, 2 * 1024.0 * GSL_DBL_EPSILON);
gsl_test(f, " choleskyc_invert 3x3 Hermitian");
s += f;
gsl_matrix_complex_view rv4 = gsl_matrix_complex_view_array(dat4, 4, 4);
f = test_choleskyc_invert_dim(&rv4.matrix, 2 * 64.0 * GSL_DBL_EPSILON);
gsl_test(f, " choleskyc_invert 4x4 Hermitian");
s += f;
return s;
}
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
