docs: public static double[][] transposeDiagonalTimes(double[][] m1, double[] d2, double[][] m3) {
      final int innerdim = d2.length;
      assert m1.length == innerdim : ERR_MATRIX_INNERDIM;
      assert m3.length == innerdim : ERR_MATRIX_INNERDIM;
      final int coldim1 = getColumnDimensionality(m1);
      final int coldim2 = getColumnDimensionality(m3);
      final double[][] r = new double[coldim1][coldim2];
      final double[] Acoli = new double[innerdim]; // Buffer
      // multiply it with each row from A
      for(int i = 0; i < coldim1; i++) {
        final double[] r_i = r[i]; // Output row
        // Make a linear copy of column i from A
        for(int k = 0; k < innerdim; k++) {
          Acoli[k] = m1[k][i] * d2[k];
        }
        for(int j = 0; j < coldim2; j++) {
          double s = 0;
          for(int k = 0; k < innerdim; k++) {
            s += Acoli[k] * m3[k][j];
          }
          r_i[j] = s;
        }
      }
      return r;
    }