Banded matrix with dirty off-diagonal triagles

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

Banded matrix with dirty off-diagonal triagles

Michael Smolsky-3
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.


Reply | Threaded
Open this post in threaded view
|

Re: Banded matrix with dirty off-diagonal triagles

Mario Storti-4

>>>>> On Wed, 17 Jun 1998 19:05:12 +0300,
>>>>>      Michael Smolsky <[hidden email]> 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: [hidden email], http://venus.unl.edu.ar/gtm-eng.html |


Reply | Threaded
Open this post in threaded view
|

Re: Banded matrix with dirty off-diagonal triagles

Mario Storti-4
In reply to this post by Michael Smolsky-3

My apologies.  In  my previous posting,  I forgot  to  include the off
diagonal terms when showing the matrix. I missed to show the following
steps, but the rest of the mail, and the renumbering is still valid!

octave> A(1,20)=-1;
octave> A(20,1)=-1;
octave> A(1,20)=-1;
octave> A
A =

   2  -1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  -1
  -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
  -1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0  -1   2

Mario


Reply | Threaded
Open this post in threaded view
|

Re: Banded matrix with dirty off-diagonal triagles

Michael Hanke-2
In reply to this post by Michael Smolsky-3
I am using UMFPACK from netlib.
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: Banded matrix with dirty off-diagonal triagles

Mario Storti-4

>>>>> On Thu, 18 Jun 1998 14:30:55 +0200 (MET DST),
>>>>>      Michael Hanke <[hidden email]> said:

> I am using UMFPACK from netlib.
> 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]       |
> +---------------------------------------------------------------+

Seems fine. However, it seems to me that their copyright notice is not
compatible with the GPL, if what we want is a sparse module coming
standard with Octave. Am I wrong?

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 |


================================================================

> UMFPACK Version 2.2:  Unsymmetric-pattern Multifrontal Package
> --------------------------------------------------------------
>
> Authors:  Timothy A. Davis and Iain S. Duff.  Copyright (C) 1997.
> Date:     July, 1997.
>
>
> ***********************************************************************
> * NOTICE:  "The UMFPACK Package may be used SOLELY for educational,   *
> * research, and benchmarking purposes by non-profit organizations and *
> * the U.S. government.  Commercial and other organizations may make   *
> * use of UMFPACK SOLELY for benchmarking purposes only.  UMFPACK may  *
> * be modified by or on behalf of the User for such use but at no time *
> * shall UMFPACK or any such modified version of UMFPACK become the    *
> * property of the User.  UMFPACK is provided without warranty of any  *
> * kind, either expressed or implied.  Neither the Authors nor their   *
> * employers shall be liable for any direct or consequential loss or   *
> * damage whatsoever arising out of the use or misuse of UMFPACK by    *
> * the User.  UMFPACK must not be sold.  You may make copies of        *
> * UMFPACK, but this NOTICE and the Copyright notice must appear in    *
> * all copies.  Any other use of UMFPACK requires written permission.  *
> * Your use of UMFPACK is an implicit agreement to these conditions."  *
> *                                                                     *


Reply | Threaded
Open this post in threaded view
|

Re: Banded matrix with dirty off-diagonal triagles

John W. Eaton-6
On 18-Jun-1998, Mario Storti <[hidden email]> wrote:

| >>>>> On Thu, 18 Jun 1998 14:30:55 +0200 (MET DST),
| >>>>>      Michael Hanke <[hidden email]> said:
|
| > I am using UMFPACK from netlib.
| > Michael

| Seems fine. However, it seems to me that their copyright notice is not
| compatible with the GPL, if what we want is a sparse module coming
| standard with Octave. Am I wrong?

No, you're not wrong.  The copyright terms of UMFPACK clearly conflict
with the GPL, so it can't be linked with Octave.

Thanks,

jwe


Reply | Threaded
Open this post in threaded view
|

Re: Banded matrix with dirty off-diagonal triagles

mhagger-2
> No, you're not wrong.  The copyright terms of UMFPACK clearly conflict
> with the GPL, so it can't be linked with Octave.
Actually if you are looking for an excellent Sparse matrix solving
package, somewhat like UMFPACK only with no restrictions, then you might
want to take a look at either SuperLU or SPOOLES (SPOOLES will run in
parallel as well, although this is a little more involved to use).  Both
packages are available via the High performance Netlib web page (a link to
this is at the bottom on the main netlib page).

I was forced to abandon using UMFPACK some time ago, partly due to the
license restrictions and partly as you had to know in advance how much
memory the sparse LU needed, which was not always very convenient when
memory was tight.

Mark

--

Note that these views in no way represent those of my employer.

--