LU decomposition: backsubstitution available?

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

LU decomposition: backsubstitution available?

Thomas Hoffmann-3
It is a known procedure for the solution of a sequence of
systems of linear equations Ax=b with only varying rhs to
make once a LU decompostion of A (function lu() of octave)
and then backsubstitute for the several b's, e.g.
  x=backsubs(l,u,b)
As e.g. RLaB has such an mechanism, I thought of similar
functionality in Octave, but did not find any.

Can anybody shed some light on this?

Thomas.


Reply | Threaded
Open this post in threaded view
|

Re: LU decomposition: backsubstitution available?

Mario Storti-4
>>>>> On Wed, 25 Mar 1998 10:27:40 +0100 (MET),
>>>>>      Thomas Hoffmann <[hidden email]> said:

> It is a known procedure for the solution of a sequence of
> systems of linear equations Ax=b with only varying rhs to
> make once a LU decompostion of A (function lu() of octave)
> and then backsubstitute for the several b's, e.g.
>   x=backsubs(l,u,b)
> As e.g. RLaB has such an mechanism, I thought of similar
> functionality in Octave, but did not find any.

> Can anybody shed some light on this?

> Thomas.

I faced this same problem  before. It would be   nice to have  someone
writing a 'backsubs' routine. It seems to me that it should be written
in fortran in order to be efficient (do to the triangular structure of
l and u).

Meanwhile, I think  that you make a  significant save in computational
effort if compute the inverse of A, say

> invA=inv(A);
> x1=invA*b1;
> x2=invA*b2;
> x3=invA*b3;
> .
> .

since  computing the  inverse requires  O(N^3) ops. and  matrix-vector
multiplication requires only O(N^2). Of course, if you know all the
rhs's at once, then you can put them incolumns in a matrix B and then
make a:

X= A \ B;

and solve all the systems at once.

Mario

%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%%%%%%<>%
Mario Alberto Storti                           | Fax: (54)(42) 55.09.44 |
Centro Internacional de Metodos Computacionales| Tel: (54)(42) 55.91.75 |
     en Ingenieria - CIMEC                     |........................|  
INTEC, Guemes 3450 - 3000 Santa Fe, Argentina                           |
Reply: [hidden email], http://venus.unl.edu.ar/gtm-eng.html |


Reply | Threaded
Open this post in threaded view
|

Re: LU decomposition: backsubstitution available?

Michael Hanke-2
On Wed, 25 Mar 1998, Mario Storti wrote:

> >>>>> On Wed, 25 Mar 1998 10:27:40 +0100 (MET),
> >>>>>      Thomas Hoffmann <[hidden email]> said:
>
> > It is a known procedure for the solution of a sequence of
> > systems of linear equations Ax=b with only varying rhs to
> > make once a LU decompostion of A (function lu() of octave)
> > and then backsubstitute for the several b's, e.g.
> >   x=backsubs(l,u,b)
> > As e.g. RLaB has such an mechanism, I thought of similar
> > functionality in Octave, but did not find any.
>
> > Can anybody shed some light on this?
>
> > Thomas.
>
> I faced this same problem  before. It would be   nice to have  someone
> writing a 'backsubs' routine. It seems to me that it should be written
> in fortran in order to be efficient (do to the triangular structure of
> l and u).
>
> Meanwhile, I think  that you make a  significant save in computational
> effort if compute the inverse of A, say
>
> > invA=inv(A);
> > x1=invA*b1;
> > x2=invA*b2;
> > x3=invA*b3;
> > .
> > .
>
> since  computing the  inverse requires  O(N^3) ops. and  matrix-vector
> multiplication requires only O(N^2). Of course, if you know all the
> rhs's at once, then you can put them incolumns in a matrix B and then
> make a:
>
> X= A \ B;
>
> and solve all the systems at once.
>
In a matlab paper, I found a better solution as a side remark for
sparse matrices: The matrix division routine checks before trying any
LU decomposition if the matrix to be "inverted" is triangular (or a
permuted triangular matrix). If so,
the simple forward/backward substitution takes place. The amount for
checking this property for full matrices is O(n^2). Maybe it is a
better idea to carry a tag for every matrix??

Michael

+---------------------------------------------------------------+
|  Michael Hanke                Royal Institute of Technology   |
|                               NADA                            |
|                               S-10044 Stockholm               |
|                               Sweden                          |
+---------------------------------------------------------------+
|  Visiting address:            Lindstedtsvaegen 3              |
|  Phone:                       + (46) (8) 790 6278             |
|  Fax:                         + (46) (8) 790 0930             |
|  Email:                       [hidden email]               |
|                               [hidden email]       |
+---------------------------------------------------------------+


Reply | Threaded
Open this post in threaded view
|

Re: LU decomposition: backsubstitution available?

John W. Eaton-6
On 25-Mar-1998, Michael Hanke <[hidden email]> wrote:

| In a matlab paper, I found a better solution as a side remark for
| sparse matrices: The matrix division routine checks before trying any
| LU decomposition if the matrix to be "inverted" is triangular (or a
| permuted triangular matrix). If so,
| the simple forward/backward substitution takes place. The amount for
| checking this property for full matrices is O(n^2). Maybe it is a
| better idea to carry a tag for every matrix??

Ugh.  Not very OO.  :-)

A better solution, and one that I would much prefer, would be to
convert the representation to be a triangular matrix object.

I have some plans to do this, but I'm not sure exactly how it will
happen.  One possibility is to convert Octave to use Blitz++ for
handling arrays and matrices (see http://monet.uwaterloo.ca/blitz/ for
more information about Blitz++).  I can see lots of reasons to do use
Blitz++, including the ability to handle special matrix classes
(including multidimensional matrices) relatively easily.  But there
are also some significant problems.  Blitz++ relies heavily on C++
templates, and gcc (even 2.8.1 and egcs) still seems to have trouble
with some of the Blitz++ code.  Also, Blitz++ is still in relatively
early stages of development, and so may not really be ready for some
time.

jwe