help-octave
[Top][All Lists]

## 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,

> 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.

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                           |