gsl_linalg_solve_symm_cyc_tridiag problem

[David Necas (Yeti) <yeti at physics dot muni dot cz>]

Hello,

the gsl_linalg_solve_symm_cyc_tridiag() solver
sometimes give quite a strange results. It gives
as a solution a vector satisfying all the
equations except the last one.  Below, I send a
demo program -- it solves the same cyclic
tridiagonal system using gauss elimination and the
specialized method and compares the results (for
reference, gsl_linalg_solve_symm_cyc_tridiag()
gives me as solution: -324, 1011, -1722, 714,
315).

Yeti


#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_permutation.h>

/* doesn't seem as an almost-singular case, but triggers the problem quite well
 */
static const double diag[] = { 4, 2, 1, 2, 4 };
static const double offdiag[] = { 1, 1, 1, 1, 1 };
static const double rhs[] = { 30, -24, 3, 21, -30 };

static const int N = sizeof(diag)/sizeof(double);

int
main(int argc, char *argv[])
{
  int i, j, sgn;
  gsl_vector_const_view diag_v = gsl_vector_const_view_array(diag, N);
  gsl_vector_const_view offdiag_v = gsl_vector_const_view_array(offdiag, N);
  gsl_vector_const_view rhs_v = gsl_vector_const_view_array(rhs, N);
  gsl_vector *sol_td_v;
  gsl_vector *sol_ge_v;
  gsl_vector *residual;
  gsl_matrix *symtd_m;
  gsl_matrix *a;
  gsl_permutation *p;

  sol_td_v = gsl_vector_alloc(N);
  sol_ge_v = gsl_vector_alloc(N);
  residual = gsl_vector_alloc(N);
  symtd_m = gsl_matrix_alloc(N, N);
  a = gsl_matrix_alloc(N, N);
  p = gsl_permutation_alloc(N);
  gsl_matrix_set_zero(symtd_m);
  for (i = 0; i < N; i++) gsl_matrix_set(symtd_m, i, i, diag[i]);
  for (i = 0; i < N; i++) gsl_matrix_set(symtd_m, (i+1)%N, i, offdiag[i]);
  for (i = 0; i < N; i++) gsl_matrix_set(symtd_m, i, (i+1)%N, offdiag[i]);
  gsl_matrix_memcpy(a, symtd_m);

  puts("matrix");
  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      printf("% g\t", gsl_matrix_get(symtd_m, i, j));
    }
    puts("");
  }

  gsl_linalg_solve_symm_cyc_tridiag(&diag_v.vector,
                                    &offdiag_v.vector,
                                    &rhs_v.vector,
                                    sol_td_v);
  gsl_linalg_LU_decomp(a, p, &sgn);
  gsl_linalg_LU_solve(a, p, &rhs_v.vector, sol_ge_v);
/*  gsl_linalg_LU_refine(symtd_m, a, p, &rhs_v.vector, sol_ge_v, residual); */
  puts("\nsymm_cyc_tridiag solution compared with LU decomposition");
  puts("i\t sol_ge\t sol_td");
  for (i = 0; i < N; i++) {
    printf("%d:\t% g\t% g\n", i,
           gsl_vector_get(sol_ge_v, i),
           gsl_vector_get(sol_td_v, i));
  }
  puts("\nrhs compared with lhs computed from each solution");
  puts("i\t rhs\t lhs_ge\t lhs_td");
  for (i = 0; i < N; i++) {
    double lhs_ge = gsl_vector_get(sol_ge_v, (i+N-1)%N) * offdiag[(i+N-1)%N]
                    + gsl_vector_get(sol_ge_v, i) * diag[i]
                    + gsl_vector_get(sol_ge_v, (i+1)%N) * offdiag[i];
    double lhs_td = gsl_vector_get(sol_td_v, (i+N-1)%N) * offdiag[(i+N-1)%N]
                    + gsl_vector_get(sol_td_v, i) * diag[i]
                    + gsl_vector_get(sol_td_v, (i+1)%N) * offdiag[i];
    printf("%d:\t% g\t% g\t% g\n", i, rhs[i], lhs_ge, lhs_td);
  }

  return 0;
}