(*********************************************************************** Mathematica-Compatible Notebook This notebook can be used on any computer system with Mathematica 4.0, MathReader 4.0, or any compatible application. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. ***********************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 34417, 979]*) (*NotebookOutlinePosition[ 35071, 1002]*) (* CellTagsIndexPosition[ 35027, 998]*) (*WindowFrame->Normal*) Notebook[{ Cell[CellGroupData[{ Cell["Packed Arrays", "Title", TextAlignment->Center], Cell[TextData[{ "An efficient internal storage format in ", StyleBox["Mathematica", FontSlant->"Italic"], " 4" }], "Subtitle", TextAlignment->Center], Cell[TextData[{ "While packed arrays may sound quite mysterious and deep, the idea is \ really quite simple. When it is clearly appropriate to represent a list of a \ single type machine numbers (integer, real, or complex) internally as an \ array, ", StyleBox["Mathematica", FontSlant->"Italic"], " does so. The basic technology is actually there in ", StyleBox["Mathematica", FontSlant->"Italic"], " 3.0, an is visible in the List operations which were added to the \ compiler. Basically what we have done for ", StyleBox["Mathematica", FontSlant->"Italic"], " 4 is to extend this representation to kernel functions where it is \ appropriate. The internal representation is essentially invisible to you as \ a user (you have to use functions in the Developer` context to even see it) \ so you do not have to worry about it. However, understanding the packed \ array representation can help you write faster more memory efficient code. \ This part of my presentation is intended to give you the information you \ need." }], "Text"], Cell[CellGroupData[{ Cell["Representations of data in Mathematica", "Section"], Cell["If I enter a list in Mathematica", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(list\ = \ {1, 2, 3, 4}\)], "Input"], Cell[BoxData[ \({1, 2, 3, 4}\)], "Output"] }, Open ]], Cell["\<\ its internal representation may be different from the OutputForm. The basic \ common representation for lists (and all expressions) is given by the \ FullForm\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(FullForm[list]\)], "Input"], Cell[BoxData[ TagBox[ StyleBox[\(List[1, 2, 3, 4]\), ShowSpecialCharacters->False, ShowStringCharacters->True, NumberMarks->True], FullForm]], "Output"] }, Open ]], Cell[TextData[{ "each integer you see in the list is itself an expression. Internally in \ the ", StyleBox["Mathematica", FontSlant->"Italic"], " C code, the data stored carries information about the type of the \ expression (in this case, list, or number which is a machine integer)" }], "Text"], Cell["If I enter the data a different way,", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(packed\ = \ Range[4]\)], "Input"], Cell[BoxData[ \({1, 2, 3, 4}\)], "Output"] }, Open ]], Cell["\<\ It is internally stored in a different way. The basic print forms are all \ the same. To all intents and puposes it is the same. \ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(packed\ === \ list\)], "Input"], Cell[BoxData[ \(True\)], "Output"] }, Open ]], Cell["\<\ You can only tell that it is stored in a different way by using a Developer` \ context function:\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Developer`PackedArrayQ\ /@ \ {list, \ packed}\)], "Input"], Cell[BoxData[ \({False, True}\)], "Output"] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell["Advantages of packed arrays", "Section"], Cell["\<\ There are two advantages of packed arrays, memory and speed. \ \>", "Text"], Cell[" The memory is the most obvious of these:", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(ByteCount[a\ = \ N[Range[1000]]]\)], "Input"], Cell[BoxData[ \(8056\)], "Output"] }, Open ]], Cell["\<\ If we change the first element to be a different type, the result cannot be \ stored as a packed array:\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(ByteCount[b\ = \ a; \ b[\([1]\)]\ = \ 1; \ b]\)], "Input"], Cell[BoxData[ \(20024\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(N[%/%%]\)], "Input"], Cell[BoxData[ \(2.4856007944389273`\)], "Output"] }, Open ]], Cell["The speed advantage typically comes in common operations:", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[Do[a + a, {1000}]]\)], "Input"], Cell[BoxData[ \({0.22000000000116415`\ Second, Null}\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[Do[b + b, {1000}]]\)], "Input"], Cell[BoxData[ \({4.120000000002619`\ Second, Null}\)], "Output"] }, Open ]], Cell["\<\ Understanding why the timing difference occurs is the key to knowing when \ packed arrays will help your code run faster.\ \>", "Text"], Cell[TextData[{ "In the first operation, a is stored as a packed array, so when Plus see \ it, it knows all of the elements are machine real numbers, so internally, a \ function which adds arrays of real numbers is added. Except for overflow and \ underflow checking which ", StyleBox["Mathematica", FontSlant->"Italic"], " does, this is done as fast as you could get by wrting a loop to do the \ addition in a C-program for two arrays of machine numbers." }], "Text"], Cell["\<\ In the second operation, b is stored as List[..], so what happens is first \ that Plus is threaded over the lists, producing (roughly) {b[[1]] + b[[1]], b[[2]] + b[[2]], .... } i.e. Plus has to process 1000 expressions. For each pair of numbers, plus \ has to use the information attached to those numbers (i.e real, integer, \ machine, bignum, etc) to determine what function to use to do the actual \ adding. This all takes time and accounts for nearly all of the difference \ here.\ \>", "Text"], Cell["\<\ There are different ways that packed arrays can speed things up, but this \ exemplifies the type which typically gives the greatest speedup.\ \>", "Text"] }, Closed]], Cell[CellGroupData[{ Cell["When do you get packed arrays?", "Section"], Cell["\<\ When trying to write efficient programs, this is a very important question. \ There are three important situtations which commonly produce packed arrays.\ \>", "Text"], Cell[CellGroupData[{ Cell["\<\ 1. Listable and structural operations with packed array input (where the \ output can be a packed array.)\ \>", "Subsubsection"], Cell["\<\ By listable operations, I mean those with the attribute Listable. \ \>", "Text"], Cell[BoxData[ \(Select[Names["\"], \ MemberQ[Attributes[#], \ Listable] &]\)], "Input"], Cell["\<\ the majority of which will produce packed arrays given packed array \ arguments.\ \>", "Text"], Cell["\<\ By structural operations, I mean commands like Part, Join, Transpose, etc. \ which act on the structure of an expression.\ \>", "Text"] }, Closed]], Cell[CellGroupData[{ Cell["\<\ 2. Functions which use an internal array format for machine numbers.\ \>", "Subsubsection"], Cell["Fourier is a good example of this:", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Developer`PackedArrayQ[Fourier[{1. , 2. , 3. , 4. }]]\)], "Input"], Cell[BoxData[ \(True\)], "Output"] }, Open ]], Cell["\<\ Whether the input is a packed array or not, if Fourier does its computation \ with machine numbers, it internally uses arrays already, so it is natural to \ give the output in terms of a packed array.\ \>", "Text"] }, Closed]], Cell[CellGroupData[{ Cell["\<\ 3. Functions which can quickly determine from simple input that the output \ could be a packed array.\ \>", "Subsubsection"], Cell["\<\ I have already used on example in this class, Range. In Range[4], it is \ pretty easy to see before generating any numbers that the result can be \ represented as a packed array.\ \>", "Text"], Cell["For something like", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[ Developer`PackedArrayQ[ data\ = \ Table[Random[], {10000}, {2}]]]\)], "Input"], Cell[BoxData[ \({0.050000000000000044`\ Second, True}\)], "Output"] }, Open ]], Cell["\<\ it is a little more difficult to determine. Commands like this are handled \ by auto-compilation. However, since trying the compiler does take some time, \ there is a (user modifiable) threshold length which must be exceeded before \ compilation is attempted.\ \>", "Text"] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell["Can packed arrays ever be slower?", "Section"], Cell["Yes. Here is an example", "Text"], Cell[BoxData[ \(f[{a_, b_}]\ := \ a^2\ - \ b^2\)], "Input"], Cell[BoxData[ \(Timing[Developer`PackedArrayQ[Map[f, data]]]\)], "Input"], Cell[BoxData[ \(Timing[\(udata\ = \ Developer`FromPackedArray[data];\)]\)], "Input"], Cell[BoxData[ \(Timing[\(Map[f, udata];\)]\)], "Input"], Cell["\<\ However, if we inline the function, the auto-compilation can work:\ \>", "Text"], Cell[BoxData[ \(Timing[ Developer`PackedArrayQ[ Map[\((\((#[\([1]\)])\)^2\ + \ \((#[\([2]\)])\)^2)\) &, \ data]]]\)], "Input"], Cell["\<\ Now it is not always the case that things are so easily patched up: some \ commands are incompatible with packed arrays. However, before complaining \ too loudly about these (we hope quite rare) cases, please take into account \ the whole time going into your computation. In the above example, we were \ not taking into account the time it too to generate the random list. Without \ auto-compilation and packed arrays, it takes\ \>", "Text"], Cell[BoxData[{ \(\(Developer`SetSystemOptions["\" -> Infinity];\)\), "\[IndentingNewLine]", \(Timing[ Developer`PackedArrayQ[ data\ = \ Table[Random[], {10000}, {2}]]]\)}], "Input"], Cell["which is quite a bit longer.", "Text"] }, Closed]], Cell[CellGroupData[{ Cell["How do I write code which takes advantage of packed arrays?", "Section"], Cell["\<\ The most important thing to do is to make sure that you do operations on \ lists of elements all at once or in groups whenever possible.\ \>", "Text"], Cell[TextData[{ "The good news is that this programming principle has not really changed \ from ", StyleBox["Mathematica", FontSlant->"Italic"], " 1.0: efficient code in ", StyleBox["Mathematica", FontSlant->"Italic"], " has always used functional and list based constructs." }], "Text"], Cell["\<\ The one other thing to keep in mind is that type needs to be consistent for \ there to even be a packed array representation.\ \>", "Text"], Cell["\<\ We have provided some aids which can help you in putting together your \ code:\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Developer`PackedArrayForm[Range[1000]]\)], "Input"], Cell[BoxData[ \("PackedArray"[Integer, "<"\[InvisibleSpace]1000\[InvisibleSpace]">"]\)], "Output"] }, Open ]], Cell[BoxData[ \(\(Developer`SetSystemOptions["\" \[Rule] True];\)\)], "Input"], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[Apply[Plus, \ N[Range[1000]]]]\)], "Input"], Cell[BoxData[ \(Developer`FromPackedArray::"punpack1" \(\(:\)\(\ \)\) "Unpacking array to level \!\(1\)."\)], "Message"], Cell[BoxData[ \({0.05000000000001137`\ Second, 500500.`}\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[Tr[N[Range[1000]]]]\)], "Input"], Cell[BoxData[ \({0.`\ Second, 500500.`}\)], "Output"] }, Open ]], Cell["\<\ Sometimes, you know your code will run more efficiently if you pack the input \ beforehand. For example\ \>", "Text"], Cell[BoxData[{ \(\(SeedRandom[1234];\)\), "\[IndentingNewLine]", \(\(data\ = \ Table[Random[Complex], {240}];\)\)}], "Input"], Cell[BoxData[ \(\(one[a_]\ := \ Max[Abs[1. \ - \ \((\ Sin[a]^2\ - \ Cos[a]^2)\)]];\)\)], "Input"], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[one[data]]\)], "Input"], Cell[BoxData[ \({0.049999999999954525`\ Second, 4.741196834993035`}\)], "Output"] }, Open ]], Cell[BoxData[ \(One[a_]\ := \ Module[{a1\ = \ Developer`ToPackedArray[a]}, one[a1]]\)], "Input"], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[One[data]]\)], "Input"], Cell[BoxData[ \({0.`\ Second, 4.741196834993035`}\)], "Output"] }, Open ]], Cell["Probably, the best thing I can do is show you an example", "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Example: LU Decomposition", "Section"], Cell[TextData[{ "Here is the code for actually doing an LU decomposition (factorization) \ from the package LinearAlgebra`GaussianElimination. Since this functionality \ exists in the kernel, the real purpose of this package is pedagogical, but I \ think the original author has gone too far with the intent to stick to \ standard C-looking code so that people familiar only with procedural code can \ understand it. The result is, in my opinion, a prime example of how you \ should not write code in ", StyleBox["Mathematica", FontSlant->"Italic"], ". Lets see why, and how we can improve it." }], "Text"], Cell[BoxData[ RowBox[{\(lufactor[aa_]\), ":=", RowBox[{"Module", "[", RowBox[{\({a = aa, pivot, ii, iip, i, ip, j, k, mpiv, m, n = Length[aa], tmp}\), ",", "\[IndentingNewLine]", RowBox[{\(pivot = Table[i, {i, n}]\), ";", "\[IndentingNewLine]", RowBox[{ StyleBox["For", FontColor->RGBColor[1, 0, 0]], StyleBox["[", FontColor->RGBColor[1, 0, 0]], RowBox[{ StyleBox[\(ii = 1\), FontColor->RGBColor[1, 0, 0]], StyleBox[",", FontColor->RGBColor[1, 0, 0]], StyleBox[\(ii \[LessEqual] n\), FontColor->RGBColor[1, 0, 0]], StyleBox[",", FontColor->RGBColor[1, 0, 0]], StyleBox[\(ii++\), FontColor->RGBColor[1, 0, 0]], ",", \( (*for\ each\ row\ do ... *) \), "\[IndentingNewLine]", \( (*find\ a\ pivot*) \), "\[IndentingNewLine]", RowBox[{\(mpiv = Abs[a[\([pivot[\([ii]\)], ii]\)]]\), ";", "\[IndentingNewLine]", \(k = ii\), ";", "\[IndentingNewLine]", RowBox[{ StyleBox["For", FontColor->RGBColor[1, 0, 0]], StyleBox["[", FontColor->RGBColor[1, 0, 0]], RowBox[{ StyleBox[\(i = ii + 1\), FontColor->RGBColor[1, 0, 0]], StyleBox[",", FontColor->RGBColor[1, 0, 0]], StyleBox[\(i \[LessEqual] n\), FontColor->RGBColor[1, 0, 0]], StyleBox[",", FontColor->RGBColor[1, 0, 0]], StyleBox[\(i++\), FontColor->RGBColor[1, 0, 0]], StyleBox[",", FontColor->RGBColor[1, 0, 0]], \(tmp = Abs[a[\([pivot[\([i]\)], ii]\)]]; \[IndentingNewLine]If[tmp > mpiv, mpiv = tmp; \ k = i]\)}], "]"}], ";", "\[IndentingNewLine]", \(tmp = pivot[\([ii]\)]\), ";", "\[IndentingNewLine]", \(pivot[\([ii]\)] = \(iip = pivot[\([k]\)]\)\), ";", "\[IndentingNewLine]", \(pivot[\([k]\)] = tmp\), ";", "\[IndentingNewLine]", \(mpiv = a[\([iip, ii]\)]\), ";", "\[IndentingNewLine]", \( (*calculate\ and\ store\ the\ \ multipliers\ and\ reduce\ the\ other\ \(\(elements\)\(.\)\)*) \), "\[IndentingNewLine]", RowBox[{ StyleBox["For", FontColor->RGBColor[1, 0, 0]], StyleBox["[", FontColor->RGBColor[1, 0, 0]], RowBox[{ StyleBox[\(i = ii + 1\), FontColor->RGBColor[1, 0, 0]], StyleBox[",", FontColor->RGBColor[1, 0, 0]], StyleBox[\(i \[LessEqual] n\), FontColor->RGBColor[1, 0, 0]], StyleBox[",", FontColor->RGBColor[1, 0, 0]], StyleBox[\(i++\), FontColor->RGBColor[1, 0, 0]], StyleBox[",", FontColor->RGBColor[1, 0, 0]], "\[IndentingNewLine]", RowBox[{\(ip = pivot[\([i]\)]\), ";", "\[IndentingNewLine]", \(m = \(a[\([ip, ii]\)] = a[\([ip, ii]\)]/mpiv\)\), ";", "\[IndentingNewLine]", RowBox[{ StyleBox["For", FontColor->RGBColor[1, 0, 0]], StyleBox["[", FontColor->RGBColor[1, 0, 0]], RowBox[{ StyleBox[\(j = ii + 1\), FontColor->RGBColor[1, 0, 0]], StyleBox[",", FontColor->RGBColor[1, 0, 0]], StyleBox[\(j \[LessEqual] n\), FontColor->RGBColor[1, 0, 0]], StyleBox[",", FontColor->RGBColor[1, 0, 0]], StyleBox[\(j++\), FontColor->RGBColor[1, 0, 0]], ",", \(a[\([ip, j]\)] -= m*a[\([iip, j]\)]\)}], "]"}]}]}], "]"}], ";"}]}], "]"}], ";", "\[IndentingNewLine]", \(LU[a, pivot]\)}]}], "]"}]}]], "Input"], Cell["\<\ Here are two test examples so that we can compare timings of different \ versions.\ \>", "Text"], Cell[BoxData[{ \(\(SeedRandom[1234];\)\), "\[IndentingNewLine]", \(\(machine\ = \ Table[Random[], {80}, {80}];\)\), "\[IndentingNewLine]", \(\(big\ = \ Table[Random[Real, {0, 1}, 20], {40}, {40}];\)\)}], "Input"], Cell[BoxData[ \(Timing[\(luf\ = \ lufactor[machine];\)]\)], "Input"], Cell[BoxData[ \(Timing[\(LUDecomposition[machine];\)]\)], "Input"], Cell[BoxData[ \(Timing[\(lufb\ = \ lufactor[big];\)]\)], "Input"], Cell[BoxData[ \(Timing[\(LUDecomposition[big];\)]\)], "Input"], Cell[TextData[{ "Ouch! We are a factor of 300 slower in the machine case and a factor of 5 \ slower in the bignum example. Before screaming an exasperated \"", StyleBox["Mathematica", FontSlant->"Italic"], " is just too slow\", lets try to improve the code." }], "Text"], Cell[TextData[{ "The first thing I saw is that this code is using For to handle the loops. \ In ", StyleBox["Mathematica", FontSlant->"Italic"], ", For is slower than alternative iterator constructs. Try to avoid using \ For in innermost loops unless it would make for more complicated code. In \ this example simply changing from For to Do makes a substantial difference." }], "Text"], Cell[BoxData[ RowBox[{\(lufactor1[aa_]\), ":=", RowBox[{"Module", "[", RowBox[{\({a = aa, pivot, ii, iip, i, ip, j, k, mpiv, m, n = Length[aa], tmp}\), ",", "\[IndentingNewLine]", RowBox[{\(pivot = Table[i, {i, n}]\), ";", "\[IndentingNewLine]", RowBox[{ StyleBox["Do", FontColor->RGBColor[0, 0, 1]], StyleBox["[", FontColor->RGBColor[0, 0, 1]], \( (*for\ each\ row\ do ... *) \), "\[IndentingNewLine]", \( (*find\ a\ pivot*) \), "\[IndentingNewLine]", RowBox[{ RowBox[{ StyleBox[\(mpiv = Abs[a[\([pivot[\([ii]\)], ii]\)]]\), FontColor->RGBColor[1, 0, 0]], StyleBox[";", FontColor->RGBColor[1, 0, 0]], StyleBox["\[IndentingNewLine]", FontColor->RGBColor[1, 0, 0]], StyleBox[\(k = ii\), FontColor->RGBColor[1, 0, 0]], StyleBox[";", FontColor->RGBColor[1, 0, 0]], StyleBox["\[IndentingNewLine]", FontColor->RGBColor[1, 0, 0]], StyleBox[ RowBox[{ StyleBox["Do", FontColor->RGBColor[0, 0, 1]], StyleBox["[", FontColor->RGBColor[0, 0, 1]], \(tmp = Abs[a[\([pivot[\([i]\)], ii]\)]]; \[IndentingNewLine]If[tmp > mpiv, mpiv = tmp; \ k = i], {i, ii + 1, n}\), "]"}], FontColor->RGBColor[1, 0, 0]], StyleBox[";", FontColor->RGBColor[1, 0, 0]], StyleBox["\[IndentingNewLine]", FontColor->RGBColor[1, 0, 0]], StyleBox[\(tmp = pivot[\([ii]\)]\), FontColor->RGBColor[1, 0, 0]], StyleBox[";", FontColor->RGBColor[1, 0, 0]], StyleBox["\[IndentingNewLine]", FontColor->RGBColor[1, 0, 0]], StyleBox[\(pivot[\([ii]\)] = \(iip = pivot[\([k]\)]\)\), FontColor->RGBColor[1, 0, 0]], StyleBox[";", FontColor->RGBColor[1, 0, 0]], StyleBox["\[IndentingNewLine]", FontColor->RGBColor[1, 0, 0]], StyleBox[\(pivot[\([k]\)] = tmp\), FontColor->RGBColor[1, 0, 0]], StyleBox[";", FontColor->RGBColor[1, 0, 0]], StyleBox["\[IndentingNewLine]", FontColor->RGBColor[1, 0, 0]], \(mpiv = a[\([iip, ii]\)]\), ";", "\[IndentingNewLine]", \( (*calculate\ and\ store\ the\ \ multipliers\ and\ reduce\ the\ \(\(other\)\(.\)\)*) \), "\[IndentingNewLine]", RowBox[{ StyleBox["Do", FontColor->RGBColor[0, 0, 1]], StyleBox["[", FontColor->RGBColor[0, 0, 1]], RowBox[{ RowBox[{\(ip = pivot[\([i]\)]\), ";", "\[IndentingNewLine]", \(m = \(a[\([ip, ii]\)] = a[\([ip, ii]\)]/mpiv\)\), ";", "\[IndentingNewLine]", RowBox[{ StyleBox["Do", FontColor->RGBColor[0, 0, 1]], StyleBox["[", FontColor->RGBColor[0, 0, 1]], RowBox[{\(a[\([ip, j]\)] -= m*a[\([iip, j]\)]\), ",", StyleBox[\({j, ii + 1, n}\), FontColor->RGBColor[0, 0, 1]]}], StyleBox["]", FontColor->RGBColor[0, 0, 1]]}]}], StyleBox[",", FontColor->RGBColor[0, 0, 1]], "\[IndentingNewLine]", StyleBox[\({i, ii + 1, n}\), FontColor->RGBColor[0, 0, 1]]}], StyleBox["]", FontColor->RGBColor[0, 0, 1]]}]}], StyleBox[",", FontColor->RGBColor[0, 0, 1]], "\[IndentingNewLine]", StyleBox[\({ii, n}\), FontColor->RGBColor[0, 0, 1]]}], StyleBox["]", FontColor->RGBColor[0, 0, 1]]}], ";", "\[IndentingNewLine]", \(LU[a, pivot]\)}]}], "]"}]}]], "Input"], Cell[BoxData[ \(Timing[lufactor1[machine]\ === \ luf]\)], "Input"], Cell[BoxData[ \(Timing[lufactor1[big]\ === \ lufb]\)], "Input"], Cell["\<\ The next change I will make will help make the code more readable. At this \ point, it won't seem to make much of a speed difference, but it does help. \ The idea is to replace the lines of code finding the correct pivot to use \ with a function.\ \>", "Text"], Cell[BoxData[ \(\(SetAttributes[FindPivot, \ HoldFirst];\)\)], "Input"], Cell[BoxData[ \(\(FindPivot[pivots_, \ a_, \ c_]\ := \ \[IndentingNewLine]Module[{k = First[MaxPosition[ Abs[a[\([Take[pivots, {c, \(-1\)}], c]\)]]]] + c - 1}, pivots[\([{c, k}]\)] = pivots[\([{k, c}]\)]; pivots[\([c]\)]];\)\)], "Input"], Cell["\<\ This finds which of the possible pivots (in column c, rows c and above) is \ largest in absolute value, and swaps with the one in row c.\ \>", "Text"], Cell["This, in turn, uses the simple function", "Text"], Cell[BoxData[ \(\(MaxPosition[x_]\ := \ First[Position[x, \ Max[x]]];\)\)], "Input"], Cell["And putting it into the main function", "Text"], Cell[BoxData[ RowBox[{\(lufactor2[aa_]\), ":=", RowBox[{"Module", "[", RowBox[{\({a = aa, pivot, ii, iip, i, ip, j, k, mpiv, m, n = Length[aa], tmp}\), ",", "\[IndentingNewLine]", RowBox[{\(pivot = Table[i, {i, n}]\), ";", "\[IndentingNewLine]", RowBox[{ "Do", "[", \( (*for\ each\ row\ do ... *) \), "\[IndentingNewLine]", RowBox[{ RowBox[{ StyleBox[\(iip\ = \ FindPivot[pivot, \ a, \ ii]\), FontColor->RGBColor[0, 0, 1]], StyleBox[";", FontColor->RGBColor[0, 0, 1]], "\[IndentingNewLine]", \(mpiv = a[\([iip, ii]\)]\), ";", "\[IndentingNewLine]", \( (*calculate\ and\ store\ the\ \ multipliers\ and\ reduce\ the\ \(\(other\)\(.\)\)*) \), "\[IndentingNewLine]", RowBox[{"Do", "[", RowBox[{ RowBox[{\(ip = pivot[\([i]\)]\), ";", "\[IndentingNewLine]", StyleBox[\(m = \(a[\([ip, ii]\)] = a[\([ip, ii]\)]/mpiv\)\), FontColor->RGBColor[0, 1, 0]], ";", "\[IndentingNewLine]", StyleBox[\(Do[ a[\([ip, j]\)] -= m*a[\([iip, j]\)], {j, ii + 1, n}]\), FontColor->RGBColor[1, 0, 0]]}], ",", "\[IndentingNewLine]", \({i, ii + 1, n}\)}], "]"}]}], ",", "\[IndentingNewLine]", \({ii, n}\)}], "]"}], ";", "\[IndentingNewLine]", \(LU[a, pivot]\)}]}], "]"}]}]], "Input"], Cell[BoxData[ \(Timing[lufactor2[machine]\ === \ luf]\)], "Input"], Cell[BoxData[ \(Timing[lufactor2[big]\ === \ lufb]\)], "Input"], Cell[TextData[{ "Still going the right direction, but no big savings. The place to look \ for big savings is, of course, the innermost loop. The next change replaces \ the innermost loop by giving ", StyleBox["Mathematica", FontSlant->"Italic"], " commands which operate on all the elements of a list:" }], "Text"], Cell[BoxData[ RowBox[{\(lufactor3[aa_]\), ":=", RowBox[{"Module", "[", RowBox[{\({a = aa, pivot, ii, iip, i, ip, j, mpiv, m, n = Length[aa]}\), ",", "\[IndentingNewLine]", RowBox[{\(pivot = Table[i, {i, n}]\), ";", "\[IndentingNewLine]", RowBox[{ "Do", "[", \( (*for\ each\ row\ do ... *) \), "\[IndentingNewLine]", RowBox[{ RowBox[{\(iip\ = \ FindPivot[pivot, \ a, ii]\), ";", "\[IndentingNewLine]", \(mpiv = a[\([iip, ii]\)]\), ";", "\[IndentingNewLine]", \( (*calculate\ and\ store\ the\ \ multipliers\ and\ reduce\ the\ \(\(other\)\(.\)\)*) \), "\[IndentingNewLine]", \(rc\ = \ Table[j, {j, ii + 1, n}]\), ";", "\[IndentingNewLine]", RowBox[{ StyleBox["Do", FontColor->RGBColor[1, 0, 0]], StyleBox["[", FontColor->RGBColor[1, 0, 0]], RowBox[{ RowBox[{ StyleBox[\(ip = pivot[\([i]\)]\), FontColor->RGBColor[1, 0, 0]], StyleBox[";", FontColor->RGBColor[1, 0, 0]], StyleBox["\[IndentingNewLine]", FontColor->RGBColor[1, 0, 0]], StyleBox[\(m = \(-\((a[\([ip, ii]\)] /= mpiv)\)\)\), FontColor->RGBColor[1, 0, 1]], StyleBox[";", FontColor->RGBColor[1, 0, 0]], StyleBox["\[IndentingNewLine]", FontColor->RGBColor[1, 0, 0]], StyleBox[\(a[\([ip, rc]\)] += m*a[\([iip, rc]\)]\), FontColor->RGBColor[0, 0, 1]]}], StyleBox[",", FontColor->RGBColor[0, 0, 1]], StyleBox["\[IndentingNewLine]", FontColor->RGBColor[1, 0, 0]], StyleBox[\({i, ii + 1, n}\), FontColor->RGBColor[1, 0, 0]]}], StyleBox["]", FontColor->RGBColor[1, 0, 0]]}]}], ",", "\[IndentingNewLine]", \({ii, n}\)}], "]"}], ";", "\[IndentingNewLine]", \(LU[a, pivot]\)}]}], "]"}]}]], "Input"], Cell[BoxData[ \(Timing[lufactor3[machine]\ === \ luf]\)], "Input"], Cell[BoxData[{ \(\(luf3\ = \ lufactor3[machine];\)\), "\[IndentingNewLine]", \({Max[Abs[First[luf3]\ - \ First[luf]]], luf3[\([2]\)]\ === \ luf[\([2]\)]}\)}], "Input"], Cell[BoxData[ \(Timing[lufactor3[big]\ === \ lufb]\)], "Input"], Cell["\<\ So replacing the innermost loop was spectacularly sucessful. There is still \ a set of nested loops which can be worked on though. \ \>", "Text"], Cell[BoxData[ RowBox[{\(lufactor4[aa_]\), ":=", RowBox[{"Module", "[", RowBox[{\({a = aa, pivot, ii, iip, i, ip, j, k, mpiv, m, n = Length[aa], row}\), ",", "\[IndentingNewLine]", RowBox[{\(pivot = Table[i, {i, n}]\), ";", "\[IndentingNewLine]", RowBox[{ "Do", "[", \( (*for\ each\ row\ do ... *) \), "\[IndentingNewLine]", RowBox[{ RowBox[{\(iip\ = \ FindPivot[pivot, \ a, \ ii]\), ";", "\[IndentingNewLine]", \(mpiv = a[\([iip, ii]\)]\), ";", "\[IndentingNewLine]", \( (*calculate\ and\ store\ the\ \ multipliers\ and\ reduce\ the\ \(\(other\)\(.\)\)*) \), "\[IndentingNewLine]", \(rc\ = \ Table[j, {j, ii + 1, n}]\), ";", "\[IndentingNewLine]", StyleBox[\(prc\ = \ pivot[\([rc]\)]\), FontColor->RGBColor[0, 0, 1]], StyleBox[";", FontColor->RGBColor[0, 0, 1]], StyleBox["\[IndentingNewLine]", FontColor->RGBColor[0, 0, 1]], StyleBox[\(m\ = \ \(-\((a[\([prc, ii]\)]\ /= \ mpiv)\)\)\), FontColor->RGBColor[0, 0, 1]], StyleBox[";", FontColor->RGBColor[0, 0, 1]], "\[IndentingNewLine]", StyleBox[\(row\ = \ a[\([iip, rc]\)]\), FontColor->RGBColor[0, 0, 1]], StyleBox[";", FontColor->RGBColor[0, 0, 1]], StyleBox["\[IndentingNewLine]", FontColor->RGBColor[0, 0, 1]], StyleBox[\(a[\([prc, rc]\)]\ += \ Map[\((row\ #)\) &, \ m]\), FontColor->RGBColor[0, 0, 1]]}], StyleBox[",", FontColor->RGBColor[0, 0, 1]], "\[IndentingNewLine]", \({ii, n}\)}], "]"}], ";", "\[IndentingNewLine]", \(LU[a, pivot]\)}]}], "]"}]}]], "Input"], Cell[BoxData[ \(Timing[lufactor4[machine]\ === \ luf3]\)], "Input"], Cell[BoxData[ \(Timing[lufactor4[big]\ === \ lufb]\)], "Input"] }, Closed]] }, Open ]] }, FrontEndVersion->"4.0 for Microsoft Windows", ScreenRectangle->{{0, 1024}, {0, 695}}, WindowSize->{1015, 668}, WindowMargins->{{0, Automatic}, {Automatic, 0}} ] (*********************************************************************** Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. ***********************************************************************) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[1739, 51, 55, 1, 105, "Title"], Cell[1797, 54, 161, 6, 65, "Subtitle"], Cell[1961, 62, 1064, 21, 109, "Text"], Cell[CellGroupData[{ Cell[3050, 87, 57, 0, 53, "Section"], Cell[3110, 89, 48, 0, 33, "Text"], Cell[CellGroupData[{ Cell[3183, 93, 56, 1, 30, "Input"], Cell[3242, 96, 46, 1, 29, "Output"] }, Open ]], Cell[3303, 100, 183, 4, 52, "Text"], Cell[CellGroupData[{ Cell[3511, 108, 47, 1, 30, "Input"], Cell[3561, 111, 192, 6, 43, "Output"] }, Open ]], Cell[3768, 120, 307, 7, 71, "Text"], Cell[4078, 129, 52, 0, 33, "Text"], Cell[CellGroupData[{ Cell[4155, 133, 54, 1, 30, "Input"], Cell[4212, 136, 46, 1, 29, "Output"] }, Open ]], Cell[4273, 140, 156, 3, 52, "Text"], Cell[CellGroupData[{ Cell[4454, 147, 52, 1, 30, "Input"], Cell[4509, 150, 38, 1, 29, "Output"] }, Open ]], Cell[4562, 154, 120, 3, 52, "Text"], Cell[CellGroupData[{ Cell[4707, 161, 79, 1, 30, "Input"], Cell[4789, 164, 47, 1, 29, "Output"] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell[4885, 171, 46, 0, 33, "Section"], Cell[4934, 173, 85, 2, 33, "Text"], Cell[5022, 177, 57, 0, 33, "Text"], Cell[CellGroupData[{ Cell[5104, 181, 66, 1, 30, "Input"], Cell[5173, 184, 38, 1, 29, "Output"] }, Open ]], Cell[5226, 188, 127, 3, 52, "Text"], Cell[CellGroupData[{ Cell[5378, 195, 80, 1, 30, "Input"], Cell[5461, 198, 39, 1, 29, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[5537, 204, 40, 1, 30, "Input"], Cell[5580, 207, 53, 1, 29, "Output"] }, Open ]], Cell[5648, 211, 73, 0, 33, "Text"], Cell[CellGroupData[{ Cell[5746, 215, 58, 1, 30, "Input"], Cell[5807, 218, 70, 1, 29, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[5914, 224, 58, 1, 30, "Input"], Cell[5975, 227, 68, 1, 29, "Output"] }, Open ]], Cell[6058, 231, 145, 3, 52, "Text"], Cell[6206, 236, 479, 9, 109, "Text"], Cell[6688, 247, 512, 11, 185, "Text"], Cell[7203, 260, 164, 3, 52, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[7404, 268, 49, 0, 33, "Section"], Cell[7456, 270, 178, 3, 52, "Text"], Cell[CellGroupData[{ Cell[7659, 277, 139, 3, 60, "Subsubsection"], Cell[7801, 282, 90, 2, 33, "Text"], Cell[7894, 286, 109, 2, 50, "Input"], Cell[8006, 290, 104, 3, 33, "Text"], Cell[8113, 295, 145, 3, 52, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[8295, 303, 101, 2, 29, "Subsubsection"], Cell[8399, 307, 50, 0, 33, "Text"], Cell[CellGroupData[{ Cell[8474, 311, 86, 1, 30, "Input"], Cell[8563, 314, 38, 1, 29, "Output"] }, Open ]], Cell[8616, 318, 224, 4, 71, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[8877, 327, 134, 3, 46, "Subsubsection"], Cell[9014, 332, 203, 4, 52, "Text"], Cell[9220, 338, 34, 0, 33, "Text"], Cell[CellGroupData[{ Cell[9279, 342, 121, 3, 50, "Input"], Cell[9403, 347, 71, 1, 29, "Output"] }, Open ]], Cell[9489, 351, 285, 5, 71, "Text"] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell[9823, 362, 52, 0, 33, "Section"], Cell[9878, 364, 40, 0, 33, "Text"], Cell[9921, 366, 65, 1, 30, "Input"], Cell[9989, 369, 77, 1, 30, "Input"], Cell[10069, 372, 89, 1, 30, "Input"], Cell[10161, 375, 59, 1, 30, "Input"], Cell[10223, 378, 90, 2, 33, "Text"], Cell[10316, 382, 156, 4, 50, "Input"], Cell[10475, 388, 456, 7, 109, "Text"], Cell[10934, 397, 237, 5, 90, "Input"], Cell[11174, 404, 44, 0, 33, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[11255, 409, 78, 0, 33, "Section"], Cell[11336, 411, 160, 3, 52, "Text"], Cell[11499, 416, 306, 9, 71, "Text"], Cell[11808, 427, 149, 3, 52, "Text"], Cell[11960, 432, 102, 3, 33, "Text"], Cell[CellGroupData[{ Cell[12087, 439, 71, 1, 30, "Input"], Cell[12161, 442, 109, 2, 29, "Output"] }, Open ]], Cell[12285, 447, 109, 2, 30, "Input"], Cell[CellGroupData[{ Cell[12419, 453, 70, 1, 30, "Input"], Cell[12492, 456, 131, 2, 42, "Message"], Cell[12626, 460, 74, 1, 29, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[12737, 466, 59, 1, 30, "Input"], Cell[12799, 469, 57, 1, 29, "Output"] }, Open ]], Cell[12871, 473, 128, 3, 52, "Text"], Cell[13002, 478, 135, 2, 50, "Input"], Cell[13140, 482, 113, 2, 30, "Input"], Cell[CellGroupData[{ Cell[13278, 488, 50, 1, 30, "Input"], Cell[13331, 491, 85, 1, 29, "Output"] }, Open ]], Cell[13431, 495, 109, 2, 50, "Input"], Cell[CellGroupData[{ Cell[13565, 501, 50, 1, 30, "Input"], Cell[13618, 504, 67, 1, 29, "Output"] }, Open ]], Cell[13700, 508, 72, 0, 33, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[13809, 513, 44, 0, 33, "Section"], Cell[13856, 515, 619, 11, 71, "Text"], Cell[14478, 528, 4984, 100, 370, "Input"], Cell[19465, 630, 106, 3, 33, "Text"], Cell[19574, 635, 237, 4, 70, "Input"], Cell[19814, 641, 73, 1, 30, "Input"], Cell[19890, 644, 70, 1, 30, "Input"], Cell[19963, 647, 70, 1, 30, "Input"], Cell[20036, 650, 66, 1, 30, "Input"], Cell[20105, 653, 283, 6, 52, "Text"], Cell[20391, 661, 398, 8, 52, "Text"], Cell[20792, 671, 4849, 100, 390, "Input"], Cell[25644, 773, 71, 1, 30, "Input"], Cell[25718, 776, 68, 1, 30, "Input"], Cell[25789, 779, 272, 5, 52, "Text"], Cell[26064, 786, 75, 1, 30, "Input"], Cell[26142, 789, 310, 6, 50, "Input"], Cell[26455, 797, 160, 3, 33, "Text"], Cell[26618, 802, 55, 0, 33, "Text"], Cell[26676, 804, 89, 1, 30, "Input"], Cell[26768, 807, 53, 0, 33, "Text"], Cell[26824, 809, 1833, 36, 250, "Input"], Cell[28660, 847, 71, 1, 30, "Input"], Cell[28734, 850, 68, 1, 30, "Input"], Cell[28805, 853, 326, 7, 52, "Text"], Cell[29134, 862, 2472, 47, 270, "Input"], Cell[31609, 911, 71, 1, 30, "Input"], Cell[31683, 914, 186, 3, 50, "Input"], Cell[31872, 919, 68, 1, 30, "Input"], Cell[31943, 922, 157, 3, 33, "Text"], Cell[32103, 927, 2140, 42, 270, "Input"], Cell[34246, 971, 72, 1, 30, "Input"], Cell[34321, 974, 68, 1, 30, "Input"] }, Closed]] }, Open ]] } ] *) (*********************************************************************** End of Mathematica Notebook file. ***********************************************************************)