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