50g/48GX vs. 71B matrix operations--getting 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 (Ax-b):
50g/48GX:
A x * b - (Z on stack)
71B:
MAT Z=A*X @ MAT Z=Z-B
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.9E-12 !=
-.00000000001 -.00000000001
.000000000009 .000000000009
.000000000006 .000000000006
.000000000015 .000000000015
.000000000005 .000000000005
.00000000001 .00000000001
-.00000000001 -.00000000001
-1.3E-12 -1.4E-12 !=
2.1E-12 2.1E-12
-.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)=.50-RND @ NEXT I
90 FOR I=1 TO N @ FOR J=1 TO N @ A(I,J)=.50-RND @ 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=Z-B
170 L$(1)="ORDER:"&STR$(N)
180 L$(2)="EST MEM USED:"&STR$(N*(N+1)*8)
190 L$(3)="MEM USED:"&STR$(M1-M3)
200 L$(4)="FREE MEM:"&STR$(M2)
210 L$(5)="ARRAY TIME:"&STR$(T2-T1)
220 L$(6)="SOLVE TIME:"&STR$(T3-T2)
230 FIX 2
240 L$(7)="FLOPS:"&STR$(N^3/(T3-T2))
250 STD
260 L$(8)="TOTAL TIME:"&STR$(T3-T1)
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$&"Ax-b"&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$&"Ax-b"&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$&"Ax-b"&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.
|
Go back to the main exhibit hall