On Mon, Jun 08, 2009 at 03:28:50PM -0400, Alain Michaud wrote:
>
> would someone know of a procedure for solving a system of linear
> equations using a SPARSE MATRIX storage model.
>
> if one uses the usual: array [1..10000;1..10000] of "mostly_zeroes" then
> this is a really big waist of memory.
>
> I am thinking of using: Tlist instead, but I would not like to reinvent
> the wheel! I would apreciate if someone could give me a tip.
Alain,
I have attached a program (matrix.pas) to compare a full matrix against
a sparse matrix implementation (using singly linked lists with pointers)
for multiplication time. The surprising thing for me was that (back then
in 1995), the turnover point is at about 75% of nonzero (i.e. 25% zero)
elements (see mat-out.txt). This was actually a programming assignment
for our second-year students.
Note, however, that compring storage needs to be done carefully
because of the overhead in pointers and indexes.
The program is written in Turbo Pascal, but it should not be hard
to run it using FPC.
Best regards,
Tom
--
E-MAIL: T.Verhoeff @ TUE.NL | Dept. of Math. & Comp. Science
PHONE: +31 40 247 41 25 | Technische Universiteit Eindhoven
FAX: +31 40 247 54 04 | PO Box 513, NL-5600 MB Eindhoven
http://www.win.tue.nl/~wstomv/ | The Netherlands
Random Matrix2 a2:
. 0.08908 . 0.37243
. . 0.12135 .
. . 0.78422 .
. 0.19221 . .
Random Matrix1 a1 = a2:
0.00000 0.08908 0.00000 0.37243
0.00000 0.00000 0.12135 0.00000
0.00000 0.00000 0.78422 0.00000
0.00000 0.19221 0.00000 0.00000
b1 = a1 * a1:
0.00000 0.07159 0.01081 0.00000
0.00000 0.00000 0.09516 0.00000
0.00000 0.00000 0.61501 0.00000
0.00000 0.00000 0.02332 0.00000
b2 = a2 * a2:
. 0.07159 0.01081 .
. . 0.09516 .
. . 0.61501 .
. . 0.02332 .
Multiplication of random 100 x 100 Matrix1
% nonzero time in s
--------- ---------
0.00 3.30
20.00 3.63
40.00 4.29
60.00 5.49
80.00 7.19
100.00 9.39
Multiplication of random 100 x 100 Matrix2
% nonzero time in s
--------- ---------
0.00 0.00
20.00 0.71
40.00 2.31
60.00 4.45
80.00 7.58
100.00 11.26
Search turnover between 0.00% and 100.00%
Fraction = 50.00%, Mult2 takes 3.25s, Mult1 takes 4.83s
Turnover between 25.00% and 100.00%
Fraction = 62.50%, Mult2 takes 4.84s, Mult1 takes 5.71s
Turnover between 43.75% and 100.00%
Fraction = 71.88%, Mult2 takes 6.21s, Mult1 takes 6.48s
Turnover between 57.81% and 100.00%
Fraction = 78.91%, Mult2 takes 7.31s, Mult1 takes 7.14s
Turnover between 57.81% and 89.45%
Fraction = 73.63%, Mult2 takes 6.37s, Mult1 takes 6.59s
Turnover between 65.72% and 89.45%
Fraction = 77.59%, Mult2 takes 7.04s, Mult1 takes 6.92s
Turnover between 65.72% and 83.52%
Fraction = 74.62%, Mult2 takes 6.54s, Mult1 takes 6.70s
Turnover between 70.17% and 83.52%
Fraction = 76.85%, Mult2 takes 7.03s, Mult1 takes 6.92s
Turnover between 70.17% and 80.18%
Fraction = 75.18%, Mult2 takes 6.59s, Mult1 takes 6.70s
Turnover between 72.67% and 80.18%
Fraction = 76.43%, Mult2 takes 6.81s, Mult1 takes 6.81s
Turnover between 72.67% and 78.31%
Fraction = 75.49%, Mult2 takes 6.76s, Mult1 takes 6.81s
Turnover between 74.08% and 78.31%
program MatrixMultCompare;
{ Tom Verhoeff, Eindhoven University of Technology, Dept. of Math/CS,
October 1995, Version 2.0 }
{ Purpose: Compare two implementations of (sparse) matrix type, especially
w.r.t. multiplication time }
uses Timers;
type
EntryType = real;
const
Zero: EntryType = 0.0;
const
MaxN = 100; { SizeOf(EntryType) * sqr(MaxN) <= 64 KB, i.e. MaxN < 105 }
var
N: 0..MaxN;
type
Index = 0..MaxN-1; { intended 0..N-1 }
Vector1 = array [Index] of EntryType;
Matrix1 = array [Index] of Vector1;
Matrix1P = ^Matrix1; { too overcome Turbo Pascal size limitation }
Vector2 = ^Vector2Rec;
Vector2Rec = record
vix: Index; { N.B. it is not necessary to introduce upperbound here }
val: EntryType;
vnext: Vector2;
end;
{ Vector2.v == v --vnext-> nil /\ increasing in vix /\ val <> Zero }
{ [[v]].i = "if v --vnext-> w /\ w^.vix = i then w^.val else Zero"
for 0 <= i < N }
Matrix2 = ^Matrix2Rec;
Matrix2Rec = record
mix: Index;
vec: Vector2;
mnext: Matrix2;
end;
{ Matrix2.m == m --mnext-> nil /\ sorted on mix /\ vec <> nil }
{ [[v]].i = "if m --mnext-> p /\ p^.mix = i then [[p^.vec]]
else Zero-vector" for 0 <= i < N }
type
Fraction = real; { Fraction.a = 0 <= a <= 1 }
procedure WriteEntry(e: EntryType);
begin
write(e:9:5)
end; { WriteEntry }
procedure WriteZero;
begin
write('.':4, ' ':5)
end; { WriteZero }
procedure WriteZeroVec2;
begin
writeln(':':4)
end; { WriteZeroVec2 }
procedure WriteMatrix1(const a: Matrix1);
var i, j: Index;
begin
for i := 0 to N-1 do begin
for j := 0 to N-1 do WriteEntry(a[i, j]) ;
writeln
end { for i }
end; { WriteMatrix1 }
procedure WriteVector2(u: Vector2);
var j: Index;
begin
for j := 0 to N-1 do
if u = nil then WriteZero
else with u^ do
if vix = j then begin
write(val:9:5) ;
u := vnext
end { then }
else WriteZero ;
writeln
end; { WriteVector2 }
procedure WriteMatrix2(a: Matrix2);
var i: Index;
begin
for i := 0 to N-1 do
if a = nil then WriteZeroVec2
else with a^ do
if mix = i then begin
WriteVector2(vec) ;
a := mnext
end { then }
else WriteZeroVec2
end; { WriteMatrix2 }
procedure NewVector2(var v: Vector2; ix: Index; x: EntryType; next: Vector2);
{ pre: [[next]] = V /\ (A i: i<ix: V.i = Zero) ;
post: [[v]].i = if i = ix then x else V.i }
begin
if x <> Zero then begin
new(v) ;
with v^ do begin
vix := ix ;
val := x ;
vnext := next
end { with }
end { if }
end; { NewVector2 }
procedure NewMatrix2(var a: Matrix2; ix: Index; v: Vector2; next: Matrix2);
{ pre: [[next]] = A /\ (A i: i<ix: A.i = 0-vector) /\ [[v]] = V ;
post: [[a]].i = if i = ix then V else A.i }
begin
if v <> nil then begin
new(a) ;
with a^ do begin
mix := ix ;
vec := v ;
mnext := next
end { with }
end { if }
end; { NewMatrix2 }
procedure RandomMatrix1(var a: Matrix1; f: Fraction);
{ pre: true ; post: [[a]] = random N x N matrix with fraction f nonzeroes }
var i, j: Index;
begin
for i := 0 to N-1 do
for j := 0 to N-1 do
if Random < f then a[i, j] := Random else a[i, j] := 0
end; { RandomMatrix1 }
procedure RandomMatrix2(var a: Matrix2; f: Fraction);
{ pre: true ; post: [[a]] = random N x N matrix with fraction f nonzeroes }
var i, j: Index; u: Vector2;
begin
a := nil ;
for i := N-1 downto 0 do begin
u := nil ;
for j := N-1 downto 0 do
if Random < f then
NewVector2(u, j, Random, u) ;
NewMatrix2(a, i, u, a)
end { for i }
end; { RandomMatrix2 }
procedure Vector2ToVector1(u: Vector2; var v: Vector1);
{ pre: [[u]] = U ; post: [[v]] = U }
var j: Index;
begin
for j := 0 to N-1 do begin
v[j] := Zero ;
if u <> nil then
with u^ do
if vix = j then begin
v[j] := val ;
u := vnext
end { if }
end { for j }
end; { Vector2ToVector1 }
procedure Matrix2ToMatrix1(a: Matrix2; var b: Matrix1);
{ pre: [[a]] = A ; post: [[b]] = A }
var i: Index;
begin
for i := 0 to N-1 do
if a = nil then Vector2ToVector1(nil, b[i])
else with a^ do
if mix = i then begin
Vector2ToVector1(vec, b[i]) ;
a := mnext
end { then }
else Vector2ToVector1(nil, b[i])
end; { Matrix2ToMatrix1 }
procedure Matrix1ToMatrix2(const a: Matrix1; var b: Matrix2);
{ pre: [[a]] = A ; post: "b is fresh" /\ [[b]] = A }
var i, j: Index; u: Vector2;
begin
b := nil ;
for i := N-1 downto 0 do begin
u := nil ;
for j := N-1 downto 0 do
NewVector2(u, j, a[i, j], u) ;
NewMatrix2(b, i, u, b)
end { for i }
end; { Matrix1ToMatrix2 }
procedure Vector1Mult(const u: Vector1; const b: Matrix1; j: Index;
var p: EntryType) ;
{ pre: [[u]] = U /\ [[b]] = B ; post: p = U . B[_, j] }
var k: Index;
begin
p := Zero ;
for k := 0 to N-1 do p := p + u[k] * b[k, j] ;
end; { Vector1Mult }
procedure Matrix1Mult(const a, b: Matrix1; var c: Matrix1);
{ pre: [[a]] = A /\ [[b]] = B ; post: [[c]] = A * B }
var i, j: Index;
begin
for i := 0 to N-1 do
for j := 0 to N-1 do
Vector1Mult(a[i], b, j, c[i, j])
end; { Matrix1Mult }
procedure ScaleAdd(e: EntryType; u: Vector2; var v: Vector2);
{ pre: [[u]] = U /\ [[v]] = V ; post: [[v]] = e*U + V }
begin
if u <> nil then
with u^ do
if v = nil then begin
ScaleAdd(e, vnext, v) ;
NewVector2(v, vix, e*val, v)
end { then }
else if vix > v^.vix then
ScaleAdd(e, u, v^.vnext)
else if vix < v^.vix then begin
ScaleAdd(e, vnext, v) ;
NewVector2(v, vix, e*val, v)
end { then }
else { vix = v^.vix } begin
v^.val := v^.val + e*val ;
ScaleAdd(e, vnext, v^.vnext)
end { else }
end; { ScaleAdd }
procedure Vec2Mat2Mult(u: Vector2; a: Matrix2; var v: Vector2);
{ pre: [[u]] = U /\ [[a]] = A /\ [[v]] = V ; post: [[v]] = V + U * A }
begin
while (u <> nil) and (a <> nil) do
with u^, a^ do
if vix < mix then u := vnext
else if vix > mix then a := mnext
else { vix = mix } begin
ScaleAdd(val, vec, v) ;
u := vnext ;
a := mnext
end { else }
end; { Vector2Mult }
procedure Matrix2Mult(a, b: Matrix2; var c: Matrix2);
{ pre: [[a]] = A /\ [[b]] = B /\ (not Def.c \/ c=nil) ;
post: "c is fresh" /\ [[c]] = A * B }
var v: Vector2;
begin
if a = nil then c := nil
else with a^ do begin
v := nil ; { [[v]] = 0-vector }
Vec2Mat2Mult(vec, b, v) ;
Matrix2Mult(mnext, b, c) ;
NewMatrix2(c, mix, v, c)
end { with }
end; { Matrix2Mult }
procedure TestMult(f: Fraction);
{ pre: true ; post: output random NxN matrix and its square computed
both as Matrix1 and Matrix2 }
var a1, b1: Matrix1P; { dynamic variables because of large size }
a2, b2: Matrix2;
p: pointer; { for Mark/Release }
begin
Mark(p) ;
new(a1) ; new(b1) ;
RandomMatrix2(a2, f) ;
writeln('Random Matrix2 a2 with ', 100*f:1:2, '% nonzeroes:') ;
WriteMatrix2(a2) ;
Matrix2ToMatrix1(a2, a1^) ;
writeln('Random Matrix1 a1 = a2:') ;
WriteMatrix1(a1^) ;
Matrix1Mult(a1^, a1^, b1^) ;
writeln('b1 = a1 * a1:') ;
WriteMatrix1(b1^) ;
Matrix2Mult(a2, a2, b2) ;
writeln('b2 = a2 * a2:') ;
WriteMatrix2(b2) ;
Release(p)
end; { TestMult }
procedure TimeMult1(f: Fraction; var t: longint);
{ pre: true ; post: t = time for multiplying f-random Matrix1 to itself }
var
a, b: Matrix1P;
p: pointer;
tim: Timer;
k: integer;
begin
Mark(p) ;
new(a) ; new(b) ;
RandomMatrix1(a^, f) ;
tim.Start ;
Matrix1Mult(a^, a^, b^) ;
tim.Stop ;
t := tim.Elapsed(false) ;
Release(p)
end; { TimeMult1 }
procedure TimesMult1(flo, fhi: Fraction; i: integer);
{ pre: flo < fhi /\ i > 0;
post: output mult. times for i intervals evenly spaced over flo..fhi }
var k: integer; f: Fraction; t: longint;
begin
writeln('Multiplication of random ', N:1, ' x ', N:1, ' Matrix1') ;
writeln(' % nonzero time in s') ;
writeln(' --------- ---------') ;
for k := 0 to i do begin
f := ((i-k)*flo + k*fhi)/i ;
write(' ':2+1, 100*f:6:2, ' ':2+5) ;
TimeMult1(f, t) ;
writeln(' ':2, t/100:5:2)
end { for k }
end; { TimesMult1 }
procedure TimeMult2(f: Fraction; var t: longint);
{ pre: true ; post: t = time for multiplying f-random Matrix2 to itself }
var
a, b: Matrix2;
p: pointer;
tim: Timer;
begin
Mark(p) ;
RandomMatrix2(a, f) ;
tim.Start ;
Matrix2Mult(a, a, b) ;
tim.Stop ;
t := tim.Elapsed(false) ;
Release(p)
end; { TimeMult2 }
procedure TimesMult2(flo, fhi: Fraction; i: integer);
{ pre: flo < fhi /\ i > 0;
post: output mult. times for i intervals evenly spaced over flo..fhi }
var k: integer; f: Fraction; t: longint;
begin
writeln('Multiplication of random ', N:1, ' x ', N:1, ' Matrix2') ;
writeln(' % nonzero time in s') ;
writeln(' --------- ---------') ;
for k := 0 to i do begin
f := ((i-k)*flo + k*fhi)/i ;
write(' ':2+1, 100*f:6:2, ' ':2+5) ;
TimeMult2(f, t) ;
writeln(' ':2, t/100:5:2)
end { for k }
end; { TimesMult2 }
procedure FindTurnOver(hi: Fraction; eps: real; verbose: boolean);
{ pre: true;
post: ouput fraction f for which times for squaring an f-random Matrix1
and an f-random Matrix2 are equal to precision eps/2 provided f <= hi
(otherwise hi is output to precision eps/2),
if verbose then also output approximations }
{ N.B. For values N < 50 times become too small; repeated multiplication
might help; but this is not done. }
var
a1, b1: Matrix1P;
a2, b2: Matrix2;
p, q: pointer;
flo, fhi, f: Fraction;
tim: Timer;
t1, t2: longint;
begin
flo := 0.0 ; fhi := hi ;
if verbose then
writeln('Search turnover between ', 100*flo:1:2, '% and ', 100*fhi:1:2,
'%') ;
while fhi - flo > eps do begin
f := (flo + fhi) / 2 ;
if verbose then write('Fraction = ', 100*f:5:2, '%, ') ;
Mark(p) ;
RandomMatrix2(a2, f) ;
Mark(q) ;
if verbose then write('Mult2 takes ') ;
tim.Start ;
Matrix2Mult(a2, a2, b2) ;
tim.Stop ;
t2 := tim.Elapsed(false) ;
if verbose then write(t2/100:1:2, 's') ;
Release(q) ;
new(a1) ; new(b1) ;
if verbose then write(', ') ;
Matrix2ToMatrix1(a2, a1^) ;
if verbose then write('Mult1 takes ') ;
tim.Start ;
Matrix1Mult(a1^, a1^, b1^) ;
tim.Stop ;
t1 := tim.Elapsed(false) ;
if verbose then write(t1/100:1:2, 's') ;
Release(p) ;
writeln ;
{ interval reduction is done a bit more slowly than halving because
randomness of matrices makes times less reliable }
if t1 > t2 then flo := (flo + f) / 2
else fhi := (f + fhi) / 2 ;
if verbose then
writeln('Turnover between ', 100*flo:1:2, '% and ', 100*fhi:1:2, '%')
end { while } ;
writeln('Turnover fraction = ', 100*(flo+fhi)/2:1:2, ' +/- ',
100*eps/2:1:2, ' s')
end; { FindTurnOver }
begin
Randomize ;
N := 4 ;
TestMult(0.5) ;
N := MaxN ;
TimesMult1(0.0, 1.0, 5) ;
{TimesMult1(0.0, 0.2, 4) ;}
TimesMult2(0.0, 1.0, 5) ;
write('Type <return> to continue: ') ; readln ;
FindTurnOver(1.0, 0.05, true)
end.
_______________________________________________
fpc-pascal maillist - fpc-pascal@lists.freepascal.org
http://lists.freepascal.org/mailman/listinfo/fpc-pascal
