(*********************************************************************** 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[ 36422, 1098]*) (*NotebookOutlinePosition[ 37293, 1129]*) (* CellTagsIndexPosition[ 37209, 1123]*) (*WindowFrame->Normal*) Notebook[{ Cell["Packed Arrays", "Title"], Cell[TextData[{ "An efficient internal storage format in ", StyleBox["Mathematica", FontSlant->"Italic"], " 4" }], "Subtitle"], Cell["\<\ Rob Knapp Numerical Algorithms Development Wolfram Research, Inc. rknapp@wolfram.com\ \>", "Subsubtitle"], 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, and is visible in the ", StyleBox["List", "MR"], " operations that were added to the compiler. Basically what has been 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 ", StyleBox["Developer`", "MR"], " 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." }], "Text"], Cell[CellGroupData[{ Cell[TextData[{ "Representations of data in ", StyleBox["Mathematica", FontSlant->"Italic"] }], "Section"], Cell[TextData[{ "If you enter a list in ", StyleBox["Mathematica", FontSlant->"Italic"], "," }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(list\ = \ {1, 2, 3, 4}\)], "Input"], Cell[BoxData[ \({1, 2, 3, 4}\)], "Output"] }, Closed]], Cell[TextData[{ "its internal representation may be different from the ", StyleBox["OutputForm", "MR"], ".", " ", "The basic common representation for lists (and all expressions) is given \ by the ", StyleBox["FullForm", "MR"], "." }], "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"] }, Closed]], 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 expression \ (in this case, list, or number which is a machine integer)." }], "Text"], Cell["If you enter the data a different way,", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(packed\ = \ Range[4]\)], "Input"], Cell[BoxData[ \({1, 2, 3, 4}\)], "Output"] }, Closed]], Cell["\<\ it is internally stored in a different way. The basic print forms are all the \ same. For all intents and purposes it is the same.\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(packed\ === \ list\)], "Input"], Cell[BoxData[ \(True\)], "Output"] }, Closed]], Cell[TextData[{ "You can only tell that it is stored in a different way by using a ", StyleBox["Developer`", "MR"], " context function:" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Developer`PackedArrayQ\ /@ \ {list, \ packed}\)], "Input"], Cell[BoxData[ \({False, True}\)], "Output"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["Advantages of packed arrays", "Section"], Cell["\<\ There are two advantages of packed arrays, memory and speed. \ \>", "Text", CellTags->"Advantages of PAs"], Cell[" The memory is the most obvious of these.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(ByteCount[a = N[Range[1000]]]\)], "Input"], Cell[BoxData[ \(TraditionalForm\`8056\)], "Output"] }, Closed]], Cell["\<\ If you 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\[LeftDoubleBracket]1\[RightDoubleBracket] = 1; b]\)], "Input"], Cell[BoxData[ \(TraditionalForm\`20024\)], "Output"] }, Closed]], Cell[CellGroupData[{ Cell[BoxData[ \(N[%\/%%]\)], "Input"], Cell[BoxData[ \(TraditionalForm\`2.4856007944389273`\)], "Output"] }, Closed]], Cell["The speed advantage typically comes in common operations.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[Do[a + a, {1000}]]\)], "Input"], Cell[BoxData[ \(TraditionalForm\`{0.6099999999999999`\ Second, Null}\)], "Output"] }, Closed]], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[Do[b + b, {1000}]]\)], "Input"], Cell[BoxData[ \(TraditionalForm\`{8.24`\ Second, Null}\)], "Output"] }, Closed]], 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, ", StyleBox["a", "MR"], " is stored as a packed array, so when ", StyleBox["Plus", "MR"], " sees it, it knows all of the elements are machine real numbers, so that \ internally, a function which adds arrays of real numbers is used. Except for \ overflow and underflow checking, which ", StyleBox["Mathematica", FontSlant->"Italic"], " does, this is done as fast as you could get by writing a loop to do the \ addition in a C-program for two arrays of machine numbers." }], "Text"], Cell[TextData[{ "In the second operation, ", StyleBox["b", "MR"], " is stored as ", StyleBox["List[\[Ellipsis]]", "MR"], ", so what happens first is ", StyleBox["Plus", "MR"], " is threaded over the lists, producing (roughly) the following." }], "Text"], Cell[BoxData[ \({b[\([1]\)] + b[\([1]\)], b[\([2]\)] + b[\([2]\)], \(\(...\)\(.\)\)}\)], "Input", Evaluatable->False], Cell[TextData[{ StyleBox["Plus", "MR"], " processes 1000 expressions. For each pair of numbers, ", StyleBox["Plus", "MR"], " uses the information attached to those numbers (", StyleBox["i.e.", FontSlant->"Italic"], ", real, integer, machine, bignum, and so forth) 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 increase in speed.\ \>", "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 situations that 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[TextData[{ "Listable operations mean those with the attribute ", StyleBox["Listable", "MR"], ". " }], "Text"], Cell[BoxData[ \(Select[Names["\"], \ MemberQ[Attributes[#], \ Listable] &]\)], "Input"], Cell["\<\ The majority of which will produce packed arrays given packed array \ arguments.\ \>", "Text"], Cell[TextData[{ "Structural operations mean commands like ", StyleBox["Part", "MR"], ", ", StyleBox["Join", "MR"], ", ", StyleBox["Transpose", "MR"], ", ", "and so forth", ", which act on the structure of an expression." }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["\<\ 2. Functions that use an internal array format for machine numbers. \ \>", "Subsubsection"], Cell[TextData[{ StyleBox["Fourier", "MR"], " is a good example of this:" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Developer`PackedArrayQ[Fourier[{1. , 2. , 3. , 4. }]]\)], "Input"], Cell[BoxData[ \(True\)], "Output"] }, Closed]], Cell[TextData[{ "Whether the input is a packed array or not, if ", StyleBox["Fourier", "MR"], " 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[TextData[{ "You have already used one example in this class, ", StyleBox["Range", "MR"], ". In ", StyleBox["Range[4]", "MR"], ", it is easy to see before generating any numbers that the result can be \ represented as a packed array." }], "Text"], Cell["\<\ For something like the following, it is a little more difficult to determine.\ \ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[ Developer`PackedArrayQ[ data\ = \ Table[Random[], {10000}, {2}]]]\)], "Input"], Cell[BoxData[ \({0.050000000000000044`\ Second, True}\)], "Output"] }, Closed]], Cell["\<\ Commands like this are handled by auto-compilation. However, since trying the \ compiler does take some time, there is a (user modifiable) threshold length \ that must be exceeded before compilation is attempted.\ \>", "Text"] }, Closed]] }, 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[CellGroupData[{ Cell[BoxData[ \(Timing[Developer`PackedArrayQ[Map[f, data]]]\)], "Input"], Cell[OutputFormData["\<\ {0.*Second, False}\ \>", "\<\ {0. Second, False}\ \>"], "Output"] }, Closed]], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[\(udata\ = \ Developer`FromPackedArray[data];\)]\)], "Input"], Cell[OutputFormData["\<\ {0.*Second, Null}\ \>", "\<\ {0. Second, Null}\ \>"], "Output"] }, Closed]], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[\(Map[f, udata];\)]\)], "Input"], Cell[OutputFormData["\<\ {0.*Second, Null}\ \>", "\<\ {0. Second, Null}\ \>"], "Output"] }, Closed]], Cell["\<\ However, if you inline the function, the auto-compilation can work.\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[ Developer`PackedArrayQ[ Map[\((\((#[\([1]\)])\)^2\ + \ \((#[\([2]\)])\)^2)\) &, \ data]]]\)], "Input"], Cell[OutputFormData["\<\ {0.*Second, False}\ \>", "\<\ {0. Second, False}\ \>"], "Output"] }, Closed]], Cell["\<\ 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 cases, please take into account the whole time going into your \ computation. In the above example, you were not taking into account the time \ it took to generate the random list. Without auto-compilation and packed \ arrays, it takes quite a bit longer to compute.\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[{ \(\(Developer`SetSystemOptions["\" -> Infinity];\)\), "\[IndentingNewLine]", \(Timing[ Developer`PackedArrayQ[ data\ = \ Table[Random[], {10000}, {2}]]]\)}], "Input"], Cell[OutputFormData["\<\ {0.16999999999984539*Second, False}\ \>", "\<\ {0.17 Second, False}\ \>"], "Output"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["How do I write code that 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["\<\ Some aids are provided that 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"] }, Closed]], 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"] }, Closed]], Cell[CellGroupData[{ Cell[BoxData[ \(Timing[Tr[N[Range[1000]]]]\)], "Input"], Cell[BoxData[ \({0.`\ Second, 500500.`}\)], "Output"] }, Closed]], Cell["\<\ Sometimes, you know your code will run more efficiently if you pack the input \ beforehand. \ \>", "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"] }, Closed]], 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"] }, Closed]], Cell["Probably, the best thing to 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 ", StyleBox["LinearAlgebra`GaussianElimination", "MR"], ". Since this functionality exists in the kernel, the real purpose of this \ package is pedagogical, but the original author seems to have 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 a prime example of \ how you should not write code in ", StyleBox["Mathematica", FontSlant->"Italic"], ". Let\[CloseCurlyQuote]s see why and how you 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 you 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! You are a factor of 300 slower in the machine case and a factor of 5 \ slower in the bignum example. Before screaming an exasperated \ \[OpenCurlyDoubleQuote]", StyleBox["Mathematica", FontSlant->"Italic"], " is just too slow,\[CloseCurlyDoubleQuote] lets try to improve the code." }], "Text"], Cell[TextData[{ "The first thing noticed is that this code is using ", StyleBox["For", "MR"], " to handle the loops. In ", StyleBox["Mathematica", FontSlant->"Italic"], ", ", StyleBox["For", "MR"], " is slower than the alternative iterator constructs. Try to avoid using ", StyleBox["For", "MR"], " in innermost loops unless it would make for more complicated code. In \ this example simply changing from ", StyleBox["For", "MR"], " to ", StyleBox["Do", "MR"], " 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 will help make the code more readable. At this point, it won\ \[CloseCurlyQuote]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[TextData[{ "This function finds which of the possible pivots (in column ", StyleBox["c", "MR"], ", rows ", StyleBox["c", "MR"], " and above) is largest in absolute value, and swaps with the one in row ", StyleBox["c", "MR"], "." }], "Text"], Cell["This, in turn, uses the simple function", "Text"], Cell[BoxData[ \(\(MaxPosition[x_]\ := \ First[Position[x, \ Max[x]]];\)\)], "Input"], Cell["and puts 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[{ "You are still going in the right direction, but there is no big savings. \ The place to look for big savings is the innermost loop. The next change \ replaces the innermost loop by giving ", StyleBox["Mathematica", FontSlant->"Italic"], " commands that 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 successful. There is still \ a set of nested loops that 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]] }, FrontEndVersion->"4.0 for Microsoft Windows", ScreenRectangle->{{0, 1024}, {0, 695}}, ScreenStyleEnvironment->"Presentation", WindowSize->{1016, 668}, WindowMargins->{{0, Automatic}, {Automatic, 0}}, StyleDefinitions -> "Demo.nb" ] (*********************************************************************** 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->{ "Advantages of PAs"->{ Cell[5176, 193, 118, 3, 45, "Text", CellTags->"Advantages of PAs"]} } *) (*CellTagsIndex CellTagsIndex->{ {"Advantages of PAs", 37105, 1116} } *) (*NotebookFileOutline Notebook[{ Cell[1717, 49, 30, 0, 111, "Title"], Cell[1750, 51, 136, 5, 81, "Subtitle"], Cell[1889, 58, 115, 5, 152, "Subsubtitle"], Cell[2007, 65, 1033, 23, 170, "Text"], Cell[CellGroupData[{ Cell[3065, 92, 114, 4, 77, "Section"], Cell[3182, 98, 114, 5, 45, "Text"], Cell[CellGroupData[{ Cell[3321, 107, 56, 1, 54, "Input"], Cell[3380, 110, 46, 1, 54, "Output"] }, Closed]], Cell[3441, 114, 257, 9, 64, "Text"], Cell[CellGroupData[{ Cell[3723, 127, 47, 1, 54, "Input"], Cell[3773, 130, 192, 6, 54, "Output"] }, Closed]], Cell[3980, 139, 303, 7, 64, "Text"], Cell[4286, 148, 54, 0, 45, "Text"], Cell[CellGroupData[{ Cell[4365, 152, 54, 1, 54, "Input"], Cell[4422, 155, 46, 1, 38, "Output"] }, Closed]], Cell[4483, 159, 154, 3, 39, "Text"], Cell[CellGroupData[{ Cell[4662, 166, 52, 1, 54, "Input"], Cell[4717, 169, 38, 1, 38, "Output"] }, Closed]], Cell[4770, 173, 154, 4, 39, "Text"], Cell[CellGroupData[{ Cell[4949, 181, 79, 1, 54, "Input"], Cell[5031, 184, 47, 1, 54, "Output"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[5127, 191, 46, 0, 49, "Section"], Cell[5176, 193, 118, 3, 45, "Text", CellTags->"Advantages of PAs"], Cell[5297, 198, 57, 0, 45, "Text"], Cell[CellGroupData[{ Cell[5379, 202, 62, 1, 54, "Input"], Cell[5444, 205, 55, 1, 55, "Output"] }, Closed]], Cell[5514, 209, 128, 3, 39, "Text"], Cell[CellGroupData[{ Cell[5667, 216, 108, 2, 54, "Input"], Cell[5778, 220, 56, 1, 55, "Output"] }, Closed]], Cell[CellGroupData[{ Cell[5871, 226, 41, 1, 64, "Input"], Cell[5915, 229, 70, 1, 55, "Output"] }, Closed]], Cell[6000, 233, 73, 0, 39, "Text"], Cell[CellGroupData[{ Cell[6098, 237, 58, 1, 54, "Input"], Cell[6159, 240, 86, 1, 55, "Output"] }, Closed]], Cell[CellGroupData[{ Cell[6282, 246, 58, 1, 48, "Input"], Cell[6343, 249, 72, 1, 55, "Output"] }, Closed]], Cell[6430, 253, 145, 3, 39, "Text"], Cell[6578, 258, 539, 12, 95, "Text"], Cell[7120, 272, 267, 8, 45, "Text"], Cell[7390, 282, 130, 3, 54, "Input", Evaluatable->False], Cell[7523, 287, 427, 10, 70, "Text"], Cell[7953, 299, 174, 3, 45, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[8164, 307, 49, 0, 49, "Section"], Cell[8216, 309, 175, 3, 45, "Text"], Cell[CellGroupData[{ Cell[8416, 316, 138, 4, 97, "Subsubsection"], Cell[8557, 322, 120, 4, 45, "Text"], Cell[8680, 328, 109, 2, 54, "Input"], Cell[8792, 332, 104, 3, 45, "Text"], Cell[8899, 337, 250, 10, 45, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[9186, 352, 101, 3, 55, "Subsubsection"], Cell[9290, 357, 88, 3, 45, "Text"], Cell[CellGroupData[{ Cell[9403, 364, 86, 1, 54, "Input"], Cell[9492, 367, 38, 1, 70, "Output"] }, Closed]], Cell[9545, 371, 258, 5, 64, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[9840, 381, 134, 4, 75, "Subsubsection"], Cell[9977, 387, 260, 7, 70, "Text"], Cell[10240, 396, 103, 3, 45, "Text"], Cell[CellGroupData[{ Cell[10368, 403, 121, 3, 54, "Input"], Cell[10492, 408, 71, 1, 54, "Output"] }, Closed]], Cell[10578, 412, 236, 4, 64, "Text"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[10863, 422, 52, 0, 49, "Section"], Cell[10918, 424, 40, 0, 45, "Text"], Cell[10961, 426, 65, 1, 54, "Input"], Cell[CellGroupData[{ Cell[11051, 431, 77, 1, 54, "Input"], Cell[11131, 434, 90, 4, 36, "Output"] }, Closed]], Cell[CellGroupData[{ Cell[11258, 443, 89, 1, 48, "Input"], Cell[11350, 446, 88, 4, 36, "Output"] }, Closed]], Cell[CellGroupData[{ Cell[11475, 455, 59, 1, 48, "Input"], Cell[11537, 458, 88, 4, 36, "Output"] }, Closed]], Cell[11640, 465, 91, 2, 39, "Text"], Cell[CellGroupData[{ Cell[11756, 471, 156, 4, 54, "Input"], Cell[11915, 477, 90, 4, 36, "Output"] }, Closed]], Cell[12020, 484, 460, 7, 89, "Text"], Cell[CellGroupData[{ Cell[12505, 495, 237, 5, 77, "Input"], Cell[12745, 502, 109, 4, 36, "Output"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[12903, 512, 77, 0, 49, "Section"], Cell[12983, 514, 160, 3, 45, "Text"], Cell[13146, 519, 305, 9, 70, "Text"], Cell[13454, 530, 149, 3, 45, "Text"], Cell[13606, 535, 95, 2, 45, "Text"], Cell[CellGroupData[{ Cell[13726, 541, 71, 1, 54, "Input"], Cell[13800, 544, 109, 2, 38, "Output"] }, Closed]], Cell[13924, 549, 109, 2, 48, "Input"], Cell[CellGroupData[{ Cell[14058, 555, 70, 1, 54, "Input"], Cell[14131, 558, 131, 2, 28, "Message"], Cell[14265, 562, 74, 1, 38, "Output"] }, Closed]], Cell[CellGroupData[{ Cell[14376, 568, 59, 1, 48, "Input"], Cell[14438, 571, 57, 1, 38, "Output"] }, Closed]], Cell[14510, 575, 116, 3, 39, "Text"], Cell[14629, 580, 135, 2, 77, "Input"], Cell[14767, 584, 113, 2, 54, "Input"], Cell[CellGroupData[{ Cell[14905, 590, 50, 1, 54, "Input"], Cell[14958, 593, 85, 1, 38, "Output"] }, Closed]], Cell[15058, 597, 109, 2, 48, "Input"], Cell[CellGroupData[{ Cell[15192, 603, 50, 1, 54, "Input"], Cell[15245, 606, 67, 1, 38, "Output"] }, Closed]], Cell[15327, 610, 70, 0, 39, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[15434, 615, 44, 0, 49, "Section"], Cell[15481, 617, 646, 12, 120, "Text"], Cell[16130, 631, 5014, 101, 445, "Input"], Cell[21147, 734, 102, 3, 45, "Text"], Cell[21252, 739, 237, 4, 100, "Input"], Cell[21492, 745, 73, 1, 54, "Input"], Cell[21568, 748, 70, 1, 54, "Input"], Cell[21641, 751, 70, 1, 54, "Input"], Cell[21714, 754, 66, 1, 54, "Input"], Cell[21783, 757, 326, 7, 70, "Text"], Cell[22112, 766, 538, 16, 95, "Text"], Cell[22653, 784, 4902, 102, 468, "Input"], Cell[27558, 888, 71, 1, 54, "Input"], Cell[27632, 891, 68, 1, 54, "Input"], Cell[27703, 894, 275, 5, 70, "Text"], Cell[27981, 901, 75, 1, 54, "Input"], Cell[28059, 904, 310, 6, 100, "Input"], Cell[28372, 912, 259, 8, 45, "Text"], Cell[28634, 922, 55, 0, 45, "Text"], Cell[28692, 924, 89, 1, 54, "Input"], Cell[28784, 927, 51, 0, 45, "Text"], Cell[28838, 929, 1833, 36, 307, "Input"], Cell[30674, 967, 71, 1, 54, "Input"], Cell[30748, 970, 68, 1, 54, "Input"], Cell[30819, 973, 331, 7, 70, "Text"], Cell[31153, 982, 2472, 47, 330, "Input"], Cell[33628, 1031, 71, 1, 54, "Input"], Cell[33702, 1034, 186, 3, 77, "Input"], Cell[33891, 1039, 68, 1, 54, "Input"], Cell[33962, 1042, 155, 3, 45, "Text"], Cell[34120, 1047, 2140, 42, 330, "Input"], Cell[36263, 1091, 72, 1, 54, "Input"], Cell[36338, 1094, 68, 1, 54, "Input"] }, Closed]] } ] *) (*********************************************************************** End of Mathematica Notebook file. ***********************************************************************)