       Correction and improvement (Re: Re: Why is Mathematica so slow ?)

• To: mathgroup at smc.vnet.net
• Subject: [mg25478] Correction and improvement (Re: [mg25461] Re: [mg25415] Why is Mathematica so slow ?)
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Mon, 2 Oct 2000 22:26:57 -0400 (EDT)
• References: <200009290507.BAA18098@smc.vnet.net> <200010010644.CAA29727@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```It appears I need to correct my own post. Thought I'd show some nice
speed improvement while I'm at it.

Daniel Lichtblau 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,
> >
> >
> > (* 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 = 0;
> >       beta = difflower[]/diffdiag[];
> >       alpha = diffdiag[];
> >
> >       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 = 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.
>
> [...]

The code I previously sent had been through several refinements and
apparently I posted an erroneous version that I had not tested. The
problem was that the diffdiag and diffupper vectors need to be
initialized. I also thought it preferable to put everything into a
Module and return the maxnorms in a list. Here is a variation that does
all this.

a0 = 0;
b0 = 1;
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;

func[pcons_Integer] := Module[
{a0=0, b0=1, meshsize, number, diffdiag, diffupper, difflower,
alpha, beta, gamma, x, y, ulist, maxnorm},
maxnorm = Table[0., {pcons}];
Do [
meshsize = h[p];
number = n[p];
(*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.*)
diffdiag = Table[0., {number}];
diffupper = Table[0., {number}];
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;
(*Calculation of alpha, beta and gamma parameters from the
original eigenvalue equation.*)
gamma = diffupper;
alpha = Table[0., {number}];
beta = Table[0., {number}];
beta[] = difflower[]/diffdiag[];
alpha[] = diffdiag[];
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[] = 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],
{p, 2, pcons}];
maxnorm]

In:= Timing[func]

Out= {288.56 Second, {0., 0.0165198, 0.00415479, 0.00105684,
0.000264318,
-6            -6
-7
>     0.0000660863, 0.000016522, 4.13053 10  , 1.03263 10  , 2.58158 10  ,
-8            -8            -9            -10
>     6.45409 10  , 1.61356 10  , 3.08023 10  , 2.81992 10   }}

Today I realized we can get it to run alot faster if we make judicious
use of the Compile function. Below is some code to do it this way, with

Clear[f,g,h,n,phi];

g = Compile[{{x,_Real}}, x^2];
f = Compile[{{x,_Real}}, (1 + 4 x + 2 x^2 - x^4)*Exp[x]];
phi = Compile[{{x,_Real}}, (1 - x^2)*Exp[x]];
h = Compile[{{p,_Integer}}, 1/(2.^p)];
n = Compile[{{p,_Integer},{a0,_Integer},{b0,_Integer}},
(b0 - a0)*2^p - 1];

compiledfunc = Compile[{{pcons,_Integer}},
Module[{a0=0, b0=1, meshsize, number, diffdiag,
diffupper, difflower, alpha, beta, gamma, x, y, ulist, maxnorm},
maxnorm = Table[0., {pcons}];
Do [
meshsize = h[p];
number = n[p,a0,b0];
diffdiag = Table[0., {number}];
diffupper = Table[0., {number}];
Do [
diffdiag[[k]] = 2/meshsize^2 + g[a0 + k*meshsize];
diffupper[[k]] = -1/meshsize^2,
{k, 1, number-1}];
diffdiag[[number]] = 2/meshsize^2 + g[a0 + number*meshsize];
difflower = diffupper;
gamma = diffupper;
alpha = Table[0., {number}];
beta = Table[0., {number}];
beta[] = difflower[]/diffdiag[];
alpha[] = diffdiag[];
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]];
y = Table[0., {number}];
y[] = f[a0 + meshsize] + 1/meshsize^2;
Do [y[[i]] = f[a0 + i*meshsize] - beta[[i]]*y[[i-1]], {i,2,number}];
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],
{p, 2, pcons}];
maxnorm],
{{f[_],_Real}, {g[_],_Real}, {h[_],_Real},
{phi[_],_Real}, {n[_],_Integer}}]

In:= Timing[compiledfunc]

Out= {3.65 Second, {0., 0.0165198, 0.00415479, 0.00105684,
0.000264318,

-6            -6
-7
>     0.0000660863, 0.000016522, 4.13053 10  , 1.03263 10  , 2.58158 10  ,

-8            -8            -9           -10
>     6.45409 10  , 1.61356 10  , 3.08005 10  , 2.8163 10   }}

That seemed like hardly a test. We'll try 18 instead.

In:= Timing[compiledfunc]

Out= {58.35 Second, {0., 0.0165198, 0.00415479, 0.00105684,
0.000264318,

-6            -6
-7
>     0.0000660863, 0.000016522, 4.13053 10  , 1.03263 10  , 2.58158 10  ,

-8            -8            -9           -10
-10
>     6.45409 10  , 1.61356 10  , 3.08005 10  , 2.8163 10   , 3.0005 10   ,

-10            -10           -9
>     4.87473 10   , 4.90901 10   , 6.6502 10  }}

Still runs in under a minute. Not so slow, eh? ("Not too shaggy", as my
son once said).

Daniel Lichtblau
Wolfram Research

```

• Prev by Date: Re: Why is Mathematica so slow ?
• Next by Date: interpolation
• Previous by thread: Re: Why is Mathematica so slow ?
• Next by thread: Re: Why is Mathematica so slow ?