Re: solvetoep

[Gorazd Brumen]

Hi David, all,

OK. I am a bit inexperienced and have only recently realized the
following points:

a) The code is a C++ rewrittement of the netlib solver. Perhaps it
would be better to use netlib fortran code and then do a
C++ wrapper around it, since then we dont have to "worry" about the
netlib development more. If some better algorithm comes into the
netlib - voila, our code changes automagically.

b) The code is O(n^2) which is better than the general O(n^3),
it is from netlib (supposedly good, have no idea how good the netlib
code is) and uses a lot less space than the general inverse 2*n, instead
of the full n^2 through

toeplitz (c,r)\b

The matter of fact is that there exist even better routines, which
do the job in O(n log ^2 (n)) time, but with a big coefficient and
are better than the proposed solver for matrices of sizes > 256.
I have not yet seen that code, but have read the papers about it.

c) I dont quite know what you mean by sparse.

Please, suggest what is to be done.

Gorazd


David Bateman wrote:
> 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);
>>
>> }
>>   
> 
> 

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