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:
>
> 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.
>
> [...]
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[[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],
{p, 2, pcons}];
maxnorm]
In[49]:= Timing[func[14]]
Out[49]= {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
comments stripped out.
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[[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]];
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}];
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[21]:= Timing[compiledfunc[14]]
Out[21]= {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[22]:= Timing[compiledfunc[18]]
Out[22]= {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
- References:
- Re: Why is Mathematica so slow ?
- From: Daniel Lichtblau <danl@wolfram.com>
- Re: Why is Mathematica so slow ?