help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Banded matrix with dirty off-diagonal triagles


From: Mario Storti
Subject: Re: Banded matrix with dirty off-diagonal triagles
Date: Wed, 17 Jun 1998 13:27:41 -0300

>>>>> On Wed, 17 Jun 1998 19:05:12 +0300,
>>>>>      Michael Smolsky <address@hidden> said:

> Dear Octave users,
>       I know, support of sparse matrices in on the waiting list. However, I'm
> interested in a particular case, which I would like to implement as an
> .oct file, as a temporary solution. Could anyone point me to a free
> fortran source of the code I need?

>       I need to solve a large sparse system with a banded matrix (i.e. a
> matrix, which is zero everywhere except for a few diagonals). Contrary
> to the well-examined case, in my problem my matrix also has "dirty"
> (non-zero) triangles in the upper-right and lower-left corners of the
> matrix. Such a system arizes in the case when a PDE equation has
> periodic boundary conditions (the "last" variable is related to the
> "first" exactly as the nth one to the n+1st).

>       I scaned gams.nist.gov, but didn't find a standard routine for
> inversion of such a matrix, or solution of a system of linear equations
> built it.

>       Thank you in advance for your help, MSi.

You can convert  this to a standard banded  problem by renumbering the
unknowns in  a  (rather tricky) way.   However, you got  a system with
twice the original bandwidth. Take for instance the  matrix for the 1d
Laplace operator on 20 nodes with cyclic b.c.'s.

octave> A=2*eye(20)-diag(ones(19,1),1)-diag(ones(19,1),-1);
octave> A
A =

   2  -1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  -1   2  -1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
   0  -1   2  -1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
   0   0  -1   2  -1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
   0   0   0  -1   2  -1   0   0   0   0   0   0   0   0   0   0   0   0   0   0
   0   0   0   0  -1   2  -1   0   0   0   0   0   0   0   0   0   0   0   0   0
   0   0   0   0   0  -1   2  -1   0   0   0   0   0   0   0   0   0   0   0   0
   0   0   0   0   0   0  -1   2  -1   0   0   0   0   0   0   0   0   0   0   0
   0   0   0   0   0   0   0  -1   2  -1   0   0   0   0   0   0   0   0   0   0
   0   0   0   0   0   0   0   0  -1   2  -1   0   0   0   0   0   0   0   0   0
   0   0   0   0   0   0   0   0   0  -1   2  -1   0   0   0   0   0   0   0   0
   0   0   0   0   0   0   0   0   0   0  -1   2  -1   0   0   0   0   0   0   0
   0   0   0   0   0   0   0   0   0   0   0  -1   2  -1   0   0   0   0   0   0
   0   0   0   0   0   0   0   0   0   0   0   0  -1   2  -1   0   0   0   0   0
   0   0   0   0   0   0   0   0   0   0   0   0   0  -1   2  -1   0   0   0   0
   0   0   0   0   0   0   0   0   0   0   0   0   0   0  -1   2  -1   0   0   0
   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  -1   2  -1   0   0
   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  -1   2  -1   0
   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  -1   2  -1
   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  -1   2


If you take the following renumbering:

octave> new=vec([(1:10); (20:-1:11)])
new =

   1
  20
   2
  19
   3
  18
   4
  17
   5
  16
   6
  15
   7
  14
   8
  13
   9
  12
  10
  11

then you got as new matrix:

octave> AA=A(new,:);
octave> AA=AA(:,new);
octave> AA
AA =

   2  -1  -1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  -1   2   0  -1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  -1   0   2   0  -1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
   0  -1   0   2   0  -1   0   0   0   0   0   0   0   0   0   0   0   0   0   0
   0   0  -1   0   2   0  -1   0   0   0   0   0   0   0   0   0   0   0   0   0
   0   0   0  -1   0   2   0  -1   0   0   0   0   0   0   0   0   0   0   0   0
   0   0   0   0  -1   0   2   0  -1   0   0   0   0   0   0   0   0   0   0   0
   0   0   0   0   0  -1   0   2   0  -1   0   0   0   0   0   0   0   0   0   0
   0   0   0   0   0   0  -1   0   2   0  -1   0   0   0   0   0   0   0   0   0
   0   0   0   0   0   0   0  -1   0   2   0  -1   0   0   0   0   0   0   0   0
   0   0   0   0   0   0   0   0  -1   0   2   0  -1   0   0   0   0   0   0   0
   0   0   0   0   0   0   0   0   0  -1   0   2   0  -1   0   0   0   0   0   0
   0   0   0   0   0   0   0   0   0   0  -1   0   2   0  -1   0   0   0   0   0
   0   0   0   0   0   0   0   0   0   0   0  -1   0   2   0  -1   0   0   0   0
   0   0   0   0   0   0   0   0   0   0   0   0  -1   0   2   0  -1   0   0   0
   0   0   0   0   0   0   0   0   0   0   0   0   0  -1   0   2   0  -1   0   0
   0   0   0   0   0   0   0   0   0   0   0   0   0   0  -1   0   2   0  -1   0
   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  -1   0   2   0  -1
   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  -1   0   2  -1
   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  -1  -1   2

octave>


This  matrix is  banded but with   bandwidth m=2 in  contrast with the
original one, which had m=1.  I think that this  is optimal. It is not
magic, it is  the result of applying the  Cuthill Mc-Kee algorithm. On
the graph.  The graph of the matrix is something like this:

     
1---2---3---4---5---6---7---8---9---10---11---12---13---14---15---16---17---18---19---20-.
     |                                                                          
              |
     |                                                                          
              |
     
.----------------------------------------------------------------------------------------.

Starting at node 1, their neighbors are 2 and 20, their neighbors 3
and 19, and so on...

BTW,  I have a  fortran routine and a c++  wrapper that  I'm using for
solving  Finite Element matrices. It's   based on the 'uactcl' skyline
solver that appears in the book  of Zienkiewicz. I  can pass you it if
it's useful to you.  Off course, I'm  interested in providing this  to
the rest but I'm brand new to c++ programming. Just  a week ago Michel
Hanke sent me his wrappers for the ode's and dae's solvers (those that
he has sent today to octave-sources), and studying them I succeeded in
writing the wrapper.

Regards,

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: address@hidden, http://venus.unl.edu.ar/gtm-eng.html |



reply via email to

[Prev in Thread] Current Thread [Next in Thread]