1# Translated to Limbo by Eric Grosse <ehg@netlib.bell-labs.com> 3/96 2# Translated to Java by Reed Wade (wade@cs.utk.edu) 2/96 3# Translated to C by Bonnie Toy 5/88 4# Will Menninger, 10/93 5# Jack Dongarra, linpack, 3/11/78. 6# Cleve Moler, linpack, 08/14/78 7 8implement linbench; 9 10include "sys.m"; 11 sys: Sys; 12 print: import sys; 13include "math.m"; 14 math: Math; 15 dot, fabs, gemm, iamax: import math; 16include "draw.m"; 17 draw: Draw; 18 Rect, Screen, Display, Image: import draw; 19include "tk.m"; 20 tk: Tk; 21 Toplevel: import tk; 22include "tkclient.m"; 23 tkclient: Tkclient; 24include "bufio.m"; 25 bufio: Bufio; 26 Iobuf: import bufio; 27include "linalg.m"; 28 linalg: LinAlg; 29 dgefa, dgesl, printmat: import linalg; 30 31ctxt: ref Draw->Context; 32top: ref Tk->Toplevel; 33buttonlist: string; 34 35 36linbench: module 37{ 38 init: fn(nil: ref Draw->Context, argv: list of string); 39}; 40 41tkcmd(arg: string): string{ 42 rv := tk->cmd(top,arg); 43 if(rv!=nil && rv[0]=='!') 44 print("tk->cmd(%s): %s\n",arg,rv); 45 return rv; 46} 47 48init(xctxt: ref Draw->Context, nil: list of string) 49{ 50 sys = load Sys Sys->PATH; 51 math = load Math Math->PATH; 52 draw = load Draw Draw->PATH; 53 tk = load Tk Tk->PATH; 54 tkclient = load Tkclient Tkclient->PATH; 55 linalg = load LinAlg LinAlg->PATH; 56 if(linalg==nil) print("couldn't load LinAlg\n"); 57 sys->pctl(Sys->NEWPGRP, nil); 58 tkclient->init(); 59 ctxt = xctxt; 60 menubut: chan of string; 61 (top, menubut) = tkclient->toplevel(ctxt, "", 62 "Linpack in Limbo", Tkclient->Appl); 63 cmd := chan of string; 64 tk->namechan(top, cmd, "cmd"); 65 66 tkcmd("pack .Wm_t -fill x"); 67 tkcmd("frame .b"); 68 tkcmd("button .b.Run -text Run -command {send cmd run}"); 69 tkcmd("entry .b.N -width 10w"); tkcmd(".b.N insert 0 200"); 70 tkcmd("pack .b.Run .b.N -side left"); 71 tkcmd("pack .b -anchor w"); 72 tkcmd("frame .d"); 73 tkcmd("listbox .d.f -width 35w -height 150 -selectmode single -yscrollcommand {.d.fscr set}"); 74 tkcmd("scrollbar .d.fscr -command {.d.f yview}"); 75 tkcmd("pack .d.f .d.fscr -expand 1 -fill y -side right"); 76 tkcmd("pack .d -side top"); 77 tkcmd("focus .b.N"); 78 tkcmd("pack propagate . 0;update"); 79 tkclient->onscreen(top, nil); 80 tkclient->startinput(top, "kbd"::"ptr"::nil); 81 82 for(;;) alt { 83 s := <-top.ctxt.kbd => 84 tk->keyboard(top, s); 85 s := <-top.ctxt.ptr => 86 tk->pointer(top, *s); 87 s := <-top.ctxt.ctl or 88 s = <-top.wreq or 89 s = <-menubut => 90 tkclient->wmctl(top, s); 91 press := <-cmd => 92 case press { 93 "run" => 94 tkcmd("cursor -bitmap cursor.wait; update"); 95 nstr := tkcmd(".b.N get"); 96 n := int nstr; 97 (mflops,secs) := benchmark(n); 98 result := sys->sprint("%8.2f Mflops %8.1f secs",mflops,secs); 99 tkcmd("cursor -default"); 100 tkcmd(".d.f insert end {" + result + "}"); 101 tkcmd(".d.f yview moveto 1; update"); 102 } 103 } 104} 105 106benchmark(n: int): (real,real) 107{ 108 math = load Math Math->PATH; 109 110 time := array [2] of real; 111 lda := 201; 112 if(n>lda) lda = n; 113 a := array [lda*n] of real; 114 b := array [n] of real; 115 x := array [n] of real; 116 ipvt := array [n] of int; 117 ops := (2*n*n*n)/3 + 2*n*n; 118 119 norma := matgen(a,lda,n,b); 120 printmat("a",a,lda,n,n); 121 printmat("b",b,lda,n,1); 122 t1 := second(); 123 dgefa(a,lda,n,ipvt); 124 time[0] = second() - t1; 125 printmat("a",a,lda,n,n); 126 t1 = second(); 127 dgesl(a,lda,n,ipvt,b,0); 128 time[1] = second() - t1; 129 total := time[0] + time[1]; 130 131 for(i := 0; i < n; i++) { 132 x[i] = b[i]; 133 } 134 printmat("x",x,lda,n,1); 135 norma = matgen(a,lda,n,b); 136 for(i = 0; i < n; i++) { 137 b[i] = -b[i]; 138 } 139 dmxpy(b,x,a,lda); 140 resid := 0.; 141 normx := 0.; 142 for(i = 0; i < n; i++){ 143 if(resid<fabs(b[i])) resid = fabs(b[i]); 144 if(normx<fabs(x[i])) normx = fabs(x[i]); 145 } 146 147 eps_result := math->MachEps; 148 residn_result := (real n)*norma*normx*eps_result; 149 if(residn_result!=0.) 150 residn_result = resid/residn_result; 151 else 152 print("can't scale residual."); 153 if(residn_result>math->sqrt(real n)) 154 print("resid/MachEps=%.3g\n",residn_result); 155 time_result := total; 156 mflops_result := 0.; 157 if(total!=0.) 158 mflops_result = real ops/(1e6*total); 159 else 160 print("can't measure time\n"); 161 return (mflops_result,time_result); 162} 163 164 165# multiply matrix m times vector x and add the r_result to vector y. 166dmxpy(y, x, m:array of real, ldm: int) 167{ 168 n1 := len y; 169 n2 := len x; 170 gemm('N','N',n1,1,n2,1.,m,ldm,x,n2,1.,y,n1); 171} 172 173second(): real 174{ 175 return(real sys->millisec()/1000.); 176} 177 178 179# generate a (fixed) random matrix and right hand side 180# a[i][j] => a[lda*i+j] 181matgen(a: array of real, lda, n: int, b: array of real): real 182{ 183 seed := 1325; 184 norma := 0.; 185 for(j := 0; j < n; j++) 186 for(i := 0; i < n; i++){ 187 seed = 3125*seed % 65536; 188 a[lda*j+i] = (real seed - 32768.0)/16384.0; 189 if(norma < a[lda*j+i]) norma = a[lda*j+i]; 190 } 191 for (i = 0; i < n; i++) 192 b[i] = 0.; 193 for (j = 0; j < n; j++) 194 for (i = 0; i < n; i++) 195 b[i] += a[lda*j+i]; 196 return norma; 197} 198