Re: Why is Mathematica so slow ?
- To: mathgroup at smc.vnet.net
- Subject: [mg25461] Re: [mg25415] Why is Mathematica so slow ?
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sun, 1 Oct 2000 02:44:45 -0400 (EDT)
- References: <200009290507.BAA18098@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Madhusudan Singh wrote:
>
> Hi
>
> I am solving a problem in a numerical linear algebra course that
> involves solving a two point boundary value problem with finite element
> method and LU factorisation. Since this is a part of a course, I am not
> permitted to use inbuilt functions in Mathematica.
> The matrix generated is a tridiagonal matrix. It is sparse. The
> order varies from 1 to 2^14. I solved this problem by using the
> following code (end of posting).
> I did the same problem on C (since I am using Mathematica purely as
> a simple programming language, it does not matter if I use C instead.)
>
> To solve for matrix orders 1 through 2^14 (14 steps) it takes about 4
> minutes to solve the problem with C, and more than 7 hours with
> Mathematica. The two codes are functionally identical. What is going on
> ?
>
> With regards,
>
> Madhusudan Singh.
>
> (* The definition and initialisation section.*)
> Clear["'*"]; Off[
> Part::"pspec"]; Off[General::"spell1"]; Off[Part::"partw"];
> a0 := 0;
> b0 := 1;
> pcons := 14;
> g[x_] := x^2;
> f[x_] := (1 + 4 x + 2 x^2 - x^4) Exp[x];
> phi[x_] := (1 - x^2) Exp[x];
> h[p_] := 1/(2^(p));
> n[p_] := (b0 - a0)/h[p] - 1;
> Array[maxnorm, pcons];
> Print["h(p) ||uh-phih|| ||uh-phih||/h2"];
>
> (*The loop over all p's."*)
> Do[{
> Clear[meshsize]; Clear[number];
> meshsize = h[p];
> number = n[p];
> Clear[a]; Clear[b]; Clear[c];
> Array[a, {number}];
> Array[b, {number}];
> Array[c, {number}];
> Clear[diffupper]; Clear[difflower]; Clear[diffdiag];
>
> (*Definition of the three diagonals.
> The matrix definition is eschewed as the matrix is sparse
> and \
> (for larger p's) memory can become an issue.*)
>
> Do[{a[k] := N[2/(meshsize)^(2) + g[a0 + k meshsize]];
> b[k + 1] := N[-1/(meshsize)^(2)];
> c[k] := N[-1/(meshsize)^(2)];}, {k, 1, number - 1}];
> a[number] := N[2/(meshsize)^(2) + g[a0 + number meshsize]];
> diffdiag := Table[a[k], {k, 1, number}];
> diffupper := Table[c[k], {k, 1, number - 1}];
> difflower := diffupper;
>
> Clear[alpha]; Clear[beta]; Clear[gamma];
>
> Array[gamma, number - 1];
> Array[alpha, number];
> Array[beta, number]; i := 1;
>
> (*Calculation of alpha,
> beta and gamma parameters from the original eigenvalue
> equation.*)
> \
>
> Do[{gamma[i] := diffupper[[i]]; }, {i, 1, number - 1}];
> beta[1] = 0;
> beta[2] = difflower[[1]]/diffdiag[[1]];
> alpha[1] = diffdiag[[1]];
>
> Do[{alpha[i] = diffdiag[[i]] - diffupper[[i - 1]] beta[i];
> beta[i + 1] = difflower[[i]]/alpha[i];}, {i, 2, number - 1}];
> alpha[number] =
> diffdiag[[number]] - diffupper[[number - 1]] beta[number];
>
> (*Forward elimination.*)
> Clear[y]; Clear[x];
> Array[y, number];
> y[1] = f[a0 + 1 meshsize] + 1/(meshsize)^(2);
> Do[{y[i] = f[a0 + i meshsize] - beta[i] y[i - 1];}, {i, 2,
> number}];
>
> (*Backward substitution.*)
> Clear[ulist];
>
> Array[x, number];
> x[number] = y[number]/alpha[number];
> Do[{x[i] = (y[i] - gamma[i] x[i + 1])/alpha[i];}, {i, number - 1,
> 1, -1}];
> ulist := Table[Abs[x[i] - phi[a0 + i meshsize]], {i, 1, number}];
> maxnorm[p] = Max[ulist];
> resultstream =
> OpenAppend["result13.dat", FormatType -> OutputForm,
> PageWidth -> Infinity];
> Write[resultstream, N[meshsize], " ", maxnorm[p], " ",
> maxnorm[p]/(meshsize)^(2)];
> Print[N[meshsize], " ", maxnorm[p], " ",
> maxnorm[p]/(meshsize)^(2)];
> Close[resultstream];}, {p, 1, pcons}];
>
> (*Plotting the results*)
>
> maxnormlist := Table[{h[p], maxnorm[p]}, {p, 1, pcons}];
> algorithmefficiency := Table[{h[p], maxnorm[p]/(h[p])^(2)}, {p, 1,
> pcons}];
>
> Display["maxnorm.eps", ListPlot[maxnormlist, PlotJoined -> True],
> "EPS"];
> Display["algoeff.eps", ListPlot[algorithmefficiency, PlotJoined ->
> True],
> "EPS"];
I did not try to run the code as written above so I cannot say with
certainty why it might be slow. It uses alot of hash table look-up when
instead it can index into numeric vectors, and this latter is alot
faster so I rewrote it to do just that.
You never appear to use b so I got rid of it. Also I got rid of a and c
in favor of the diffXXX vectors. I replaced the file I/O with a Print
statement and added another, just to get quick visual feedback as to
whether things appear to progress as they should.
a0 = 0;
b0 = 1;
pcons = 14;
g[x_] := x^2;
f[x_] := (1 + 4 x + 2 x^2 - x^4) Exp[x];
phi[x_] := (1 - x^2) Exp[x];
h[p_] := 1/(2^p);
n[p_] := (b0 - a0)/h[p] - 1;
Timing[Do [
meshsize = h[p];
number = n[p];
Print["number is ", number];
(*Definition of the three diagonals.
The matrix definition is eschewed as the matrix is sparse
and (for larger p's) memory can become an issue.*)
Do [
diffdiag[[k]] = N[2/(meshsize)^(2) + g[a0 + k meshsize]];
diffupper[[k]] = N[-1/(meshsize)^(2)],
{k, 1, number-1}];
diffdiag[[number]] = N[2/(meshsize)^(2) + g[a0 + number meshsize]];
difflower = diffupper;
i = 1;
(*Calculation of alpha, beta and gamma parameters from the
original eigenvalue equation.*)
gamma = diffupper;
alpha = Table[0., {number}];
beta = Table[0., {number}];
beta[[2]] = difflower[[1]]/diffdiag[[1]];
alpha[[1]] = diffdiag[[1]];
Do [
alpha[[i]] = diffdiag[[i]] - diffupper[[i - 1]] beta[[i]];
beta[[i + 1]] = difflower[[i]]/alpha[[i]],
{i, 2, number-1}];
alpha[[number]] =
diffdiag[[number]] - diffupper[[number - 1]] beta[[number]];
(*Forward elimination.*)
y = Table[0., {number}];
y[[1]] = f[a0 + meshsize] + 1/(meshsize)^(2);
Do [y[[i]] = f[a0 + i meshsize] - beta[[i]] y[[i - 1]], {i, 2,
number}];
(*Backward substitution.*)
x = Table[0., {number}];
x[[number]] = y[[number]]/alpha[[number]];
Do[x[[i]] = (y[[i]] - gamma[[i]] x[[i + 1]])/alpha[[i]],
{i, number-1, 1, -1}];
ulist = Abs[x - Table[phi[a0 + i meshsize], {i, 1, number}]];
maxnorm[p] = Max[ulist];
Print[N[meshsize], " ", maxnorm[p], " ", maxnorm[p]/(meshsize)^(2)],
{p, 1, pcons}];]
When I run this I see that the p==1 case is a bit problematic. The rest
are fine. It takes a few minutes on my 300Mhz machine, running a
development version of Mathematica under Linux.
number is 1
Set::partw: Part 2 of {0.} does not exist.
55 Sqrt[E] 55
Sqrt[E]
4 + ---------- 4 +
----------
-3 Sqrt[E] 16 -3 Sqrt[E] 16
0.5 Abs[---------- + --------------] 4 Abs[---------- +
--------------]
4 8.25 + 0. List 4 8.25 + 0.
List
number is 3
0.25 0.0165198 0.264317
number is 7
0.125 0.00415479 0.265907
number is 15
0.0625 0.00105684 0.270552
number is 31
0.03125 0.000264318 0.270662
number is 63
0.015625 0.0000660863 0.27069
number is 127
0.0078125 0.000016522 0.270696
number is 255
-6
0.00390625 4.13053 10 0.270698
number is 511
-6
0.00195313 1.03263 10 0.270699
number is 1023
-7
0.000976563 2.58158 10 0.270699
number is 2047
-8
0.000488281 6.45409 10 0.270704
number is 4095
-8
0.000244141 1.61356 10 0.27071
number is 8191
-9
0.00012207 3.08023 10 0.206711
number is 16383
-10
0.0000610352 2.81992 10 0.0756968
Out[80]= {266.12 Second, Null}
The code below likewise behaves just fine.
maxnormlist = Table[{h[p], maxnorm[p]}, {p, 1, pcons}];
algorithmefficiency = Table[{h[p], maxnorm[p]/(h[p])^(2)}, {p, 1,
pcons}];
ListPlot[maxnormlist, PlotJoined -> True]
ListPlot[algorithmefficiency, PlotJoined -> True]
I have not looked hard at the equation solving part either in regards to
correctness or efficiency. In general as regards the latter I'd not be
surprised if there are further small improvements to be made without
resorting to built-in linear algebra functionality.
As a general remark, I doubt ordinary LU factorization is appropriate
for finite difference or finite element type problems. When possible
you'd want a dedicated banded (e.g. tridiagonal) solver, as above. More
generally a sparse solver seems like a better way to go than dense LU
factorization. In version 4 of Mathematica there is a sparse variant,
Developer`SparseLinearSolve, which can be useful for such problems.
Daniel Lichtblau
Wolfram Research
- Follow-Ups:
- Correction and improvement (Re: Re: Why is Mathematica so slow ?)
- From: Daniel Lichtblau <danl@wolfram.com>
- Correction and improvement (Re: Re: Why is Mathematica so slow ?)