BLAS & LAPACK
Performance
in
Mathematica 5.0

Schoeller Porter
Software Technology
Wolfram Research, Inc.

Needs["Graphics`MultipleListPlot`"];
Needs["Graphics`Legend`"];
Needs["Graphics`"];
SetDirectory[$HomeDirectory  <> "/Developers Conference \
2003"];

SetOptions[MultipleListPlot, PlotRange -> Automatic, 
        PlotJoined -> True, SymbolShape -> None ,  
        Frame -> 
          True, GridLines  ->  {{Automatic, #}, {Automatic, #}} &[
          GrayLevel[0.5]], 
        PlotStyle -> 
          Map[{Thickness[0.004], #} &, {RGBColor[1, 0.4, 0.4], 
              RGBColor[0.4, 0.4, 1], RGBColor[0.4, 0.8, 0.4],  
              RGBColor[0.8, 0.4, 0.8]}], ImageSize -> 850,  
        RotateLabel -> False, LegendSize -> {1.4, 0.14}, 
        LegendOrientation -> Horizontal, 
        LegendPosition -> {-0.58, -0.8}];
$TextStyle = {FontSize -> 16,  FontFamily -> "Helvetica",  
          FontWeight -> "Bold"};

Get["DC03Data.mx"];

Introduction

Advantages

Significant performance improvement

Support for multi-processor platforms

A weath of vendor-supported and free, high performance libraries

Topics

Dot[]

LinearSolve[]

Eigenvalues[]

SingularValueDecomposition[]

A Comparison with MATLAB

BLAS & LAPACK Integration with Mathematica

Windows, Linux IA32, and Linux Itanium

Intel Math Kernel Library 5.2sp1

• Equivalent performance on Windows and Linux on IA32
• SSE2-optimizations for the Pentium 4 and Itanium processors
• Automatic, runtime processor detection
• Multi-threaded BLAS

Mac OS X

Apple vecLib framework

• Included with Mac OS 10.2
• AltiVec optimizations for the G4 processor
• Automatic, runtime processor detection
• Multi-threaded BLAS

BLAS & LAPACK Integration with Mathematica (Continued . . . )

AIX

ATLAS 3.4.0

HP-UX PA-RISC

HP MLIB 8.1

Linux Alpha

ATLAS 3.4.0

Solaris

ATLAS 3.4.0

Tru64

ATLAS 3.4.0

A Close Look at Dot[ ] . . .

Description:

Dot[a,b,c] gives the products of vectors, matrices, and tensors.

BLAS/LAPACK Implementation:

[d,z]gemm()

An  Example in Mathematica 4.2.1:

SeedRandom[1];
matrix = 
        Table[Random[Real], {500}, {500}];
{Timing[Dot[matrix, matrix];][[1]], $Version}

An Example in Mathematica 5.0 Beta 2

SeedRandom[1];
matrix = 
        Table[Random[Real], {500}, {500}];
{Timing[Dot[matrix, matrix];][[1]], $Version}

. . . Dot[ ] on Windows

MultipleListPlot[
            Table[Map[{#[[1]], 
                    10^-6*2*#[[1]]^3/#[[2]]} &, #[[i, 
                    2]]], {i, 2, 5}],  
            FrameLabel -> {"Vector   Length",  "MFLOPS"},  
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 2, 5}],  
            PlotLabel ->  "MFLOPS  on a  2.4 GHz Pentium 4"] &[
        DotData["Windows"]];

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 2, 5}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 2, 5}],  
            PlotLabel -> "Evaluation Time on a 2.4 GHz Pentium 4"] &[
        DotData["Windows"]];

. . . Dot[ ] on Mac OS X

MultipleListPlot[
            Table[Map[{#[[1]], 
                    10^-6*2*#[[1]]^3/#[[2]]} &, #[[i, 
                    2]]], {i, 1, 3}],  
            FrameLabel -> {"Vector   Length",  "MFLOPS"},  
            PlotLegend -> Transpose[#][[1]],  
            PlotLabel -> "MFLOPS  on a 866 MHz G4"] &[
        DotData["Darwin"]];

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 1, 3}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> Transpose[#][[1]],  
            PlotLabel -> "Evaluation Time on a 866 MHz G4"] &[
        DotData["Darwin"]];

. . . Dot[ ] on Linux IA32

MultipleListPlot[
            Table[Map[{#[[1]], 
                    10^-6*2*#[[1]]^3/#[[2]]} &, #[[i, 
                    2]]], {i, 1, 4}], 
            FrameLabel -> {"Vector   Length", "MFLOPS"}, 
            PlotLegend -> Transpose[#][[1]],  
            PlotLabel ->  "MFLOPS  on a 850 MHz Athlon"] &[
        DotData["Linux"]];

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 1, 4}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> Transpose[#][[1]],  
            PlotLabel -> "Evaluation Time on a 850 MHz Athlon"] &[
        DotData["Linux"]];

. . . Dot[ ] on Linux Itanium

MultipleListPlot[
            Table[Map[{#[[1]], 
                    10^-6*2*#[[1]]^3/#[[2]]} &, #[[i, 
                    2]]], {i, 1, 4}],  
            FrameLabel -> {"Vector   Length",  "MFLOPS"},  
            PlotLegend -> Transpose[#][[1]],  
            PlotLabel -> "MFLOPS  on a 900 MHz 4-way  Itanium 2"] &[
        DotData["Linux-IA64"]];

\(\(\(MultipleListPlot[
            Table[Map[{#[\([1]\)], Log[10, #[\([2]\)]]} &, #[\([i, 
                    2]\)]], {i, 1, 4}], \ 
            FrameLabel -> {"Vector   Length", \ "Time (s)"}, \ 
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> \(Transpose[#]\)[\([1]\)], 
            PlotLabel -> "Evaluation Time on a 900 MHz 4-way  Itanium 2"] &\)[DotData["Linux-IA64"]];\)\)

A Close Look at LinearSolve[ ] . . .

Description:

LinearSolve[m,b] finds a solution, x, to the matrix equation  [m]·x= b.

BLAS/LAPACK Implementation:

[D,Z]GESVX( )

An Example in Mathematica 4.2.1

SeedRandom[1];
matrix = 
        Table[Random[Real], {500}, {500}];
vector  =  Table[Random[Real], {500}];
{Timing[LinearSolve[matrix, 
              vector];][[1]], $Version}

An Example in Mathematica 5.0 Beta 2

SeedRandom[1];
matrix = 
        Table[Random[Real], {500}, {500}];
vector  =  Table[Random[Real], {500}];
{Timing[LinearSolve[matrix, 
              vector];][[1]], $Version}

. . . LinearSolve[ ] on Windows

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 2, 5}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 2, 5}],  
            PlotLabel -> "Evaluation Time on a 2.4 GHz Pentium 4"] &[
        LinearSolveData["Windows"]];

. . . LinearSolve[ ] on Mac OS X

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 1, 3}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 1, 3}],  
            PlotLabel -> "Evaluation Time on a 866 MHz G4"] &[
        LinearSolveData["Darwin"]];

. . . LinearSolve[ ] on Linux IA32

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 1, 4}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 1, 4}],  
            PlotLabel -> "Evaluation Timne on a 850 MHz Athlon"] &[
        LinearSolveData["Linux"]];

A Close Look at Eigensystem[ ] . . .

Description:

Eigensystem[m] gives a list of eigenvalues and eigenvectors of the square matrix m.

BLAS/LAPACK Implementation:

[D,Z]GEEV( )

An Example in Mathematica 4.2

SeedRandom[1];
a = Table[Random[Real], {350}, {350}];
b = 
        DiagonalMatrix[
          Table[Random[Real], {350}]];
matrix = a  .  b  .  Inverse[a];
{Timing[Eigensystem[matrix];][[1]], $Version}

An Example in Mathematica 5.0 Beta 2

SeedRandom[1];
a = Table[Random[Real], {350}, {350}];
b = 
        DiagonalMatrix[
          Table[Random[Real], {350}]];
matrix = a  .  b  .  Inverse[a];
{Timing[Eigensystem[matrix];][[1]], $Version}

. . . Eigensystem[ ] on Windows

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 2, 5}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 2, 5}],  
            PlotLabel -> "Evaluation Time on a 2.4 GHz Pentium 4"] &[
        EigensystemData["Windows"]];

. . . Eigensystem[ ] on Mac OS X

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 1, 3}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 1, 3}],  
            PlotLabel -> "Evaluation Time on a 866 MHz G4"] &[
        EigensystemData["Darwin"]];

. . . Eigensystem[ ] on Linux IA32

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 1, 4}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 1, 4}],  
            PlotLabel -> "Evaluation Time on a 850 MHz Athlon"] &[
        EigensystemData["Linux"]];

A Close Look at SingularValueDecomposition[ ] . . .

Description:

SingularValueDecomposition[m] gives the singular value decomposition for a matrix m, resulting in a list of matrices {u,v,w} where w is diagonal and m can be written as Conjugate[Transpose[u]] · w · v.

BLAS/LAPACK Implementation:

[D,Z]GESDD( )

An Example in Mathematica 4.2.1

SeedRandom[1];
matrix = 
        Table[Random[Real], {300}, {300}];
{Timing[SingularValues[
              matrix];][[1]], $Version}

An Example in Mathematica 5.0 Beta 2

SeedRandom[1];
matrix = 
        Table[Random[Real], {300}, {300}];
{Timing[SingularValueDecomposition[
              matrix];][[1]], $Version}

. . . SingularValueDecomposition[ ] on Windows

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 2, 5}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 2, 5}],  
            PlotLabel -> "Evaluation Time on a 2.4 GHz Pentium 4"] &[
        SVDData["Windows"]];

. . . SingularValueDecomposition[ ] on Mac OS X

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 1, 3}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 1, 3}],  
            PlotLabel -> "Evaluation Time on a 866 MHz G4"] &[
        SVDData["Darwin"]];

. . . SingularValueDecomposition[ ] on Linux IA32

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 1, 4}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 1, 4}],  
            PlotLabel -> "Evaluation Time on a 850 MHz Athlon"] &[
        SVDData["Linux"]];

A Comparison with MATLAB 6.5r13 . . .

Summary

The performance of Dot[] is nearly equivalent to  MATLAB's "*" function with the Pentium Pro library. 
With the Pentium 4 library Mathematica 5.0 is nearly 2 times faster.

LinearSolve[] starts 2 times slower than MATLAB's mldivide() function, but tends toward faster
 performance for large matrix sizes, eventually performing faster than MATLAB with the Pentium 4 library.

Eigensystem[] is within 20-30% of MATLAB's eig() function, for all matrix sizes.

SingularValueDecomposition[] is consistently 3 times faster than MATLAB's svd() function with the
Pentium Pro library, and 3 - 4 times faster with the Pentium 4 library.

. . . Dot[ ] vs *

MultipleListPlot[
            Table[Map[{#[[1]], 
                    10^-6*2*#[[1]]^3/#[[2]]} &, #[[i, 
                    2]]], {i,  1,  3}],  
            FrameLabel -> {"Vector   Length",  "MFLOPS"},  
            PlotLegend -> Transpose[#][[1]],  
            PlotLabel -> "MFLOPS  on a  2.4 GHz Pentium 4"] &[
        DotData["Windows"]];

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 1, 3}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> Transpose[#][[1]],  
            PlotLabel  ->  "Evaluation Time on a 2.4 GHz Pentium \
4"] &[DotData["Windows"]];

. . . LinearSolve[ ] vs mldivide( )

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 1, 3}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 1, 3}], 
            PlotLabel  ->  "Evaluation Time on a 2.4 GHz Pentium \
4"] &[LinearSolveData["Windows"]];

. . . Eigensystem[ ] vs eig( )

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 1, 3}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 1, 3}], 
            PlotLabel  ->  "Evaluation Time on a 2.4 GHz Pentium \
4"] &[EigensystemData["Windows"]];

. . . SingularValueDecomposition[ ] vs svd( )

MultipleListPlot[
            Table[Map[{#[[1]], Log[10, #[[2]]]} &, #[[i, 
                    2]]], {i, 1, 3}],  
            FrameLabel -> {"Vector   Length",  "Time (s)"},  
            FrameTicks -> {Automatic, LogScale}, 
            PlotLegend -> 
              Table[Transpose[#][[1, i]], {i, 1, 3}], 
            PlotLabel  ->  "Evaluation Time on a 2.4 GHz Pentium \
4"] &[SVDData["Windows"]];