Re: Why is Mathematica so slow ?
- To: mathgroup at smc.vnet.net
- Subject: [mg25781] Re: Why is Mathematica so slow ?
- From: Robert Knapp <rknapp at wolfram.com>
- Date: Wed, 25 Oct 2000 03:53:42 -0400 (EDT)
- Organization: Wolfram Research, Inc.
- References: <8r1a9q$i50@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
If I surmise correctly that the main complexity in your code is solving a tridiagonal linear system, then it can be greatly speeded up, even if you cannot use the Mathematica SparseLinearSolve command. I showed how to do this in a talk last year. You can find the notebook at http://members.wri.com/rknapp/Talks/DeveloperConference99/Tridiagonal.nb 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"];