multiply() - poc

lib/function/arithmetic/multiply.js

@@ -15,7 +15,6 @@ function factory (type, config, load, typed) {
 
   var DenseMatrix = type.DenseMatrix;
   var SparseMatrix = type.SparseMatrix;
-  var Spa = type.Spa;
 
   /**
    * Multiply two values, `x * y`. The result is squeezed.
@@ -183,15 +182,23 @@ function factory (type, config, load, typed) {
 
     // a dense
     var adata = a._data;
+    var adt = a._datatype;
     // b dense
     var bdata = b._data;
+    var bdt = b._datatype;
 
-    // result
-    var c = 0;
+    // process data types
+    var dt = adt && bdt && adt === bdt ? adt : undefined;
+    // multiply & add scalar implementations
+    var mf = dt ? multiplyScalar.signatures[dt + ',' + dt] || multiplyScalar : multiplyScalar;
+    var af = dt ? addScalar.signatures[dt + ',' + dt] || addScalar : addScalar;
+
+    // result (do not initialize it with zero)
+    var c = mf(adata[0], bdata[0]);
     // loop data
-    for (var i = 0; i < n; i++) {
+    for (var i = 1; i < n; i++) {
       // multiply and accumulate
-      c = addScalar(c, multiply(adata[i], bdata[i]));
+      c = af(c, mf(adata[i], bdata[i]));
     }
     return c;
   };
@@ -225,24 +232,32 @@ function factory (type, config, load, typed) {
     // a dense
     var adata = a._data;
     var asize = a._size;
+    var adt = a._datatype;
     // b dense
     var bdata = b._data;
     var bsize = b._size;
+    var bdt = b._datatype;
     // rows & columns
     var alength = asize[0];
     var bcolumns = bsize[1];
 
+    // process data types
+    var dt = adt && bdt && adt === bdt ? adt : undefined;
+    // multiply & add scalar implementations
+    var mf = dt ? multiplyScalar.signatures[dt + ',' + dt] || multiplyScalar : multiplyScalar;
+    var af = dt ? addScalar.signatures[dt + ',' + dt] || addScalar : addScalar;
+
     // result
     var c = new Array(bcolumns);
 
     // loop matrix columns
     for (var j = 0; j < bcolumns; j++) {
-      // sum
-      var sum = 0;      
+      // sum (do not initialize it with zero)
+      var sum = mf(adata[0], bdata[0][j]);      
       // loop vector
-      for (var i = 0; i < alength; i++) {
+      for (var i = 1; i < alength; i++) {
         // multiply & accumulate
-        sum = addScalar(sum, multiply(adata[i], bdata[i][j]));
+        sum = af(sum, mf(adata[i], bdata[i][j]));
       }
       c[j] = sum;
     }
@@ -254,7 +269,8 @@ function factory (type, config, load, typed) {
     // return matrix
     return new DenseMatrix({
       data: c,
-      size: [bcolumns]
+      size: [bcolumns],
+      datatype: dt
     });
   };
 
@@ -320,12 +336,20 @@ function factory (type, config, load, typed) {
     // a dense
     var adata = a._data;
     var asize = a._size;
+    var adt = a._datatype;
     // b dense
     var bdata = b._data;
+    var bdt = b._datatype;
     // rows & columns
     var arows = asize[0];
     var acolumns = asize[1];
 
+    // process data types
+    var dt = adt && bdt && adt === bdt ? adt : undefined;
+    // multiply & add scalar implementations
+    var mf = dt ? multiplyScalar.signatures[dt + ',' + dt] || multiplyScalar : multiplyScalar;
+    var af = dt ? addScalar.signatures[dt + ',' + dt] || addScalar : addScalar;
+
     // result
     var c = new Array(arows);
 
@@ -333,12 +357,12 @@ function factory (type, config, load, typed) {
     for (var i = 0; i < arows; i++) {
       // current row
       var row = adata[i];
-      // sum
-      var sum = 0;
+      // sum (do not initialize it with zero)
+      var sum = mf(row[0], bdata[0]);
       // loop matrix a columns
-      for (var j = 0; j < acolumns; j++) {
+      for (var j = 1; j < acolumns; j++) {
         // multiply & accumulate
-        sum = addScalar(sum, multiply(row[j], bdata[j]));
+        sum = af(sum, mf(row[j], bdata[j]));
       }
       c[i] = sum;
     }
@@ -349,7 +373,8 @@ function factory (type, config, load, typed) {
     // return matrix
     return new DenseMatrix({
       data: c,
-      size: [arows]
+      size: [arows],
+      datatype: dt
     });
   };
 
@@ -374,6 +399,7 @@ function factory (type, config, load, typed) {
     var arows = asize[0];
     var acolumns = asize[1];
     var bcolumns = bsize[1];
+
     // process data types
     var dt = adt && bdt && adt === bdt ? adt : undefined;
     // multiply & add scalar implementations
@@ -408,7 +434,8 @@ function factory (type, config, load, typed) {
     // return matrix
     return new DenseMatrix({
       data: c,
-      size: [arows, bcolumns]
+      size: [arows, bcolumns],
+      datatype: dt
     });
   };
 
@@ -429,6 +456,9 @@ function factory (type, config, load, typed) {
     var bindex = b._index;
     var bptr = b._ptr;
     var bsize = b._size;
+    // validate b matrix
+    if (!bvalues)
+      throw new Error('Cannot multiply Dense Matrix times Pattern only Matrix');
     // rows & columns
     var arows = asize[0];
     var bcolumns = bsize[1];
@@ -481,21 +511,34 @@ function factory (type, config, load, typed) {
     var avalues = a._values;
     var aindex = a._index;
     var aptr = a._ptr;
+    var adt = a._datatype;
+    // validate a matrix
+    if (!avalues)
+      throw new Error('Cannot multiply Pattern only Matrix times Dense Matrix');
     // b dense
     var bdata = b._data;
+    var bdt = b._datatype;
     // rows & columns
     var arows = a._size[0];
     var brows = b._size[0];
     // result
     var cvalues = [];
     var cindex = [];
-    var cptr = [];
+    var cptr = new Array(2);
 
-    // create sparse accumulator
-    var spa = new Spa(arows);
+    // process data types
+    var dt = adt && bdt && adt === bdt ? adt : undefined;
+    // multiply & add scalar implementations
+    var mf = dt ? multiplyScalar.signatures[dt + ',' + dt] || multiplyScalar : multiplyScalar;
+    var af = dt ? addScalar.signatures[dt + ',' + dt] || addScalar : addScalar;
+
+    // workspace
+    var x = new Array(arows);
+    // vector with marks indicating a value x[i] exists in a given column
+    var w = new Array(arows);
 
     // update ptr
-    cptr.push(0);
+    cptr[0] = 0;
     // rows in b
     for (var ib = 0; ib < brows; ib++) {
       // b[ib]
@@ -506,18 +549,31 @@ function factory (type, config, load, typed) {
         for (var ka0 = aptr[ib], ka1 = aptr[ib + 1], ka = ka0; ka < ka1; ka++) {
           // a row
           var ia = aindex[ka];
-          // accumulate
-          spa.accumulate(ia, multiply(vbi, avalues[ka]));
+          // check value exists in current j
+          if (!w[ia]) {
+            // ia is new entry in j
+            w[ia] = true;
+            // add i to pattern of C
+            cindex.push(ia);
+            // x(ia) = A
+            x[ia] = mf(vbi, avalues[ka]);
+          }
+          else {
+            // i exists in C already
+            x[ia] = af(x[ia], mf(vbi, avalues[ka]));
+          }
         }
       }
     }
-    // process spa
-    spa.forEach(0, arows - 1, function (x, v) {
-      cindex.push(x);
-      cvalues.push(v);
-    });
+    // copy values from x to column jb of c
+    for (var p1 = cindex.length, p = 0; p < p1; p++) {
+      // row
+      var ic = cindex[p];
+      // copy value
+      cvalues[p] = x[ic];
+    }
     // update ptr
-    cptr.push(cvalues.length);
+    cptr[1] = cindex.length;
 
     // check we need to squeeze the result into a scalar
     if (arows === 1)
@@ -528,7 +584,8 @@ function factory (type, config, load, typed) {
       values : cvalues,
       index: cindex,
       ptr: cptr,
-      size: [arows, 1]
+      size: [arows, 1],
+      datatype: dt
     });
   };
 
@@ -556,11 +613,13 @@ function factory (type, config, load, typed) {
     var arows = a._size[0];
     var brows = b._size[0];
     var bcolumns = b._size[1];
+
     // process data types
     var dt = adt && bdt && adt === bdt ? adt : undefined;
     // multiply & add scalar implementations
     var mf = dt ? multiplyScalar.signatures[dt + ',' + dt] || multiplyScalar : multiplyScalar;
     var af = dt ? addScalar.signatures[dt + ',' + dt] || addScalar : addScalar;
+
     // result
     var cvalues = [];
     var cindex = [];

lib/type/matrix/DenseMatrix.js

@@ -29,23 +29,28 @@ function factory (type, config, load) {
         // clone data & size
         this._data = object.clone(data._data);
         this._size = object.clone(data._size);
+        this._datatype = data._datatype;
       }
       else {
         // build data from existing matrix
         this._data = data.toArray();
         this._size = data.size();
+        this._datatype = data._datatype;
       }
     }
     else if (data && isArray(data.data) && isArray(data.size)) {
       // initialize fields from JSON representation
       this._data = data.data;
       this._size = data.size;
+      this._datatype = data.datatype;
     }
     else if (isArray(data)) {
       // replace nested Matrices with Arrays
       this._data = preprocess(data);
       // verify the size of the array, TODO: compute size while processing array
       this._size = array.size(this._data);
+      // data type unknown
+      this._datatype = undefined;
     }
     else if (data) {
       // unsupported type
@@ -55,6 +60,7 @@ function factory (type, config, load) {
       // nothing provided
       this._data = [];
       this._size = [0];
+      this._datatype = undefined;
     }
   }
   
@@ -421,7 +427,8 @@ function factory (type, config, load) {
   DenseMatrix.prototype.clone = function () {
     var m = new DenseMatrix({
       data: object.clone(this._data),
-      size: object.clone(this._size)
+      size: object.clone(this._size),
+      datatype: this._datatype
     });
     return m;
   };
@@ -529,7 +536,8 @@ function factory (type, config, load) {
     return {
       mathjs: 'DenseMatrix',
       data: this._data,
-      size: this._size
+      size: this._size,
+      datatype: this._datatype
     };
   };
   

test/type/matrix/DenseMatrix.test.js

@@ -103,7 +103,8 @@ describe('DenseMatrix', function() {
         {
           mathjs: 'DenseMatrix',
           data: [[1, 2], [3, 4]],
-          size: [2, 2]
+          size: [2, 2],
+          datatype: undefined
         });
     });
   });

tools/matrixmarket.js

@@ -195,6 +195,7 @@ var _importFromStream = function (stream, deferred) {
           var values = [];
           var index = [];
           var ptr = [];
+          var datatype = mm.datatype === 'real' ? 'number' : undefined;
           // mm data & pointer
           var d = mm.data;
           var p = -1;
@@ -238,7 +239,8 @@ var _importFromStream = function (stream, deferred) {
             values: values,
             index: index,
             ptr: ptr,
-            size: [mm.rows, mm.columns]
+            size: [mm.rows, mm.columns],
+            datatype: datatype
           }));   
           break;
         case 'array':