Computer Science Experimentation

Friday, February 18, 2011

F# Linear Algebra with Blas and Lapack

I am testing the F# MathProvider library (http://mathprovider.codeplex.com/).
Mathprovider site has a very interesting section called "About the Performance" that I recommend.


For now, I tried the following sample code with no problems.
// reference F# PowerPack & MathProvider
#r @"FSharp.PowerPack.dll"
//#r @"D:\FSharp\math-provider\MathProvider\MathProvider\bin\Release\MathProvider.dll"
#r @"C:\cs\Others\MathProvider\Net 4.0\MathProvider.dll"
#time

// rename the module name
module L = MathProvider.LinearAlgebra

// set the directory of the Lapack runtime (blas.dll & lapack.dll)
System.Environment.CurrentDirectory <- @"C:\cs\Others\MathProvider\LapackRuntime"
//System.Environment.CurrentDirectory <- @"D:\FSharp\math-provider\LapackRuntime\Netlib" 
//System.Environment.CurrentDirectory <- @"D:\FSharp\math-provider\LapackRuntime\MKL"

// start the native provider, otherwise (implementent) F# implementation will be used
L.startProvider()

// create two 3x3 matrices
let A = matrix [ [12.; -51.; 4.; ]; [6.; 167.; -68.;]; [-4.; 24.; -41.; ] ]
let B = matrix [ [ 2.; 1.; 1.;] ; [ 1.; 2.; 1.;]; [ 1.; 1.; 2.;] ]

// determinant 
let det = L.det A
// inverse
let inv = L.inv A
// qr decomposition
let q, r  = L.qr  A
// lu decomposition
let p, l, u  = L.lu  A
// cholesky decomposition
let ch  = L.chol B
// svd decomposition
let v, s, ut = L.svd A
// eigen decomposition for symetric matrix
let a, b = L.cov A A  |> L.eigenSym 

(* result *)

(*
val A : matrix = matrix [[12.0; -51.0; 4.0]
                         [6.0; 167.0; -68.0]
                         [-4.0; 24.0; -41.0]]
val B : matrix = matrix [[2.0; 1.0; 1.0]
                         [1.0; 2.0; 1.0]
                         [1.0; 1.0; 2.0]]
val det : float = -85750.0
val inv : matrix = matrix [[0.06081632653; 0.02326530612; -0.03265306122]
                           [-0.006040816327; 0.005551020408; -0.009795918367]
                           [-0.009469387755; 0.0009795918367; -0.02693877551]]
val r : matrix = matrix [[-14.0; -21.0; 14.0]
                         [0.0; -175.0; 70.0]
                         [0.0; 0.0; -35.0]]
val q : matrix = matrix [[-0.8571428571; 0.3942857143; 0.3314285714]
                         [-0.4285714286; -0.9028571429; -0.03428571429]
                         [0.2857142857; -0.1714285714; 0.9428571429]]
val u : matrix = matrix [[12.0; -51.0; 4.0]
                         [0.0; 192.5; -70.0]
                         [0.0; 0.0; -37.12121212]]
val p : (int -> int)
val l : matrix = matrix [[1.0; 0.0; 0.0]
                         [0.5; 1.0; 0.0]
                         [-0.3333333333; 0.03636363636; 1.0]]
val ch : matrix = matrix [[1.414213562; 0.7071067812; 0.7071067812]
                          [0.0; 1.224744871; 0.4082482905]
                          [0.0; 0.0; 1.154700538]]
val v : matrix = matrix [[-0.2543778627; -0.5139835724; -0.81921474]
                         [0.9464104305; 0.0419982359; -0.3202240549]
                         [0.1989954776; -0.8567712854; 0.4757559925]]
val ut : matrix = matrix [[0.00960262784; 0.922507462; -0.385859783]
                          [-0.07574450365; 0.3854399898; 0.9196188256]
                          [-0.9970810196; -0.0203960004; -0.0735761066]]
val s : Vector = vector [|190.5672437; 32.85688323; 13.69492038|]
val b : Matrix = matrix [[0.9970810196; -0.07574450365; -0.00960262784]
                                [0.0203960004; 0.3854399898; -0.922507462]
                                [0.0735761066; 0.9196188256; 0.385859783]]
val a : vector = vector [|187.5508443; 1079.574776; 36315.87438|]
*)

No comments:

Post a Comment