This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: gsl_linalg_solve_symm_cyc_tridiag problem
- From: Brian Gough <bjg at network-theory dot co dot uk>
- To: David Necas (Yeti) <yeti at physics dot muni dot cz>
- Cc: gsl-discuss at sources dot redhat dot com
- Date: Sun, 14 Apr 2002 23:16:58 +0100 (BST)
- Subject: Re: gsl_linalg_solve_symm_cyc_tridiag problem
- References: <20020406141543.B7885@physics.muni.cz>
David Necas (Yeti) writes:
> the gsl_linalg_solve_symm_cyc_tridiag() solver
> sometimes give quite a strange results.
I think the following patch should fix the problem with the symmetric
cyclic tridiagonal solver. There were some typographical errors in
the textbook that was referenced.
Brian
Index: tridiag.c
===================================================================
RCS file: /cvs/gsl/gsl/linalg/tridiag.c,v
retrieving revision 1.10
diff -c -r1.10 tridiag.c
*** tridiag.c 19 Nov 2001 22:32:06 -0000 1.10
--- tridiag.c 14 Apr 2002 22:16:11 -0000
***************
*** 164,172 ****
}
alpha[N - 2] = diag[d_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * gamma[N - 3];
gamma[N - 2] = (offdiag[o_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * delta[N - 3]) / alpha[N - 2];
! alpha[N - 1] = diag[d_stride * (N - 1)] - sum - offdiag[o_stride * (N - 2)] * gamma[N - 2] * gamma[N - 2];
/* update */
z[0] = b[0];
for (i = 1; i < N - 1; i++)
{
--- 164,173 ----
}
alpha[N - 2] = diag[d_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * gamma[N - 3];
gamma[N - 2] = (offdiag[o_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * delta[N - 3]) / alpha[N - 2];
! alpha[N - 1] = diag[d_stride * (N - 1)] - sum - alpha[(N - 2)] * gamma[N - 2] * gamma[N - 2];
/* update */
+
z[0] = b[0];
for (i = 1; i < N - 1; i++)
{
***************
*** 178,183 ****
--- 179,185 ----
sum += delta[i] * z[i];
}
z[N - 1] = b[b_stride * (N - 1)] - sum - gamma[N - 2] * z[N - 2];
+
for (i = 0; i < N; i++)
{
c[i] = z[i] / alpha[i];