The Gauss application example

The Gauss application example uses a text-book Gauss-elimination algorithm to solve a linear equation system. Two versions of the Gauss elimination algorithm are presented:

Both versions of the algorithm are compiled to C++ code. The performance is measured, both including and without MathLink overhead.

Initialization

Load the MathCode package

First put a Needs statement, to make sure that the MathCode application is loaded:

Needs["MathCode`"]

Check current directory

Check the current directory, since a number of files will be placed there during the code generation process.

Directory[]

You might want to set the directory to some place where all generated files are put, e.g. the one below, or some other directory.

SetDirectory[$MCRoot<>$PathnameSeparator<>"Demos"<>$PathnameSeparator<>"Gauss"];

Check whether this file is in the current directory.

ExistFileQ["Gauss.nb"];

Off[General::spell1]

Start the Gauss Package

BeginPackage["Gauss`",{MathCodeContexts}];

Define exported symbols

Begin["Gauss`"];

GaussSolveMatlab;
GaussSolveForLoops;

End[];

Set compilation options

Possibly set compilation options. In this case no extra options are needed.

The body of the Gauss package

Begin["`Private`"];

The GaussSolveMatlab function

GaussSolveMatlab[Real[n_,n_]@ ain_, Real[n_,m_]@ bin_ ,
 Integer@ iterations_ ]-> Real[n,m] :=
  Module [ {
    Real[n]@ 	dumc,
    Real[n,n]@ 	a,
    Real[n,m]@ 	b,
    Integer[n]@ {ipiv, indxr, indxc},
    Integer@ 	{i,k,l,irow,icol},
    Real@   	{pivinv, amax, tmp},
    Integer@ 	{beficol, afticol,  count}
   },
    
For[count=1,count<=iterations,count=count+1,(
   a=ain;  (* The arguments are always generated as const references
              and therefore cannot be changed.
              Actuall we get waste of space and time here when we copy them.   *)
   b=bin;
   For [k=1,k<=n, k=k+1,ipiv[[k]]=0];
   
   (* ipiv=0; - goes well with code generator but fails in Generated Mathematica mode *)
   For [i=1,i<=n, i=i+1,
 
     (* Algorithm first finds absolutely largest matrix element *)
     amax=0.0;
      For [k=1,k<=n, k=k+1,
       If [ ipiv[[k]]==0 ,
        For [l=1,l<=n, l=l+1,
          If [ ipiv[[l]]==0 ,
            If [ Abs[a[[k,l]]] > amax,
                    amax= Abs[a[[k,l]]];
                    irow=k;
                    icol=l
              ]]]]];
     (* Abs is generated as LIGHTABS , a macro in lightmat.h *)


     ipiv[[icol]]=ipiv[[icol]]+1; (* shortcut wanted, like   ipiv[[icol]]+=1  *)
     If[ipiv[[icol]]>1,
        "*** Gauss2 input data error ***">>""; 
        Break];
     If[ irow!=icol ,
        For [k=1,k<=n, k=k+1,
               tmp=a[[irow,k]] ;  (* Swap possibility wanted; have we? *)        
               a[[irow,k]]=a[[icol,k]]; 
               a[[icol,k]]=tmp];
        For [k=1,k<=m, k=k+1,
               tmp=b[[irow,k]] ;          
               b[[irow,k]]=b[[icol,k]]; 
               b[[icol,k]]=tmp]
        ];
      indxr[[i]]=irow;
      indxc[[i]]=icol;
      If [ a[[icol,icol]]==0,
            Print["*** Gauss2 input data error 2 ***"];
            Break];
      
      pivinv=1.0 / a[[icol,icol]];
      a[[icol,icol]]=1.0;

      a[[icol,_]]=a[[icol,_]]*pivinv;
      b[[icol,_]]=b[[icol,_]]*pivinv;
      dumc=a[[_,icol]];
      For [k=1,k<=n, k=k+1, a[[k,icol]]=0]; (* This cannot be simplified... *) 

      a[[icol,icol]]= pivinv;

      For [k=1,k<=n, k=k+1,  
       If [ k != icol,
              a[[k,_]]= a[[k,_]]- dumc[[k]]*a[[icol,_]];
              b[[k,_]]= b[[k,_]]- dumc[[k]]*b[[icol,_]]
          ]]


(* Note: OuterProd is not available in the current implementation.
 THIS IS A COMMENT, to be converted to code when OuterProd is available *)
 (*      beficol=icol-1;
       afticol=icol+1;
       a[[1| beficol, _]]= 
          a[[1| beficol,_]] - OuterProd[dumc[[1 | beficol]],a[[icol,_]]];
       b[[1| beficol, _]]= 
          b[[1| beficol, _]] -  OuterProd[dumc[[1 | beficol]],b[[icol,_]]];
       a[[afticol |_ ,_]]= 
          a[[afticol |_ ,_]] - OuterProd[dumc[[afticol | _ ]],a[[icol,_]]];
       b[[afticol |_ ,_]]= 
          b[[afticol  _, _]] - OuterProd[dumc[[afticol |_ ]],b[[icol,_]]]
*)
   ];
   (* Extra semicolons above are harmful *)

(* This fragment is present in the original algorithm, but it is
   unclear why this is really needed *)   
For [l=n,l>=1, l=l-1,
    For [k=1,k<=n, k=k+1,
     tmp= a[[k,indxr[[l]]]];
     a[[k,indxr[[l]]]]=a[[k,indxc[[l]]]];
     a[[k,indxc[[l]]]]=tmp
        ]]
   )];
  b
];

The GaussSolveForLoops function

GaussSolveForLoops[Real[n_,n_]@ ain_, Real[n_,m_]@ bin_ ,
 Integer@ iterations_ ]-> Real[n,m] :=
 Module [ {
    Real[n]@ 	dumc,
    Real[n,n]@ 	a,
    Real[n,m]@ 	b,
    Integer[n]@ {ipiv, indxr, indxc},
    Integer@ 	{i, k, l, irow, icol},
    Real@		{pivinv, amax, tmp},
    Integer@ 	{beficol, afticol, count}
  },
For[count=1,count<=iterations,count=count+1,(
   a=ain;  (* The arguments are always generated as const references
              and therefore cannot be changed.
              Actuall we get waste of space and time here when we copy them.   *)
   b=bin;
   For [k=1,k<=n, k=k+1,ipiv[[k]]=0];
   
   (* ipiv=0; - goes well with code generator but fails in Generated Mathematica mode *)
   For [i=1,i<=n, i=i+1,
 
     (* Algorithm first finds absolutely largest matrix element *)
     amax=0.0;
      For [k=1,k<=n, k=k+1,
       If [ ipiv[[k]]==0 ,
        For [l=1,l<=n, l=l+1,
          If [ ipiv[[l]]==0 ,
            If [ Abs[a[[k,l]]] > amax,
                    amax= Abs[a[[k,l]]];
                    irow=k;
                    icol=l
              ]]]]];
     (* Abs is generated as LIGHTABS , a macro in lightmat.h *)


     ipiv[[icol]]=ipiv[[icol]]+1; (* shortcut wanted, like   ipiv[[icol]]+=1  *)
     If[ipiv[[icol]]>1,
        "*** Gauss2 input data error ***">>""; 
        Break];
     If[ irow!=icol ,
        For [k=1,k<=n, k=k+1,
               tmp=a[[irow,k]] ;  (* Swap possibility wanted; have we? *)        
               a[[irow,k]]=a[[icol,k]]; 
               a[[icol,k]]=tmp];
        For [k=1,k<=m, k=k+1,
               tmp=b[[irow,k]] ;          
               b[[irow,k]]=b[[icol,k]]; 
               b[[icol,k]]=tmp]
        ];
      indxr[[i]]=irow;
      indxc[[i]]=icol;
      If [ a[[icol,icol]]==0,
            Print["*** Gauss2 input data error 2 ***"];
            Break];
      
      pivinv=1.0 / a[[icol,icol]];
      a[[icol,icol]]=1.0;

      For [k=1,k<=n, k=k+1,  
          a[[icol,k]]=a[[icol,k]]*pivinv];
      For [k=1,k<=m, k=k+1,  
          b[[icol,k]]=b[[icol,k]]*pivinv];
      For [k=1,k<=n, k=k+1,  
          dumc[[k]]=a[[k,icol]];a[[k,icol]]=0];

      
      a[[icol,icol]]= pivinv;

      For [k=1,k<=n, k=k+1,  
       If [ k != icol,
         For [l=1,l<=n, l=l+1,
              a[[k,l]]= a[[k,l]]- dumc[[k]]* a[[icol,l]]];  
             For [l=1,l<=m, l=l+1,
              b[[k,l]]= b[[k,l]]- dumc[[k]]* b[[icol,l]]] 
          ]]
   ];
   (* Extra semicolons above are harmful *)

(* Ths fragment is present in the original algorithm, but it is
   unclear why this is really needed *)   
For [l=n,l>=1, l=l-1,
    For [k=1,k<=n, k=k+1,
     tmp= a[[k,indxr[[l]]]];
     a[[k,indxr[[l]]]]=a[[k,indxc[[l]]]];
     a[[k,indxc[[l]]]]=tmp
        ]]
   )]; 
  b
];

End[];
EndPackage[];

Executing the interpreted version in Mathematica

Create two ramdom matrices.

a=Table[Random[],{10},{10}];
b=Table[Random[],{10},{2}];

Run the Matlab version

s=(c=GaussSolveMatlab[a,b,1];)//Timing;
meval=s[[1]];
Print["TIMING FOR NON-COMPILED VERSION= ",meval];
Dot[a,c] - b // MatrixForm

Run the For-loop version

s=(c=GaussSolveForLoops[a,b,1];)//Timing;
mevalFor=s[[1]];
Print["TIMING FOR NON-COMPILED VERSION (FOR-LOOPS)=",mevalFor];
Dot[a,c] - b // MatrixForm

Generation of C++ code

The command CompilePackage translates the package to C++ code:

CompilePackage["Gauss"]

Building the executable

Call MakeBinary to compile the generated C++ code and build the executable(s).

MakeBinary[]

Installing compiled code

Interpreted versions are removed, and compiled ones are used instead: (?? This should be updated. Use Activate and UnActivate instead. should not need to remove interpreted versions )

InstallCode["Gauss"]

External Execution

Preperations

SetAttributes[AbsTime,HoldFirst];
AbsTime[x_] := Module[{start,res},
		start = AbsoluteTime[];
		res=x;
		{(AbsoluteTime[]-start) Second,res}
	];

Note: Too short loop time for reliable measurements, only 3- 4 seconds combined with a resolution of 1 second.

loops=8000;
externaleval=((c=GaussSolveMatlab[a,b,loops];  )//AbsTime)[[1]]/loops

loops=750;
externalevalPass=((Do[c=GaussSolveMatlab[a,b,1],{loops}];)//AbsTime)[[1]]/loops

loops=5000;
externalevalFor=((c=GaussSolveForLoops[a,b,loops];  )//AbsTime)[[1]]/loops

loops=500;
externalevalPassFor=((Do[c=GaussSolveForLoops[a,b,1],{loops}];)//AbsTime)[[1]]/loops

loops=32000;
internalEval=((Do[c=LinearSolve[a,b],{loops}];)//AbsTime)[[1]]/loops

Dot[a,c] - b // MatrixForm

UninstallCode["Gauss"];
Map[DeleteFile,FileNames["*.o"]];Map[DeleteFile,FileNames["*.obj"]];Map[DeleteFile,FileNames["*.exe"]];

Inline C++ Compilation

Inline compilation will make matrix operations inline. This will improve performance of computation. The drawbacks are increased compilation time and size of the binary code.

MakeBinary[InlineFlag->True]

InstallCode["Gauss"];

loops=8000;
ILexternaleval=((c=GaussSolveMatlab[a,b,loops];  )//AbsTime)[[1]]/loops

loops=750;
ILexternalevalPass=((Do[c=GaussSolveMatlab[a,b,1],{loops}];)//AbsTime)[[1]]/loops

loops=5000;
ILexternalevalFor=((c=GaussSolveForLoops[a,b,loops];  )//AbsTime)[[1]]/loops

loops=500;
ILexternalevalPassFor=((Do[c=GaussSolveForLoops[a,b,1],{loops}];)//AbsTime)[[1]]/loops

UninstallCode["Gauss"];
Map[DeleteFile,FileNames["*.o"]];Map[DeleteFile,FileNames["*.obj"]];Map[DeleteFile,FileNames["*.exe"]];

Inline, without range check C++ Compilation

Inline compilation will make matrix operations inline. This will improve performance of computation. The drawbacks are increased compilation time and size of the binary code.

MakeBinary[InlineFlag->True,RangeCheckFlag->False]

InstallCode["Gauss"];

loops=32000;
RLexternaleval=((c=GaussSolveMatlab[a,b,loops];  )//AbsTime)[[1]]/loops

loops=1000;
RLexternalevalPass=((Do[c=GaussSolveMatlab[a,b,1],{loops}];)//AbsTime)[[1]]/loops

loops=32000;
RLexternalevalFor=((c=GaussSolveForLoops[a,b,loops];  )//AbsTime)[[1]]/loops

loops=500;
RLexternalevalPassFor=((Do[c=GaussSolveForLoops[a,b,1],{loops}];)//AbsTime)[[1]]/loops

a . c - b

Performance Comparison

TableForm[(({{meval, meval/RLexternaleval}, {mevalFor, mevalFor/RLexternaleval}, {"N/A", "N/A"}, {externalevalPass, externalevalPass/RLexternaleval}, {externalevalPassFor, \(externalevalPassFor/RLexternaleval  \)}, {externaleval, externaleval/RLexternaleval}, {externalevalFor, externalevalFor/RLexternaleval}, {ILexternalevalPass, ILexternalevalPass/RLexternaleval}, {ILexternalevalPassFor, ILexternalevalPassFor/RLexternaleval}, {ILexternaleval, ILexternaleval/RLexternaleval}, {ILexternalevalFor, \(ILexternalevalFor/RLexternaleval
	
                \)}, {RLexternalevalPass, RLexternalevalPass/RLexternaleval}, {RLexternalevalPassFor, RLexternalevalPassFor/RLexternaleval}, {RLexternaleval, RLexternaleval/RLexternaleval}, {RLexternalevalFor, RLexternalevalFor/RLexternaleval}, {internalEval, internalEval/RLexternaleval}})) ,\(TableHeadings -> {{"Gauss, Standard Mathematica", 
              "- GaussLoop with For loops", "Compile[]", 
              "Gauss, C++ with MathLink + call overhead", 
              "- GaussLoop with For loops + call overhead", 
              "GaussMatLab, C++ via MathLink", \ 
              "- GaussLoop with For loops", 
              "Gauss, C++ with MathLink + call overhead, Inline", 
              "- GaussLoop with For loops + call overhead, Inline", 
              "GaussMatLab, C++ via MathLink, Inline", 
              "- GaussLoop with For loops, Inline", 
			
              "Gauss, C++ with MathLink + call overhead, No range check", 
              "- GaussLoop with For loops + call overhead, No range check", "GaussMatLab, C++ via MathLink, No range check", 
              "- GaussLoop with For loops, No range check", 
              "LinearSolve"}, {"Time consumed", "Relative"}}\)]

Cleanup

UninstallCode["Gauss"]
CleanMathCodeFiles[Confirm->False]