Implementation:Haifengl Smile DenseMatrix Operations

Overview

The DenseMatrix Operations implementation provides matrix arithmetic methods on the DenseMatrix class and the Matrix interface. These methods delegate to native BLAS routines (cblas_dgemm, cblas_dgemv, cblas_daxpy, cblas_dscal, cblas_dger, and their single-precision counterparts) via the Java Foreign Function & Memory API. The implementation automatically dispatches to optimized symmetric variants (SYMM, SYMV) when the matrix is flagged as symmetric.

Type

API Doc

Source Location

• Matrix multiplication and arithmetic: base/src/main/java/smile/tensor/DenseMatrix.java:L195-921
• Matrix-vector products: base/src/main/java/smile/tensor/Matrix.java:L391-532

Import Statements

import smile.tensor.DenseMatrix;
import smile.tensor.Matrix;
import smile.tensor.Vector;
import smile.linalg.Transpose;
import smile.linalg.UPLO;

External Dependencies

Dependency	Native Functions	Integration
BLAS Level 3	cblas_dgemm, cblas_sgemm, cblas_dsymm, cblas_ssymm	FFM via smile.linalg.blas.cblas_h
BLAS Level 2	cblas_dgemv, cblas_sgemv, cblas_dsymv, cblas_ssymv, cblas_dtrmv, cblas_strmv	FFM via smile.linalg.blas.cblas_h
BLAS Level 1	cblas_daxpy, cblas_saxpy, cblas_dscal, cblas_sscal, cblas_dger, cblas_sger	FFM via smile.linalg.blas.cblas_h

API Signatures

public abstract class DenseMatrix implements Matrix, Serializable {

    // Transpose
    public DenseMatrix transpose()

    // Scaling
    public DenseMatrix scale(double alpha)

    // Matrix addition and subtraction
    public DenseMatrix add(DenseMatrix B)          // A += B
    public DenseMatrix sub(DenseMatrix B)          // A -= B
    public DenseMatrix axpy(double alpha, DenseMatrix x)  // A += alpha * x
    public void add(double alpha, DenseMatrix A, double beta, DenseMatrix B)  // C = alpha*A + beta*B

    // Matrix-matrix multiplication
    public DenseMatrix mm(DenseMatrix B)           // C = A * B
    public DenseMatrix tm(DenseMatrix B)           // C = A' * B
    public DenseMatrix mt(DenseMatrix B)           // C = A * B'
    public DenseMatrix ata()                       // C = A' * A
    public DenseMatrix aat()                       // C = A * A'

    // Static GEMM (full control)
    public static void mm(double alpha, Transpose transA, DenseMatrix A,
                          Transpose transB, DenseMatrix B, double beta, DenseMatrix C)

    // Rank-1 update
    public void ger(double alpha, Vector x, Vector y)  // A += alpha * x * y'

    // Matrix inverse
    public DenseMatrix inverse()
}

// From the Matrix interface:
public interface Matrix extends Tensor {
    // Matrix-vector multiplication
    Vector mv(Vector x)                            // y = A * x
    void mv(Vector x, Vector y)                    // y = A * x (preallocated)
    void mv(Transpose trans, double alpha, Vector x, double beta, Vector y)

    Vector tv(Vector x)                            // y = A' * x
    void tv(Vector x, Vector y)                    // y = A' * x (preallocated)

    // Quadratic form
    double xAx(Vector x)                           // x' * A * x
}

Implementation Details

GEMM Dispatch Logic

The static mm() method implements a three-way dispatch strategy:

1. If $A$ is symmetric and $B$ is not transposed: calls cblas_dsymm with Side=LEFT
2. If $B$ is symmetric and $A$ is not transposed: calls cblas_dsymm with Side=RIGHT
3. Otherwise: calls cblas_dgemm (general case)

// Symmetric dispatch example (simplified from source)
if (A.isSymmetric() && transB == NO_TRANSPOSE && B.order() == C.order()) {
    cblas_dsymm(C.order().blas(), LEFT.blas(), A.uplo().blas(), m, n,
        alpha, A.memory(), A.ld(), B.memory(), B.ld(),
        beta, C.memory(), C.ld());
} else {
    // Handle transpose flipping if layouts differ
    int k = transA == NO_TRANSPOSE ? A.ncol() : A.nrow();
    cblas_dgemm(C.order().blas(), transA.blas(), transB.blas(), m, n, k,
        alpha, A.memory(), A.ld(), B.memory(), B.ld(),
        beta, C.memory(), C.ld());
}

GEMV Dispatch Logic

The mv() method in DenseMatrix dispatches to the appropriate BLAS Level 2 routine:

Matrix Property	BLAS Routine	Condition
General	cblas_dgemv / cblas_sgemv	uplo == null
Symmetric	cblas_dsymv / cblas_ssymv	uplo != null && diag == null
Triangular (special case)	cblas_dtrmv / cblas_strmv	uplo != null && diag != null && alpha == 1 && beta == 0 && x == y

Scalar Type Safety

All arithmetic methods validate scalar type compatibility before operating:

if (scalarType() != x.scalarType()) {
    throw new IllegalArgumentException("Incompatible ScalarType: "
        + scalarType() + " != " + x.scalarType());
}

ATA/AAT Symmetry

The ata() and aat() methods automatically flag the result as symmetric:

public DenseMatrix ata() {
    DenseMatrix C = zeros(ncol(), ncol()).withUplo(LOWER);
    mm(1.0, TRANSPOSE, this, NO_TRANSPOSE, this, 0.0, C);
    return C;
}

This enables downstream operations (Cholesky, symmetric eigendecomposition) to use the optimized symmetric code path.

Usage Examples

Matrix Multiplication

import smile.tensor.DenseMatrix;
import smile.tensor.ScalarType;

// Create two matrices
DenseMatrix A = DenseMatrix.randn(ScalarType.Float64, 100, 50);
DenseMatrix B = DenseMatrix.randn(ScalarType.Float64, 50, 30);

// Standard multiplication: C = A * B
DenseMatrix C = A.mm(B);
System.out.println(C.nrow() + " x " + C.ncol()); // 100 x 30

// Transpose multiplication: C = A' * A (Gram matrix, symmetric)
DenseMatrix gram = A.ata();
System.out.println(gram.isSymmetric()); // true

// Full GEMM: C = 0.5 * A * B + 1.0 * C (accumulate into C)
DenseMatrix C2 = DenseMatrix.zeros(ScalarType.Float64, 100, 30);
DenseMatrix.mm(0.5, Transpose.NO_TRANSPOSE, A,
               Transpose.NO_TRANSPOSE, B, 1.0, C2);

Matrix-Vector Products

import smile.tensor.DenseMatrix;
import smile.tensor.Vector;
import smile.tensor.ScalarType;

DenseMatrix A = DenseMatrix.randn(ScalarType.Float64, 100, 50);
Vector x = Vector.column(new double[50]);
for (int i = 0; i < 50; i++) x.set(i, 1.0);

// y = A * x
Vector y = A.mv(x);
System.out.println(y.size()); // 100

// y = A' * x  (transpose multiply)
Vector z = A.tv(y);
System.out.println(z.size()); // 50

// Quadratic form: x' * A' * A * x
DenseMatrix AtA = A.ata();
double quad = AtA.xAx(Vector.column(new double[]{1.0, 0.0, 0.0}));

Addition, Subtraction, and Scaling

import smile.tensor.DenseMatrix;
import smile.tensor.ScalarType;

DenseMatrix A = DenseMatrix.randn(ScalarType.Float64, 10, 10);
DenseMatrix B = DenseMatrix.randn(ScalarType.Float64, 10, 10);

// A += B (in-place addition)
A.add(B);

// A -= B (in-place subtraction)
A.sub(B);

// Scale: A = 0.5 * A (in-place)
A.scale(0.5);

// AXPY: A += 2.0 * B (in-place)
A.axpy(2.0, B);

Rank-1 Update

import smile.tensor.DenseMatrix;
import smile.tensor.Vector;
import smile.tensor.ScalarType;

DenseMatrix A = DenseMatrix.zeros(ScalarType.Float64, 3, 3);
Vector x = Vector.column(new double[]{1.0, 2.0, 3.0});
Vector y = Vector.column(new double[]{4.0, 5.0, 6.0});

// A += 1.0 * x * y' (rank-1 outer product update)
A.ger(1.0, x, y);
// A is now: [[4, 5, 6], [8, 10, 12], [12, 15, 18]]

Error Handling

Condition	Exception	Method
Dimension mismatch in mm()	IllegalArgumentException	mm(DenseMatrix B)
Dimension mismatch in add()/sub()	IllegalArgumentException	axpy(double, DenseMatrix)
Incompatible ScalarType	IllegalArgumentException	All arithmetic methods
Non-square matrix for inverse()	IllegalArgumentException	inverse()
LAPACK GESV/SYSV failure	ArithmeticException	inverse()

Related

• Principle:Haifengl_Smile_Matrix_Arithmetic
• Implementation:Haifengl_Smile_DenseMatrix_Factory
• Implementation:Haifengl_Smile_DenseMatrix_Decomposition

Environment

• Environment:Haifengl_Smile_Java_25_Runtime
• Environment:Haifengl_Smile_Native_BLAS_LAPACK_ARPACK

Heuristics

• Heuristic:Haifengl_Smile_Platform_Native_Library_Selection