solvetoep

[Gorazd Brumen]

Hi all,

I have rewritten a C++ implementation of the tsld fortran routine
to solve the toeplitz matrices, which I attach. It needs a lot
of twinking still, but it is kind of OK. I even included a GPL
preambule at the beginning.

Is anybody interested?

Regards,
Gorazd



-- 
Gorazd Brumen
Mail: address@hidden
WWW: http://www.gorazdbrumen.net/
PGP: Key at http://pgp.mit.edu, ID BCC93240

/* 
   Copyright (C) 2006 Gorazd Brumen
   
   This file is part of Octave.
   
   Octave is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.
   
   Octave is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
   FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
   for more details.
   
   You should have received a copy of the GNU General Public License
   along with Octave; see the file COPYING.  If not, write to the Free
   Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
   02110-1301, USA.
   

*/
#include <octave/oct.h>
/*
    Parameters:
  
    Input, integer N, the order of the matrix.
    
    Input, real A(2*N-1), the first row of the Toeplitz matrix, followed by
    the first column of the Toeplitz matrix beginning with the second element.
    
    Input, real B(N) the right hand side vector.
    
    Output, real X(N), the solution vector.  X and B may share the
    same storage.
  
*/
DEFUN_DLD (solvetoep, args, , 
           "- Loadable Function: x = solvetoep (n, a, b) \n \
solves the toeplitz where a is the first row, followed by \n \
the first column, b is the rhs of the system, n is the \n \
system size.")
{

  int result;
  int n = args(0).int_value();
  RowVector x(n);
  RowVector a( args(1).vector_value() );
  RowVector b( args(2).vector_value() );

  RowVector c1 (n-1);
  RowVector c2 (n-1);
  int i, nsub;
  double r1, r2, r3, r5, r6;


  if ( n < 1) 
    return octave_value(0);

  r1 = a(0);
  x(0) = b(0)/r1;
  
  if ( n == 1 ) 
    return octave_value (x);

  for (nsub =2; nsub <= n; nsub = nsub + 1) {

    r5 = a(n+nsub-2);
    r6 = a(nsub-1);

    if (nsub > 2) {

      c1(nsub-2) = r2;
      for (i=1; i<= nsub-2; i= i+1) {
        r5 = r5 + a(n+i-1) * c1(nsub-i-1);
        r6 = r6 + a(i) * c2(i-1);
      }
    }

    r2 = -r5/r1;
    r3 = -r6/r1;
    r1 = r1 + r5 * r3;
    
    if (nsub > 2) {

      r6 = c2(0);
      c2(nsub-2) = 0;

      for (i=2; i<=nsub-1; i= i+1) {
        r5 = c2(i-1);
        c2(i-1) = c1(i-1) * r3 + r6;
        c1(i-1) = c1(i-1) + r6 * r2;
        r6 = r5;
      }
    }

    c2(0) = r3;
    r5 = 0;
    for (i=1; i<= nsub-1; i= i+1) 
      r5 = r5 + a(n+i-1) * x(nsub-i-1);

    r6 = ( b(nsub-1) - r5 )/r1;
    
    for (i=1; i<=nsub-1; i=i+1) 
      x(i-1) = x(i-1) + c2(i-1) * r6;
    
    x(nsub-1) = r6;
  }

  return octave_value (x);

}