Re: PentadiagonalSolve?
- To: mathgroup at smc.vnet.net
- Subject: [mg22362] Re: PentadiagonalSolve?
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Fri, 25 Feb 2000 21:14:16 -0500 (EST)
- Organization: Wolfram Research, Inc.
- References: <88lbad$bso@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
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]]];
Flatten[mat]
]
Now we'll run an example of dimension 1000.
SeedRandom[1111];
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
below.
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]]
-12
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