Re: solvetoep

[David Bateman]

One thing I forgot, but for this to be added to the \ solver, then the
code to probe for toeplitz matrices will need to be very efficient. We
don't want to pay a cost of probing for a toeplitz matrix for general
matrix problems...

D.


Gorazd Brumen wrote:
> 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
>
>
>
>   
> ------------------------------------------------------------------------
>
> /* 
>    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);
>
> }
>   


-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary