docs: public static double[][] times(final double[][] m1, final double[][] m2) {
      final int rowdim1 = m1.length, coldim1 = getColumnDimensionality(m1);
      final int coldim2 = getColumnDimensionality(m2);
      // Optimized implementation, exploiting the storage layout
      assert m2.length == coldim1 : ERR_MATRIX_INNERDIM;
      final double[][] r2 = new double[rowdim1][coldim2];
      // Optimized ala Jama. jik order.
      final double[] Bcolj = new double[coldim1];
      for(int j = 0; j < coldim2; j++) {
        // Make a linear copy of column j from B
        for(int k = 0; k < coldim1; k++) {
          Bcolj[k] = m2[k][j];
        }
        // multiply it with each row from A
        for(int i = 0; i < rowdim1; i++) {
          final double[] Arowi = m1[i];
          double s = 0;
          for(int k = 0; k < coldim1; k++) {
            s += Arowi[k] * Bcolj[k];
          }
          r2[i][j] = s;
        }
      }
      return r2;
    }