[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Banded matrix with dirty offdiagonal triagles
From: 
Mario Storti 
Subject: 
Re: Banded matrix with dirty offdiagonal 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 wellexamined case, in my problem my matrix also has "dirty"
> (nonzero) triangles in the upperright and lowerleft 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 McKee algorithm. On
the graph. The graph of the matrix is something like this:
1234567891011121314151617181920.




..
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 octavesources), 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/gtmeng.html 