2d:null 2e:["$","title",null,{"children":"Linear algebra benchmarks | Simon Frost | Observable"}] 33:T6ff,runEmlapack = (wasm, n) => { const wasmModule = emlapack_wasm const asmjsModule = emlapack_asmjs const emlapack = wasm ? wasmModule : asmjsModule const label = wasm ? 'emlapack-wasm' : 'emlapack-asmjs' const ptransa = emlapack._malloc(1) const ptransb = emlapack._malloc(1) const pm = emlapack._malloc(4) const pn = emlapack._malloc(4) const pk = emlapack._malloc(4) const palpha = emlapack._malloc(8) const pa = emlapack._malloc(8 * n * n) const plda = emlapack._malloc(4) const pb = emlapack._malloc(8 * n * n) const pldb = emlapack._malloc(4) const pbeta = emlapack._malloc(8) const pc = emlapack._malloc(8 * n * n) const pldc = emlapack._malloc(4) emlapack.setValue(ptransa, 'N'.charCodeAt(0), 'i8') emlapack.setValue(ptransb, 'N'.charCodeAt(0), 'i8') emlapack.setValue(pm, n, 'i32') emlapack.setValue(pn, n, 'i32') emlapack.setValue(pk, n, 'i32') emlapack.setValue(palpha, 1.0, 'double') emlapack.setValue(plda, n, 'i32') emlapack.setValue(pldb, n, 'i32') emlapack.setValue(pbeta, 0.0, 'double') emlapack.setValue(pldc, n, 'i32') const a = new Float64Array(emlapack.HEAPF64.buffer, pa, n * n) const b = new Float64Array(emlapack.HEAPF64.buffer, pb, n * n) const c = new Float64Array(emlapack.HEAPF64.buffer, pc, n * n) for (let i = 0; i < n; ++i) { for (let j = 0; j < n; ++j) { a[i * n + j] = Math.random() b[i * n + j] = Math.random() c[i * n + j] = 0.0 } } const dgemm = emlapack.cwrap('f2c_dgemm', null, ['number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number']) return new Benchmark.Suite() .add(label, () => { dgemm(ptransa, ptransb, pm, pn, pk, palpha, pa, plda, pb, pldb, pbeta, pc, pldc) }); }34:T52de,gpu_matrix = { function gpu_matrix() {}; gpu_matrix.prototype = { __gpumatrix__: undefined, renderbuffer : undefined, framebuffer : undefined, supportsOES_TEXTURE_FLOAT_EXTENSION : true, // multiply this into a passed matrix multiply: function(matrix) { // check conditions if (this.c != matrix.c || this.data===undefined || matrix.data===undefined) { throw "error: non-matching matrix or missing data"; } // run the web gl stuff return this._multiply(this,matrix); }, // return string representation of a matrix toString: function(radix) { if (this.data===undefined) return "undefined gpu_matrix"; var radix = radix || 10; var s = "["; var idx = 0, n = this.data.length; do { s += (idx>=this.c && idx %this.c==0) ? "\n" : ""; s += this.data[idx].toString(radix) + " "; } while(++idx=c1) { col = 0; row++; } } while(--texelcount); } return texels; }, // SUPPORTS FLOAT MATRIX -> RGBA BYTE // ### DOES NOT SUPPORT NON-SQUARE MATRICES _texelsFromMatrix: function(m,r,c) { // dimensions var texelcount = m.r*m.c; var buffer = new ArrayBuffer(4*texelcount); // get texel data (rgba) as a Float32Array var texels = new Float32Array(buffer); // copy data to Float32Array, ... var dst = 0, src1=0, src2=0; do { texels[dst++] = m.data[src1++]; } while(--texelcount); // ..., then return as IEEE754 bytes return new Uint8Array(buffer); }, // bind TEXTURE0 with m1 in .r and m2 in .g (FLOAT RGB texture) _bindDualSrcTexture: function(gl,renderer,m1,m2) { // get float array data for texture to multiply var texels = this._texelsFromMatrices(m1,m2); // create the texture from our floats var texture = gl.createTexture(); gl.activeTexture(gl.TEXTURE0); gl.bindTexture( gl.TEXTURE_2D, texture); gl.texImage2D( gl.TEXTURE_2D, /*level*/0, gl.RGB, Math.max(m1.c,m2.c), Math.max(m1.r,m2.r), 0, gl.RGB, gl.FLOAT, texels); // clamp to edge to support non-power of two textures gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE); // don't interpolate when getting data from texture gl.texParameteri( gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST); gl.texParameteri( gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST); var sampler = gl.getUniformLocation(renderer,"usampler"); gl.uniform1i(sampler,0); return texture; }, // bind passed textureUNIT to passed matrix _bindSingleSrcTexture: function(gl,renderer,m,textureUNIT,sampler) { // get float array data for texture to multiply var texels = this._texelsFromMatrix(m); // create the texture from our 4 bytes/texel (IEEE754) var texture = gl.createTexture(); gl.activeTexture(textureUNIT); gl.bindTexture( gl.TEXTURE_2D, texture); gl.texImage2D( gl.TEXTURE_2D, /*level*/0, gl.RGBA, m.c, m.r, 0, gl.RGBA, gl.UNSIGNED_BYTE, texels); // clamp to edge to support non-power of two textures gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE); // don't interpolate when getting data from texture gl.texParameteri( gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST); gl.texParameteri( gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST); // set up the appropriate sampler var sampler = gl.getUniformLocation(renderer,sampler); gl.uniform1i(sampler,textureUNIT-gl.TEXTURE0); return texture; }, // get the canvas to render to _createRenderCanvas: function(gl,m1,m2) { // Safari readPixels will not work from an 'off-screen' canvas // This code probably needs re-factoring, since the canvas should be the // first thing we grab, then the webgl context var rendercanvas = document.getElementById("matrixglcanvas"); return rendercanvas; }, // bind destination texture _createDstTexture: function(gl,rendercanvas) { // create and bind texture to render to var dstTex = gl.createTexture(); gl.activeTexture(gl.TEXTURE2); gl.bindTexture(gl.TEXTURE_2D,dstTex); gl.texImage2D(gl.TEXTURE_2D,/*level*/0,gl.RGBA,gl.RGBA,gl.UNSIGNED_BYTE,rendercanvas); return dstTex; }, // bind a framebuffer, renderbuffer, texture _bindFramebuffer: function(gl,dstTex,m1,m2) { // create and bind renderbuffer this.renderbuffer = this.renderbuffer || gl.createRenderbuffer(); gl.bindRenderbuffer(gl.RENDERBUFFER,null); gl.bindRenderbuffer(gl.RENDERBUFFER,this.renderbuffer); gl.renderbufferStorage(gl.RENDERBUFFER,gl.DEPTH_COMPONENT16,m2.c,m1.r); // create and bind framebuffer this.framebuffer = this.framebuffer || gl.createFramebuffer(); gl.bindFramebuffer(gl.FRAMEBUFFER, this.framebuffer); gl.framebufferTexture2D(gl.FRAMEBUFFER,gl.COLOR_ATTACHMENT0,gl.TEXTURE_2D,dstTex,/*level*/0); gl.framebufferRenderbuffer(gl.FRAMEBUFFER,gl.DEPTH_ATTACHMENT,gl.RENDERBUFFER,this.renderbuffer); return this.framebuffer; }, // build the glslang program to do the matrix multiply _buildRenderer: function(gl) { // get compiled shaders var vertShader = this._getShader( gl, "x-shader/x-vertex", this._vertexShader ); if (this.supportsOES_TEXTURE_FLOAT_EXTENSION) { var fragShader = this._getShader( gl, "x-shader/x-fragment", this._shader_FLOAT ); } else { // use special shader that builds own ieee754 var fragShader = this._getShader( gl, "x-shader/x-fragment", this._shader_RGBA ); } // link into a program var renderer = gl.createProgram(); gl.attachShader( renderer, vertShader ); gl.attachShader( renderer, fragShader ); gl.linkProgram( renderer ); gl.useProgram( renderer ); return renderer; }, // setup required to draw a square to our vertex shader and have // fragment shader called for each pixel _bindVertices: function(gl, renderer) { // bind vertices var aPos = gl.getAttribLocation(renderer,"aPos"); var vertexBuffer = gl.createBuffer(); gl.bindBuffer(gl.ARRAY_BUFFER, vertexBuffer); var vertices = [-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 1.0, 1.0, 0.0, -1.0, 1.0, 0.0]; gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(vertices), gl.STATIC_DRAW); gl.vertexAttribPointer(aPos, /*item size*/3, gl.FLOAT, false, 0, 0); gl.enableVertexAttribArray(aPos); // bind texture cords var aTex = gl.getAttribLocation(renderer,"aTex"); var texCoords = gl.createBuffer(); gl.bindBuffer(gl.ARRAY_BUFFER, texCoords); var textureCoords = [0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0, ]; gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(textureCoords), gl.STATIC_DRAW); gl.vertexAttribPointer(aTex, /*item size*/2, gl.FLOAT, false, 0, 0); gl.enableVertexAttribArray(aTex); // index to vertices var indices = gl.createBuffer(); gl.bindBuffer(gl.ELEMENT_ARRAY_BUFFER, indices); var vertexIndices = [0, 1, 2, 0, 2, 3]; gl.bufferData(gl.ELEMENT_ARRAY_BUFFER, new Uint16Array(vertexIndices), gl.STATIC_DRAW); }, // set up vars for the shaders _bindUniforms: function(gl,renderer,m1,m2) { // get var locations var length = gl.getUniformLocation(renderer, "uLength"); var outR = gl.getUniformLocation(renderer, "uOutRows"); var outC = gl.getUniformLocation(renderer, "uOutCols"); var stepS = gl.getUniformLocation(renderer, "uStepS"); var stepT = gl.getUniformLocation(renderer, "uStepT"); // bind length of one multiply run gl.uniform1i(length,m1.c); // bind output size // 3x1 x 1x2 -> 3x2 input and output canvas/texture // [2] x [1 1] = [2 2] called for each point in *output* texture // [3] [3 3] // [5] [5 5] gl.uniform1f(outR,m1.r); gl.uniform1f(outC,m2.c); // bind step size for input texture // 3x10 x 10x2 -> 3x2 output, but requires 10x10 *input* texture gl.uniform1f(stepS, 1./Math.max(m1.c,m2.c)); gl.uniform1f(stepT, 1./Math.max(m1.r,m2.r)); }, // multiply m1 x m2 _multiply: function(m1,m2) { // get the basics up and running var rawbuffer = new ArrayBuffer(m2.c*m1.r*4); var gl = this._checkinit(m1.r,m2.c); var renderer = this._buildRenderer(gl); if (this.supportsOES_TEXTURE_FLOAT_EXTENSION) { this._bindDualSrcTexture(gl,renderer,m1,m2); } else { this._bindSingleSrcTexture(gl,renderer,m1,gl.TEXTURE0,"usampler1"); this._bindSingleSrcTexture(gl,renderer,m2,gl.TEXTURE1,"usampler2"); } this._bindUniforms(gl,renderer,m1,m2); this._bindVertices(gl,renderer,m1,m2); // Not sure why I had this in // gl.enable(gl.DEPTH_TEST); // gl.depthFunc(gl.LEQUAL); // gl.clearDepth(1.0); // gl.clearColor(0, 0, 0, 0); // Following code was enabled, is it needed? // // Apparently not. Probably just some testing code to make sure the calls were valid. // Interestingly, leaving this code in, completely locks up Firefox 28.0, but other // older versions are just fine. NOTE: It locks up on the second press of 'Test' // Leaving it in, just in case someone wants to try and understand why FF 28.0 has an issue. // Unbind framebuffer and draw // draw to default frame buffer // gl.activeTexture(gl.TEXTURE0); // gl.activeTexture(gl.TEXTURE1); // gl.bindFramebuffer(gl.FRAMEBUFFER,null); // gl.drawElements(gl.TRIANGLES, /*num items*/6, gl.UNSIGNED_SHORT, 0); // Bind framebuffer including texture ad draw // draw to result frame buffer var rendercanvas = this._createRenderCanvas(gl,m1,m2); var dstTex = this._createDstTexture(gl,rendercanvas); var fbuffer = this._bindFramebuffer(gl,dstTex,m1,m2); // use the 6 vertex indices to draw triangles // gl.bindFramebuffer(gl.FRAMEBUFFER,null); gl.bindFramebuffer(gl.FRAMEBUFFER,fbuffer); if( gl.checkFramebufferStatus(gl.FRAMEBUFFER) != gl.FRAMEBUFFER_COMPLETE) console.log("Error: bound framebuffer is not complete."); gl.drawElements(gl.TRIANGLES, /*num items*/6, gl.UNSIGNED_SHORT, 0); // extract the product and return in new matrix var prod = new Uint8Array(rawbuffer); gl.readPixels(0,0,m2.c,m1.r,gl.RGBA,gl.UNSIGNED_BYTE,prod); return gpu_matrix.create(m1.r,m2.c,new Float32Array(rawbuffer)); }, // get shader from script tag _getShader: function( gl, shadertype, str ){ console.log("\n" +str+ "\n") // create appropriate type of shader var shader; if ( shadertype == "x-shader/x-fragment" ) shader = gl.createShader (gl.FRAGMENT_SHADER); else if ( shadertype == "x-shader/x-vertex" ) shader = gl.createShader(gl.VERTEX_SHADER); else { throw "unknown shader type " + shadertype; return null; } gl.shaderSource(shader, str); gl.compileShader(shader); if (gl.getShaderParameter(shader, gl.COMPILE_STATUS) == 0) console.log(gl.getShaderInfoLog(shader)); return shader; }, _vertexShader : " // vertex shader for a single quad \n"+ " // work is performed based on the texels being passed \n"+ " // through to the texture shader. \n"+ "#ifdef GL_ES \n"+ " precision highp float; \n"+ "#endif \n"+ " attribute vec3 aPos; \n"+ " attribute vec2 aTex; \n"+ " varying vec2 vTex; \n"+ " void main(void) \n"+ " { \n"+ " // just pass the position and texture coords \n"+ " gl_Position = vec4(aPos, 1.0); \n"+ " vTex = aTex; \n"+ " }", _shader_FLOAT : " // fragment shader that calculates the sum of the passed row and " + " // column (texture coord). \n" + " // we loop over the row and column and sum the product. \n" + " // product is then rendered to 32-bit IEEE754 floating point in the \n" + " // output RGBA canvas. \n" + " // readPixel is used to read the bytes. \n" + "#ifdef GL_ES \n" + " precision highp float; \n" + "#endif \n" + "\n" + " varying vec2 vTex; // row, column to calculate \n" + " uniform sampler2D usampler; // left in .r, right in .g \n" + " uniform int uLength; // r1xc1.r2xc2 => product has r2 (or c1) terms \n" + " uniform float uStepS; // increment across source texture \n" + " uniform float uStepT; // increment down source texture \n" + " uniform float uOutRows; // size of output in rows \n" + " uniform float uOutCols; // size of output in columns \n" + " \n" + " // sum row r x col c \n" + " float sumrowcol(float row, float col) { \n" + " float sum = 0.; // sum \n" + " float ss = 0.; // column on source texture \n" + " float tt = 0.; // row on source texture \n" + " float r = row*uStepT; // moving texture coordinate \n" + " float c = col*uStepS; // moving texture coordinate \n" + " for (int pos=0 ; pos<2048 ; ++pos) { \n" + " if(pos>=uLength) break; // stop when we multiple a row by a column \n" + " float m1 = texture2D(usampler,vec2(ss,r)).r; \n" + " float m2 = texture2D(usampler,vec2(c,tt)).g; \n" + " sum += (m1*m2); \n" + " ss += uStepS; \n" + " tt += uStepT; \n" + " } \n" + " return sum; \n" + " } \n" + " \n" + " void main(void) { \n" + " \n" + " // get the implied row and column from .s and .t of passed texel \n" + " float col = floor((vTex.s*uOutRows)); \n" + " float row = floor((vTex.t*uOutCols)); \n" + "\n" + " // sum row x col for the passed pixel \n" + " float v = sumrowcol(row,col); \n" + "\n" + " // Render to IEEE 754 Floating Point \n" + " if (v==0.) { \n" + " gl_FragColor = vec4(0.,0.,0.,0.); \n" + " return; \n" + " } \n" + " float a = abs(v); // encode absolute value + sign \n" + " float exp = floor(log2(a)); // number of powers of 2 \n" + " float mant = (a * pow(2.,23.-exp)); // multiply to fill 24 bits (implied leading 1) \n" + " float mant1 = floor(mant / 256. / 256.); // first 8 bits of mantissa \n" + " float mant2 = mod(floor(mant / 256.),256.); // second 8 bits \n" + " float mant3 = mod(mant,256.); // third 8 bits \n" + " \n" + " highp float sign = 128.-128.*(a/v); // sign bit is 256 or 0 \n" + " highp float e = (sign+exp+127.)/510.; // exponent and sign \n" + " highp float m1 = (mant1-(128.*(1.-mod(exp+127.,2.))))/255.; // handle leading bit \n" + " highp float m2 = (mant2)/255.; // middle part \n" + " highp float m3 = (mant3+.5)/255.; // scale to 0 - 255 \n" + " gl_FragColor = vec4(m3,m2,m1,e); // output an IEEE754 32-bit floating point number \n" + " } ", _shader_RGBA : " // EXPERIMENTAL: READ FLOAT DATA FROM RGBA BYTES IN IEEE754 \n" + " // fragment shader that calculates the sum of the passed row and \n" + " // column (texture coord). \n" + " // we loop over the row and column and sum the product. \n" + " // product is then rendered to 32-bit IEEE754 floating point in the \n" + " // output RGBA canvas. \n" + " // readPixel is used to read the bytes. \n" + "#ifdef GL_ES \n" + " precision highp float; \n" + "#endif \n" + "\n" + " varying vec2 vTex; // row, column to calculate \n" + " uniform sampler2D usampler1; // LEFT \n" + " uniform sampler2D usampler2; // RIGHT \n" + " uniform int uLength; // r1xc1.r2xc2 => product has r2 (or c1) terms \n" + " uniform float uStepS; // increment across source texture \n" + " uniform float uStepT; // increment down source texture \n" + " uniform float uOutRows; // size of output in rows \n" + " uniform float uOutCols; // size of output in columns \n" + "\n" + " /* \n" + " // javascript decrypt ieee754 2013/02/08 \n" + " mant3 = prod[src+0]; \n" + " mant2 = prod[src+1]; \n" + " bit = Math.floor(prod[src+2]/128.)*128; \n" + " mant1 = prod[src+2]-(bit-128.); \n" + " exp = ((prod[src+3] % 128)*2) + bit/128.; \n" + " sgn = Math.floor(prod[src+3]/128.); \n" + " f = mant1*256*256 + mant2*256 + mant3; \n" + " f = f * Math.pow(2,(exp-150))*(1-2*sgn); \n" + " */ \n" + "\n" + " float toIEEE754(vec4 bytes) { \n" + " // RETURN AN IEEE754 FLOAT FROM 4 BYTES \n" + " // GET BYTES \n" + " float byte0 = bytes.r*255.; \n" + " float byte1 = bytes.g*255.; \n" + " float byte2 = bytes.b*255.; \n" + " float byte3 = bytes.a*255.; \n" + " // COMPUTE \n" + " float mant3 = byte0; \n" + " float mant2 = byte1; \n" + " float bitv = floor(byte2/128.)*128.; \n" + " float mant1 = byte2-(bitv-128.); \n" + " float expv = (mod(byte3,128.))*2. + bitv/128.; \n" + " float sgnv = floor(byte3/128.); \n" + " float f = (mant1*256.*256.) + (mant2*256.) + mant3; \n" + " f = f * pow(2.,(expv-150.))*(1.-2.*sgnv); \n" + " return f; \n" + " } \n" + "\n" + " // sum row r x col c \n" + " float sumrowcol(float row, float col) { \n" + " float sum = 0.; // sum \n" + " float ss = 0.; // column on source texture \n" + " float tt = 0.; // row on source texture \n" + " float r = row*uStepT; // moving texture coordinate \n" + " float c = col*uStepS; // moving texture coordinate \n" + " for (int pos=0 ; pos<2048 ; ++pos) { \n" + " if(pos>=uLength) break; // stop when we multiple a row by a column \n" + " float m1 = toIEEE754(texture2D(usampler1,vec2(ss,r))); \n" + " float m2 = toIEEE754(texture2D(usampler2,vec2(c,tt))); \n" + "// used for verifying correct sampling of texture" + "// return m1; \n" + "// return float(texture2D(usampler1,vec2(ss,r)).r*255.); \n" + " sum += (m1*m2); \n" + " ss += uStepS; \n" + " tt += uStepT; \n" + " } \n" + " return sum; \n" + " } \n" + " \n" + " void main(void) { \n" + " \n" + " // get the implied row and column from .s and .t of passed texel \n" + " float col = floor((vTex.s*uOutRows)); \n" + " float row = floor((vTex.t*uOutCols)); \n" + "\n" + " // sum row x col for the passed pixel \n" + " float v = sumrowcol(row,col); \n" + "\n" + " // Render to IEEE 754 Floating Point \n" + " if (v==0.) { \n" + " gl_FragColor = vec4(0.,0.,0.,0.); \n" + " return; \n" + " } \n" + " float a = abs(v); // encode absolute value + sign \n" + " float exp = floor(log2(a)); // number of powers of 2 \n" + " float mant = (a * pow(2.,23.-exp)); // multiply to fill 24 bits (implied leading 1) \n" + " float mant1 = floor(mant / 256. / 256.); // first 8 bits of mantissa \n" + " float mant2 = mod(floor(mant / 256.),256.); // second 8 bits \n" + " float mant3 = mod(mant,256.); // third 8 bits \n" + " \n" + " highp float sign = 128.-128.*(a/v); // sign bit is 256 or 0 \n" + " highp float e = (sign+exp+127.)/510.; // exponent and sign \n" + " highp float m1 = (mant1-(128.*(1.-mod(exp+127.,2.))))/255.; // handle leading bit \n" + " highp float m2 = (mant2)/255.; // middle part \n" + " highp float m3 = (mant3+.5)/255.; // scale to 0 - 255 \n" + " gl_FragColor = vec4(m3,m2,m1,e); // output an IEEE754 32-bit floating point number \n" + " } " }; gpu_matrix.create = function(r,c,data) { var M = new gpu_matrix(); return M.setData(r,c,data); }; return gpu_matrix; }35:[] 31:{"id":"01569b2ec00fbd9a","slug":"linear-algebra-benchmarks","title":"Linear algebra benchmarks","thumbnail":"fbd1bdee690fa18a26084fcb4bb0595eed75bf9e66e973b08f7dd87898377d07","access_level":"public","create_time":"$D2018-03-21T11:04:40.935Z","update_time":"$D2018-08-14T08:11:26.040Z","view_time":null,"first_public_time":"$D2018-04-16T15:15:34.866Z","has_unpublished_edits":null,"fork":false,"stars":2,"starred":false,"owner":{"type":"team","id":"436b81191a843e22","login":"sdwfrost","name":"Simon Frost","avatar_url":"https://avatars-next.observableusercontent.com/bea598fb922216b55e8eda076a68ee10b3040a351ebe9c327ba5addd32522b64","tier":"starter_2024"},"version":663,"publish_version":663,"role":null,"authored":false,"authors":[{"type":"individual","id":"58a8b1e5ea6fc945","login":"sdwfrost","name":"Simon Frost","avatar_url":"https://avatars-next.observableusercontent.com/bea598fb922216b55e8eda076a68ee10b3040a351ebe9c327ba5addd32522b64","tier":null}],"members":[{"type":"individual","id":"58a8b1e5ea6fc945","login":"sdwfrost","name":"Simon Frost","avatar_url":"https://avatars-next.observableusercontent.com/bea598fb922216b55e8eda076a68ee10b3040a351ebe9c327ba5addd32522b64","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":7,"collection":null,"type":"notebook","body":{"title":"Linear algebra benchmarks","cells":[{"id":0,"value":"md`# Linear algebra benchmarks`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":221,"value":"{ \n var algs = mresults.map((result)=>{return result[0].name});\n var times = mresults.map((result)=>{return result[0].stats.mean*1000});\n var trace = {\n type: \"bar\",\n name: 'y',\n x: algs,\n y: times,\n }\n \n var data = [trace];\n \n var layout = {\n title: 'Time for '+m+' by '+m+' matrix',\n width: width,\n yaxis : {\n type : 'log',\n autorange : true\n },\n showlegend: false,\n font: {\n size: 16\n }\n };\n \n const div = DOM.element('div');\n Plotly.newPlot(div, data, layout);\n return div;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":89,"value":"{ \n var algs = nresults.map((result)=>{return result[0].name});\n var times = nresults.map((result)=>{return result[0].stats.mean*1000});\n var trace = {\n type: \"bar\",\n name: 'y',\n x: algs,\n y: times,\n }\n \n var data = [trace];\n \n var layout = {\n title: 'Time for '+n+ ' by '+n+' matrix',\n width: width,\n yaxis : {\n type : 'log',\n autorange : true\n },\n showlegend: false,\n font: {\n size: 16\n }\n };\n \n const div = DOM.element('div');\n Plotly.newPlot(div, data, layout);\n return div;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":87,"value":"md`## Running benchmarks`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":215,"value":"m = 2","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":218,"value":"mresults = {\n var bms = [runMathjs(m),runScijs(m),runNumjs(m),runNumericjs(m),runMatrixMultiplication(m),runLinalg(m),runLinalgAsm(m),runBluemath(m),runGpumatrix(m),runWeblas(m),runEmlapack(true,m),runEmlapack(false,m)];\n var result = [];\n for(var i=0;i {\n const a = mathjs.zeros(n, n)\n const b = mathjs.zeros(n, n)\n for (let i = 0; i < n; ++i) {\n for (let j = 0; j < n; ++j) {\n a.set([i, j], Math.random())\n b.set([i, j], Math.random())\n }\n }\n\n return new Benchmark.Suite().add('mathjs', () => {\n mathjs.multiply(a, b)\n });\n \n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":585,"value":"runNumjs = (n) => {\n const a = Create2DArray(n);\n const b = Create2DArray(n);\n for (let i = 0; i < n; ++i) {\n for (let j = 0; j < n; ++j) {\n a[i][j]=Math.random();\n b[i][j]=Math.random();\n }\n }\n const x = nj.array(a);\n const y = nj.array(b);\n \n return new Benchmark.Suite().add('numjs', () => {\n nj.dot(x, y)\n });\n \n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":186,"value":"runNumericjs = (n) => {\n const a = Create2DArray(n);\n const b = Create2DArray(n);\n for (let i = 0; i < n; ++i) {\n for (let j = 0; j < n; ++j) {\n a[i][j]=Math.random();\n b[i][j]=Math.random();\n }\n }\n\n return new Benchmark.Suite().add('numericjs', () => {\n numeric.dot(a, b)\n });\n \n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":295,"value":"runMatrixMultiplication = (n) => {\n const mul = mm()(n);\n var a = [];\n var b = [];\n for (let i = 0; i < n*n; ++i) {\n a[i]=Math.random();\n b[i]=Math.random();\n }\n\n return new Benchmark.Suite().add('matrix-multiplication', () => {\n mul(a, b)\n });\n \n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":254,"value":"runLinalg = (n) => {\n const a = Create2DArray(n);\n const b = Create2DArray(n);\n for (let i = 0; i < n; ++i) {\n for (let j = 0; j < n; ++j) {\n a[i][j]=Math.random();\n b[i][j]=Math.random();\n }\n }\n\n return new Benchmark.Suite().add('linalg', () => {\n linalg.dot(a, b)\n });\n \n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":619,"value":"runLinalgAsm = (n) => {\n \n const heap = new ArrayBuffer(24 * n * n)\n const linalgModule = linalg_asm(window, null, heap)\n const a = new Float64Array(heap, 0, n * n)\n const b = new Float64Array(heap, 8 * n * n, n * n)\n const c = new Float64Array(heap, 16 * n * n, n * n)\n \n for (let i = 0; i < n; ++i) {\n for (let j = 0; j < n; ++j) {\n a[i * n + j] = Math.random()\n b[i * n + j] = Math.random()\n c[i * n + j] = 0.0\n }\n }\n \n return new Benchmark.Suite().add('linalg-asm', () => {\n linalgModule.dgemm(0, 0, n, n, n, 1.0, a.byteOffset, n, b.byteOffset, n, 0.0, c.byteOffset, n);\n return c;\n });\n \n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":353,"value":"runBluemath = (n) => {\n const a = Create2DArray(n);\n const b = Create2DArray(n);\n for (let i = 0; i < n; ++i) {\n for (let j = 0; j < n; ++j) {\n a[i][j]=Math.random();\n b[i][j]=Math.random();\n }\n }\n \n const aa = bluemath.common.arr(a);\n const bb = bluemath.common.arr(b);\n return new Benchmark.Suite().add('bluemath', () => {\n bluemath.linalg.matmul(aa, bb)\n });\n \n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":55,"value":"runScijs = (n) => {\n const a = zeros([n, n]);\n const b = zeros([n, n]);\n const c = zeros([n, n]);\n\n ops.random(a);\n ops.random(b);\n\n return new Benchmark.Suite().add('scijs', () => {\n gemm(c, a, b);\n });\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":229,"value":"runGpumatrix = (n) => {\n const a = matrixRandom(n);\n const b = matrixRandom(n);\n \n return new Benchmark.Suite().add('gpumatrix', () => {\n a.multiply(b);\n });\n \n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":126,"value":"runGpujsGPU = (n) => {\n const a = zeros([n, n]);\n const b = zeros([n, n]);\n\n ops.random(a);\n ops.random(b);\n\n return new Benchmark.Suite().add('gpujsGPU', () => {\n gpuMultiplyMatrix(a, b, n);\n });\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":132,"value":"runGpujsCPU = (n) => {\n const a = zeros([n, n]);\n const b = zeros([n, n]);\n\n ops.random(a);\n ops.random(b);\n\n return new Benchmark.Suite().add('gpujsCPU', () => {\n cpuMultiplyMatrix(a, b, n);\n });\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":471,"value":"runWeblas = (n) => {\n var h1 = n, w1 = n, h2 = n, w2 = n;\n \n var A = new Float32Array(h1 * w1);\n for(var i=0;i {\n weblas.sgemm(M, N, K, alpha, A, B, beta, C); \n });\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":234,"value":"$33","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":13,"value":"md`## API testing matrices`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":207,"value":"a= [[1,2],[3,4]]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":367,"value":"av = [].concat.apply([], a)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":256,"value":"b = [[5,6],[7,8]]","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":369,"value":"bv = [].concat.apply([], b)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":659,"value":"af = fortranArrComplex32(av)()","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":662,"value":"bf = fortranArrComplex32(bv)()","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":580,"value":"md`## Greenlet testing`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":536,"value":"numericGreenlet = {\n let importScripts;\n var g = greenlet(async (n) => {\n importScripts('https://cdnjs.cloudflare.com/ajax/libs/numeric/1.2.6/numeric.js');\n function Create2DArray(rows) {\n var arr = [];\n for (var i=0;i {\n importScripts('https://cdnjs.cloudflare.com/ajax/libs/mathjs/4.1.2/math.min.js');\n function Create2DArray(rows) {\n var arr = [];\n for (var i=0;i {\n const mt = mathjs.transpose(m)\n const values = mt.valueOf()\n const string = values.map((row) => {\n return row.map((d) => (Math.round(d * 1000) / 1000)).join(' & ')\n }).join(\" \\\\\\\\ \")\n \n const s = `\\\\left(\\\\begin{matrix}\n ${string}\n \\\\end{matrix}\\\\right)`\n \n const ss = s.replace('/\\/\\g', '\\\\')\n \n return ss\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":396,"value":"function generateMatrices(matSize) {\n let A = [];\n let B = [];\n for (let n = 0; n < matSize * matSize; n++) {\n const randA = Math.random();\n const randB = Math.random();\n A.push(randA);\n B.push(randB);\n }\n\n A = splitArray(A, matSize);\n B = splitArray(B, matSize);\n\n function splitArray(array, part) {\n var tmp = [];\n for (var i = 0; i < array.length; i += part) {\n tmp.push(array.slice(i, i + part));\n }\n return tmp;\n }\n return {\n A,\n B\n };\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":122,"value":"function gpuMultiplyMatrix(a,b,n)\n {\n var kernel = gpu.createKernel(function(a, b, n) {\n var sum = 0;\n for (var i = 0; i < n; i++) {\n sum += a[this.thread.y][i] * b[i][this.thread.x];\n }\n return sum;\n }).setOutput([n, n]);\n \n return kernel;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":130,"value":"function cpuMultiplyMatrix(a,b,n)\n {\n var kernel = cpu.createKernel(function(a, b, n) {\n var sum = 0;\n for (var i = 0; i < n; i++) {\n sum += a[this.thread.y][i] * b[i][this.thread.x];\n }\n return sum;\n }).setOutput([n, n]);\n \n return kernel;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":393,"value":"webglChecker = {let canvas = document.createElement('canvas');\n let ctx = canvas.getContext('webgl') || canvas.getContext('experimental-webgl');\n return ctx!==null;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":15,"value":"md`## Dependencies`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":21,"value":"Benchmark = require('benchmark')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":388,"value":"greenlet = require('https://bundle.run/greenlet')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":31,"value":"md`### mathjs`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2,"value":"mathjs = require('https://bundle.run/mathjs')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":520,"value":"testmathjs = {\n var x = mathjs.matrix(a);\n var y = mathjs.matrix(b);\n return mathjs.multiply(x,y)._data;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":33,"value":"md`### scijs/ndarray`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":562,"value":"ndarray = require('https://bundle.run/ndarray')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":50,"value":"zeros = require('https://bundle.run/zeros@1.0.0')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":41,"value":"ops = require('https://bundle.run/ndarray-ops@1.2.2')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":35,"value":"gemm = require('https://bundle.run/ndarray-gemm@1.0.0')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":560,"value":"testscijs = {\n var x = ndarray(new Float64Array(av), [2,2]);\n var y = ndarray(new Float64Array(bv), [2,2]);\n var z = zeros([2,2]);\n gemm(z,x,y);\n return z;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":565,"value":"md`### numjs`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":574,"value":"nj = require('https://cdn.rawgit.com/nicolaspanel/numjs/master/dist/numjs.js')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":582,"value":"testnumjs = {\n var x = nj.array(a);\n var y = nj.array(b);\n return nj.dot(x,y);\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":112,"value":"md`### gpujs`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":114,"value":"GPU = require('https://bundle.run/gpu.js')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":116,"value":"gpu = new GPU({mode: 'gpu'})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":118,"value":"cpu = new GPU({mode: 'cpu'})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":593,"value":"cpuMultiplyMatrix(a,b,2)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":138,"value":"md`### emlapack`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":140,"value":"emlapack_asmjs = require('emlapack@0.0.3/asmjs.js').catch(() => window.Module)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":241,"value":"emlapack_wasm = require('https://unpkg.com/emlapack@0.0.3/wasm.js').catch(()=>window.Module)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":144,"value":"md`### bluemath`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":146,"value":"bluemath = require('https://bundle.run/bluemath@0.5.0')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":150,"value":"md`### matrix-multiplication`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":152,"value":"mm = require('https://bundle.run/matrix-multiplication')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":158,"value":"md`### gl-matrix`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":160,"value":"glmatrix = require('gl-matrix@2.4.0/dist/gl-matrix.js')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":154,"value":"md`### glmw`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":156,"value":"glmw = require('https://bundle.run/glmw@0.1.3')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":174,"value":"md`### gpu-matrix\n\nFrom [https://github.com/watmough/webgl-matrix-demo](https://github.com/watmough/webgl-matrix-demo)\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":327,"value":"matrixglcanvas = html`

`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":164,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":223,"value":"function matrixRandom(n) {\n var m = [];\n for(var i=0;i