MathGroup Archive 2000

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: PentadiagonalSolve?

  • To: mathgroup at
  • Subject: [mg22362] Re: PentadiagonalSolve?
  • From: Daniel Lichtblau <danl at>
  • Date: Fri, 25 Feb 2000 21:14:16 -0500 (EST)
  • Organization: Wolfram Research, Inc.
  • References: <88lbad$>
  • Sender: owner-wri-mathgroup at

Kurt Taretto wrote:
> Solving a system of linear equations with a matrix with 5 diagonals is
> done pretty well by LinearSolve.  But since I work with large matrices,
> and all of them have only 5 nonzero diagonals, I don't have to use
> memory for the nonzero entries in order to give LinearSolve a "good and
> square" matrix.
> Does anybody know how could I, for example, do some transformation from
> 5 to 3 diagonals and then use TridiagonalSolve?  I would
> appreciate any help.  Thanks!

If you require only machine arithmetic this can be done effectively
using the sparse linear solver in the Developer context (version 4
only). Below is an example.

First we create a diagonally-dominant random pentadiagonal sparse
matrix. It has the form { {row,col}->value} ... }.

random5Mat[n_] := Module[{mat},
	mat = Table[{j+If[k>0,k,0],j+If[k<0,-k,0]}->Random[],
	  {k,-2,2}, {j,n-Abs[k]}];
	mat[[3]] = Map[#[[1]]->#[[2]]+1&, mat[[3]]];

Now we'll run an example of dimension 1000.

m = 1000;
sp5mat = random5Mat[m];
rhs = Table[Random[],{m}];

On my machine (300 mhz, Linux OS):

In[198]:= Timing[sol1 = Developer`SparseLinearSolve[sp5mat, rhs];]
Out[198]= {0.06 Second, Null}

We can compare to LinearSolve by converting to an ordinary matrix, as

denseMat[smat_] := Module[{dim, res},
	dim = Max[Map[#[[1]]&,smat]];
	res = Table[{j,k}, {j,dim}, {k,dim}] /. smat;
	Map[If[ListQ[#],0,#]&, res, {2}]

In[203]:= Timing[dens5mat = denseMat[sp5mat];]
Out[203]= {56. Second, Null}

In[204]:= Timing[sol2 = LinearSolve[dens5mat, rhs];]
Out[204]= {3.36 Second, Null}

Finally we check that the results are in fairly close agreement.

In[205]:= Max[Abs[sol2-sol1]]
Out[205]= 4.1922 10

When I do the sparse solving at dimension 10000, I get

Out[210]= {1.87 Second, Null}

Daniel Lichtblau
Wolfram Research

  • Prev by Date: Re: Building lists
  • Next by Date: Re: Building lists
  • Previous by thread: PentadiagonalSolve?
  • Next by thread: Database Access Kit