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>