Github user jkbradley commented on a diff in the pull request:
https://github.com/apache/spark/pull/2451#discussion_r17765173
--- Diff: mllib/src/main/scala/org/apache/spark/mllib/linalg/BLAS.scala ---
@@ -197,4 +201,452 @@ private[mllib] object BLAS extends Serializable {
throw new IllegalArgumentException(s"scal doesn't support vector
type ${x.getClass}.")
}
}
+
+ // For level-3 routines, we use the native BLAS.
+ private def nativeBLAS: NetlibBLAS = {
+ if (_nativeBLAS == null) {
+ _nativeBLAS = NativeBLAS
+ }
+ _nativeBLAS
+ }
+
+ /**
+ * C := alpha * A * B + beta * C
+ * @param transA whether to use the transpose of matrix A (true), or A
itself (false).
+ * @param transB whether to use the transpose of matrix B (true), or B
itself (false).
+ * @param alpha a scalar to scale the multiplication A * B.
+ * @param A the matrix A that will be left multiplied to B. Size of m x
k.
+ * @param B the matrix B that will be left multiplied by A. Size of k x
n.
+ * @param beta a scalar that can be used to scale matrix C.
+ * @param C the resulting matrix C. Size of m x n.
+ */
+ def gemm(
+ transA: Boolean,
+ transB: Boolean,
+ alpha: Double,
+ A: Matrix,
+ B: Matrix,
+ beta: Double,
+ C: DenseMatrix): Unit = {
+ if (alpha == 0.0) {
+ logDebug("gemm: alpha is equal to 0. Returning C.")
+ } else {
+ A match {
+ case sparse: SparseMatrix =>
+ B match {
+ case dB: DenseMatrix => gemm(transA, transB, alpha, sparse,
dB, beta, C)
+ case sB: SparseMatrix =>
+ throw new IllegalArgumentException(s"gemm doesn't support
sparse-sparse matrix " +
+ s"multiplication")
+ case _ =>
+ throw new IllegalArgumentException(s"gemm doesn't support
matrix type ${B.getClass}.")
+ }
+ case dense: DenseMatrix =>
+ B match {
+ case dB: DenseMatrix => gemm(transA, transB, alpha, dense, dB,
beta, C)
+ case sB: SparseMatrix => gemm(transA, transB, alpha, dense,
sB, beta, C)
+ case _ =>
+ throw new IllegalArgumentException(s"gemm doesn't support
matrix type ${B.getClass}.")
+ }
+ case _ =>
+ throw new IllegalArgumentException(s"gemm doesn't support matrix
type ${A.getClass}.")
+ }
+ }
+ }
+
+ /**
+ * C := alpha * A * B + beta * C
+ *
+ * @param alpha a scalar to scale the multiplication A * B.
+ * @param A the matrix A that will be left multiplied to B. Size of m x
k.
+ * @param B the matrix B that will be left multiplied by A. Size of k x
n.
+ * @param beta a scalar that can be used to scale matrix C.
+ * @param C the resulting matrix C. Size of m x n.
+ */
+ def gemm(
+ alpha: Double,
+ A: Matrix,
+ B: Matrix,
+ beta: Double,
+ C: DenseMatrix): Unit = {
+ gemm(false, false, alpha, A, B, beta, C)
+ }
+
+ /**
+ * C := alpha * A * B + beta * C
+ * For `DenseMatrix` A.
+ */
+ private def gemm(
+ transA: Boolean,
+ transB: Boolean,
+ alpha: Double,
+ A: DenseMatrix,
+ B: DenseMatrix,
+ beta: Double,
+ C: DenseMatrix): Unit = {
+ val mA: Int = if (!transA) A.numRows else A.numCols
+ val nB: Int = if (!transB) B.numCols else B.numRows
+ val kA: Int = if (!transA) A.numCols else A.numRows
+ val kB: Int = if (!transB) B.numRows else B.numCols
+ val tAstr = if (!transA) "N" else "T"
+ val tBstr = if (!transB) "N" else "T"
+
+ require(kA == kB, s"The columns of A don't match the rows of B. A:
$kA, B: $kB")
+ require(mA == C.numRows, s"The rows of C don't match the rows of A. C:
${C.numRows}, A: $mA")
+ require(nB == C.numCols,
+ s"The columns of C don't match the columns of B. C: ${C.numCols}, A:
$nB")
+
+ nativeBLAS.dgemm(tAstr, tBstr, mA, nB, kA, alpha, A.values, A.numRows,
B.values, B.numRows,
+ beta, C.values, C.numRows)
+ }
+
+ /**
+ * C := alpha * A * B + beta * C
+ * For `SparseMatrix` A.
+ */
+ private def gemm(
+ transA: Boolean,
+ transB: Boolean,
+ alpha: Double,
+ A: SparseMatrix,
+ B: DenseMatrix,
+ beta: Double,
+ C: DenseMatrix): Unit = {
+ val mA: Int = if (!transA) A.numRows else A.numCols
+ val nB: Int = if (!transB) B.numCols else B.numRows
+ val kA: Int = if (!transA) A.numCols else A.numRows
+ val kB: Int = if (!transB) B.numRows else B.numCols
+
+ require(kA == kB, s"The columns of A don't match the rows of B. A:
$kA, B: $kB")
+ require(mA == C.numRows, s"The rows of C don't match the rows of A. C:
${C.numRows}, A: $mA")
+ require(nB == C.numCols,
+ s"The columns of C don't match the columns of B. C: ${C.numCols}, A:
$nB")
+
+ val Avals = A.values
+ val Arows = if (!transA) A.rowIndices else A.colPtrs
+ val Acols = if (!transA) A.colPtrs else A.rowIndices
+
+ // Slicing is easy in this case. This is the optimal multiplication
setting for sparse matrices
+ if (transA){
+ var colCounterForB = 0
+ if (!transB) { // Expensive to put the check inside the loop
+ while (colCounterForB < nB) {
+ var rowCounterForA = 0
+ val Cstart = colCounterForB * mA
+ val Bstart = colCounterForB * kA
+ while (rowCounterForA < mA) {
+ var i = Arows(rowCounterForA)
+ val indEnd = Arows(rowCounterForA + 1)
+ var sum = 0.0
+ while (i < indEnd) {
+ sum += Avals(i) * B.values(Bstart + Acols(i))
+ i += 1
+ }
+ val Cindex = Cstart + rowCounterForA
+ C.values(Cindex) = beta * C.values(Cindex) + sum * alpha
+ rowCounterForA += 1
+ }
+ colCounterForB += 1
+ }
+ } else {
+ while (colCounterForB < nB) {
+ var rowCounter = 0
+ val Cstart = colCounterForB * mA
+ while (rowCounter < mA) {
+ var i = Arows(rowCounter)
+ val indEnd = Arows(rowCounter + 1)
+ var sum = 0.0
+ while (i < indEnd) {
+ sum += Avals(i) * B(colCounterForB, Acols(i))
+ i += 1
+ }
+ val Cindex = Cstart + rowCounter
+ C.values(Cindex) = beta * C.values(Cindex) + sum * alpha
+ rowCounter += 1
+ }
+ colCounterForB += 1
+ }
+ }
+ } else {
+ // Scale matrix first if `beta` is not equal to 0.0
+ if (beta != 0.0){
+ f2jBLAS.dscal(C.values.length, beta, C.values, 1)
+ }
+ // Perform matrix multiplication and add to C. The rows of A are
multiplied by the columns of
+ // B, and added to C.
+ var colCounterForB = 0 // the column to be updated in C
+ if (!transB) { // Expensive to put the check inside the loop
+ while (colCounterForB < nB) {
+ var colCounterForA = 0 // The column of A to multiply with the
row of B
+ val Bstart = colCounterForB * kB
+ val Cstart = colCounterForB * mA
+ while (colCounterForA < kA) {
+ var i = Acols(colCounterForA)
+ val indEnd = Acols(colCounterForA + 1)
+ val Bval = B.values(Bstart + colCounterForA) * alpha
+ while (i < indEnd){
+ C.values(Cstart + Arows(i)) += Avals(i) * Bval
+ i += 1
+ }
+ colCounterForA += 1
+ }
+ colCounterForB += 1
+ }
+ } else {
+ while (colCounterForB < nB) {
+ var colCounterForA = 0 // The column of A to multiply with the
row of B
+ val Cstart = colCounterForB * mA
+ while (colCounterForA < kA){
+ var i = Acols(colCounterForA)
+ val indEnd = Acols(colCounterForA + 1)
+ val Bval = B(colCounterForB, colCounterForA) * alpha
+ while (i < indEnd){
+ C.values(Cstart + Arows(i)) += Avals(i) * Bval
+ i += 1
+ }
+ colCounterForA += 1
+ }
+ colCounterForB += 1
+ }
+ }
+ }
+ }
+
+ /**
+ * C := alpha * A * B + beta * C
+ * For `DenseMatrix` A and `SparseMatrix` B.
+ */
+ private def gemm(
+ transA: Boolean,
+ transB: Boolean,
+ alpha: Double,
+ A: DenseMatrix,
+ B: SparseMatrix,
+ beta: Double,
+ C: DenseMatrix): Unit = {
+ val mA: Int = if (!transA) A.numRows else A.numCols
+ val nB: Int = if (!transB) B.numCols else B.numRows
+ val kA: Int = if (!transA) A.numCols else A.numRows
+ val kB: Int = if (!transB) B.numRows else B.numCols
+
+ require(kA == kB, s"The columns of A don't match the rows of B. A:
$kA, B: $kB")
+ require(mA == C.numRows, s"The rows of C don't match the rows of A. C:
${C.numRows}, A: $mA")
+ require(nB == C.numCols,
+ s"The columns of C don't match the columns of B. C: ${C.numCols}, A:
$nB")
+
+ val Bvals = B.values
+ val Brows = if (!transB) B.rowIndices else B.colPtrs
+ val Bcols = if (!transB) B.colPtrs else B.rowIndices
+
+ if (transA){
+ var colCounterForB = 0
+ if (!transB){ // Expensive to put the check inside the loop
+ while (colCounterForB < nB) {
+ var rowCounterForA = 0
+ val Cstart = colCounterForB * mA
+ val indEnd = Bcols(colCounterForB + 1)
+ while (rowCounterForA < mA) {
+ var i = Bcols(colCounterForB)
+ val Astart = rowCounterForA * kA
+ var sum = 0.0
+ while (i < indEnd) {
+ sum += Bvals(i) * A.values(Astart + Brows(i))
+ i += 1
+ }
+ val Cindex = Cstart + rowCounterForA
+ C.values(Cindex) = beta * C.values(Cindex) + sum * alpha
+ rowCounterForA += 1
+ }
+ colCounterForB += 1
+ }
+ } else {
+ var rowCounterForA = 0
+ while (rowCounterForA < mA) {
+ var colCounterForA = 0
+ val Astart = rowCounterForA * kA
+ while (colCounterForA < kA) {
+ var i = Brows(colCounterForA)
+ val indEnd = Brows(colCounterForA + 1)
+ while (i < indEnd){
+ val Cindex = Bcols(i) * mA + rowCounterForA
+ C.values(Cindex) += A.values(Astart + colCounterForA) *
Bvals(i) * alpha
+ i += 1
+ }
+ colCounterForA += 1
+ }
+ rowCounterForA += 1
+ }
+ }
+ } else {
+ // Scale matrix first if `beta` is not equal to 0.0
+ if (beta != 0.0){
+ nativeBLAS.dscal(C.values.length, beta, C.values, 1)
+ }
+ if (!transB) { // Expensive to put the check inside the loop
+
+ // Loop over the columns of B, pick non-zero row in B, select
corresponding column in A,
+ // and update the whole column in C by looping over rows in A.
+ var colCounterForB = 0 // the column to be updated in C
+ while (colCounterForB < nB) {
+ var i = Bcols(colCounterForB)
+ val indEnd = Bcols(colCounterForB + 1)
+ while (i < indEnd) {
+ var rowCounterForA = 0
+ val Bval = Bvals(i)
+ val Cstart = colCounterForB * mA
+ val Astart = mA * Brows(i)
+ while (rowCounterForA < mA){
+ C.values(Cstart + rowCounterForA) += A.values(Astart +
rowCounterForA) * Bval * alpha
+ rowCounterForA += 1
+ }
+ i += 1
+ }
+ colCounterForB += 1
+ }
+ } else {
+ var colCounterForA = 0
+ while (colCounterForA < kA) {
+ var rowCounterForA = 0
+ val Astart = mA * colCounterForA
+ val indEnd = Brows(colCounterForA + 1)
+ while (rowCounterForA < mA) {
+ var i = Brows(colCounterForA)
+ while (i < indEnd){
+ val Cindex = Bcols(i) * mA + rowCounterForA
+ C.values(Cindex) += A.values(Astart + rowCounterForA) *
Bvals(i) * alpha
+ i += 1
+ }
+ rowCounterForA += 1
+ }
+ colCounterForA += 1
+ }
+ }
+ }
+ }
+
+ /**
+ * y := alpha * A * x + beta * y
+ * @param trans whether to use the transpose of matrix A (true), or A
itself (false).
+ * @param alpha a scalar to scale the multiplication A * x.
+ * @param A the matrix A that will be left multiplied to x. Size of m x
n.
+ * @param x the vector x that will be left multiplied by A. Size of n x
1.
+ * @param beta a scalar that can be used to scale vector y.
+ * @param y the resulting vector y. Size of m x 1.
+ */
+ def gemv(
+ trans: Boolean,
+ alpha: Double,
+ A: Matrix,
+ x: DenseVector,
+ beta: Double,
+ y: DenseVector): Unit = {
+
+ val mA: Int = if (!trans) A.numRows else A.numCols
+ val nx: Int = x.size
+ val nA: Int = if (!trans) A.numCols else A.numRows
+
+ require(nA == nx, s"The columns of A don't match the number of
elements of x. A: $nA, x: $nx")
+ require(mA == y.size,
+ s"The rows of A don't match the number of elements of y. A: $mA,
y:${y.size}}")
+ if (alpha == 0.0) {
+ logDebug("gemv: alpha is equal to 0. Returning y.")
+ } else {
+ A match {
+ case sparse: SparseMatrix =>
+ gemv(trans, alpha, sparse, x, beta, y)
+ case dense: DenseMatrix =>
+ gemv(trans, alpha, dense, x, beta, y)
+ case _ =>
+ throw new IllegalArgumentException(s"gemv doesn't support matrix
type ${A.getClass}.")
+ }
+ }
+ }
+
+ /**
+ * y := alpha * A * x + beta * y
+ *
+ * @param alpha a scalar to scale the multiplication A * x.
+ * @param A the matrix A that will be left multiplied to x. Size of m x
n.
+ * @param x the vector x that will be left multiplied by A. Size of n x
1.
+ * @param beta a scalar that can be used to scale vector y.
+ * @param y the resulting vector y. Size of m x 1.
+ */
+ def gemv(
+ alpha: Double,
+ A: Matrix,
+ x: DenseVector,
+ beta: Double,
+ y: DenseVector): Unit = {
+ gemv(false, alpha, A, x, beta, y)
+ }
+
+ /**
+ * y := alpha * A * x + beta * y
+ * For `DenseMatrix` A.
+ */
+ private def gemv(
+ trans: Boolean,
+ alpha: Double,
+ A: DenseMatrix,
+ x: DenseVector,
+ beta: Double,
+ y: DenseVector): Unit = {
+ val tStrA = if (!trans) "N" else "T"
+ nativeBLAS.dgemv(tStrA, A.numRows, A.numCols, alpha, A.values,
A.numRows, x.values, 1, beta,
+ y.values, 1)
+ }
+
+ /**
+ * y := alpha * A * x + beta * y
+ * For `SparseMatrix` A.
+ */
+ private def gemv(
+ trans: Boolean,
+ alpha: Double,
+ A: SparseMatrix,
+ x: DenseVector,
+ beta: Double,
+ y: DenseVector): Unit = {
+
+ val mA: Int = if(!trans) A.numRows else A.numCols
+ val nA: Int = if(!trans) A.numCols else A.numRows
+
+ val Avals = A.values
+ val Arows = if (!trans) A.rowIndices else A.colPtrs
+ val Acols = if (!trans) A.colPtrs else A.rowIndices
+
+ // Slicing is easy in this case. This is the optimal multiplication
setting for sparse matrices
--- End diff --
Comment belongs inside "if (trans)"
---
If your project is set up for it, you can reply to this email and have your
reply appear on GitHub as well. If your project does not have this feature
enabled and wishes so, or if the feature is enabled but not working, please
contact infrastructure at infrastruct...@apache.org or file a JIRA ticket
with INFRA.
---
---------------------------------------------------------------------
To unsubscribe, e-mail: reviews-unsubscr...@spark.apache.org
For additional commands, e-mail: reviews-h...@spark.apache.org
