50g/48GX vs. 71B matrix operationsgetting slightly different results Message #1 Posted by Egan Ford on 24 Apr 2007, 11:32 p.m.
The 50g, 48GX, and the 71B all have the same Saturn processor (50g emulated). All 3 have the same machine epsilon (the smallest n where n+1=1, n>0) and the same random number generator.
First generate a large (order 20+) random square matrix and vector with values ranging from 1/2 to 1/2.
Next solve for x (Ax=b):
50g/48GX:
b A / (x on stack)
71B:
MAT x=SYS(A,b)
Compare x on both systems. They are exactly the same, as are A and b (assuming the random seeds are the same).
Next check the accuracy of the matrix operations (Axb):
50g/48GX:
A x * b  (Z on stack)
71B:
MAT Z=A*X @ MAT Z=ZB
Compare Z on both systems. Some of the values differ by 1 ULP. Why?
Output of my test:
50g/48GX 71B
 
.000000000001 .000000000001
.000000000007 .000000000007
.000000000002 .000000000002
.000000000009 .000000000009
.000000000007 .000000000007
.000000000005 .000000000005
.000000000003 2.9E12 !=
.00000000001 .00000000001
.000000000009 .000000000009
.000000000006 .000000000006
.000000000015 .000000000015
.000000000005 .000000000005
.00000000001 .00000000001
.00000000001 .00000000001
1.3E12 1.4E12 !=
2.1E12 2.1E12
.000000000017 .000000000016 !=
.000000000012 .000000000012
.000000000001 .000000000001
.000000000004 .000000000004
50g/48GX code (put N on stack. A, b, x will be saved):
%%HP: T(3)A(R)F(.);
\<< MEM DROP MEM TICKS 0. 0. 0. 0. 1. \> N M1 T1 M2 T2 T3 S E
\<< 1. RDZ 1. N
FOR I .5 RAND 
NEXT N \>ARRY 1. N
FOR I 1. N
FOR J .5 RAND 
NEXT N \>ARRY
NEXT N ROW\> MEM 'M2' STO 'A' STO 'b' STO TICKS 'T2' STO b A / 'x' STO TICKS 'T3' STO x OBJ\> OBJ\> DROP \>LIST \GSLIST 'S' STO 0. FIX "ORDER: " N + "EST MEM USED: " N N 1. + * 8. * + "MEM USED: " M1 M2  + "FREE MEM: " M2 + 2. FIX "ARRAY TIME: " T2 T1  B\>R 8192. / + "SOLVE TIME: " T3 T2  B\>R 8192. / + "FLOPS: " 2. 3. / N 3. ^ * 3. 2. / N 2. ^ * + T3 T2  B\>R 8192. / / \>NUM + "TOTAL TIME: " T3 T1  B\>R 8192. / + STD "SUM(x): " S + "eps: "
WHILE E 1. + \>NUM 1. \=/
REPEAT E 2. / \>NUM 'E' STO
END E + 6. FIX "RESIDUAL 1: " A x * b  RNRM E A CNRM N * * / + "RESIDUAL 2: " A x * b  RNRM E A CNRM x CNRM * * / + "RESIDUAL 3: " A x * b  RNRM E A RNRM x RNRM * * / +
\>> 13. \>LIST 1.
\<< 10. CHR +
\>> DOSUBS \GSLIST 'RM.TXT' STO STD
\>>
71B code (Enter N when prompted. A, B, X, Z in RAM):
10 STD @ DESTROY ALL @ OPTION BASE 1 @ RANDOMIZE 1
20 H=13 @ DIM L$(H)[80]
30 M=INT(SQR(MEM/2/8))1
40 INPUT "N=",STR$(M);N
50 M1=MEM @ REAL A(N,N),B(N,1) @ M2=MEM @ REAL X(N,1) @ M3=MEM
60 DISP "INITIALIZING ARRAYS..."
70 T1=TIME
80 FOR I=1 TO N @ B(I,1)=.50RND @ NEXT I
90 FOR I=1 TO N @ FOR J=1 TO N @ A(I,J)=.50RND @ NEXT J @ NEXT I
100 DISP "SOLVE..."
110 T2=TIME @ MAT X=SYS(A,B) @ T3=TIME @ BEEP
120 S=0 @ FOR I=1 TO N @ S=S+X(I,1) @ NEXT I
130 E=1
140 IF E+1#1 THEN E=E/2 @ GOTO 140
150 REAL Z(N,1)
160 MAT Z=A*X @ MAT Z=ZB
170 L$(1)="ORDER:"&STR$(N)
180 L$(2)="EST MEM USED:"&STR$(N*(N+1)*8)
190 L$(3)="MEM USED:"&STR$(M1M3)
200 L$(4)="FREE MEM:"&STR$(M2)
210 L$(5)="ARRAY TIME:"&STR$(T2T1)
220 L$(6)="SOLVE TIME:"&STR$(T3T2)
230 FIX 2
240 L$(7)="FLOPS:"&STR$(N^3/(T3T2))
250 STD
260 L$(8)="TOTAL TIME:"&STR$(T3T1)
270 L$(9)="SUM(x):"&STR$(S)
280 L$(10)="eps:"&STR$(E)
290 FIX 6 ! STD ! FIX 6
300 V$=CHR$(124)&CHR$(124)
310 U$=CHR$(95)
320 L$(11)=V$&"Axb"&V$&U$&"oo/(eps*"&V$&"A"&V$&U$&"1*N)= "
330 L$(11)=L$(11)&STR$(RNORM(Z)/(E*CNORM(A)*N))
340 L$(12)=V$&"Axb"&V$&U$&"oo/(eps*"&V$&"A"&V$&U$&"1*"&V$&"x"&V$&U$&"1)= "
350 L$(12)=L$(12)&STR$(RNORM(Z)/(E*CNORM(A)*CNORM(X)))
360 L$(13)=V$&"Axb"&V$&U$&"oo/(eps*"&V$&"A"&V$&U$&"oo*"&V$&"x"&V$&U$&"oo)="
370 L$(13)=L$(13)&STR$(RNORM(Z)/(E*RNORM(A)*RNORM(X)))
380 STD
385 FOR I=1 TO H @ DISP L$(I) @ NEXT I @ GOTO 400
390 CALL SCROLLW(L$(),H)
400 ! DESTROY A,B,X,Z,L
410 END
Edited: 25 Apr 2007, 12:16 a.m.
