[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: solvetoep
From: |
Gorazd Brumen |
Subject: |
Re: solvetoep |
Date: |
Tue, 28 Nov 2006 14:55:29 +0100 |
User-agent: |
Thunderbird 1.5.0.8 (X11/20061116) |
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
- solvetoep, Gorazd Brumen, 2006/11/28
- Re: solvetoep, David Bateman, 2006/11/28
- Re: solvetoep,
Gorazd Brumen <=
- Re: solvetoep, David Bateman, 2006/11/28
- Re: solvetoep, David Bateman, 2006/11/28
- Re: solvetoep, Gorazd Brumen, 2006/11/28
- Re: solvetoep, David Bateman, 2006/11/28
- Re: solvetoep, David Bateman, 2006/11/28
- Re: solvetoep, Gorazd Brumen, 2006/11/28