Skip to content

Commit

Permalink
Added arithmetic ops for sparse matrix
Browse files Browse the repository at this point in the history
  • Loading branch information
StaticBeagle committed Jan 12, 2025
1 parent e81bbe7 commit 3f5a2d6
Show file tree
Hide file tree
Showing 2 changed files with 237 additions and 29 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@

import java.util.Arrays;
import java.util.Iterator;
// TODO add arithmetic operations

import static com.wildbitsfoundry.etk4j.math.linearalgebra.ColumnCounts.adjust;

public class MatrixSparse extends Matrix {
Expand Down Expand Up @@ -415,7 +415,7 @@ public static MatrixSparse from2DArray(double[][] array) {
return from2DArray(array, ConstantsETK.DOUBLE_EPS);
}

public static MatrixSparse from2DArray(double array[][], double tol) {
public static MatrixSparse from2DArray(double[][] array, double tol) {
int nonzero = 0;
for (int i = 0; i != array.length; i++)
for (int j = 0; j != array[i].length; j++) {
Expand Down Expand Up @@ -489,20 +489,152 @@ private void incrementColumn() {
};
}

public MatrixSparse multiply(MatrixSparse A) {
MatrixSparse C = new MatrixSparse(this.rows, A.cols);
mult(this, A, C);
/**
* Performs matrix addition:<br>
* C = this + B
*
* @param B Matrix
* @return this + B
*/
public MatrixSparse add(MatrixSparse B) {
MatrixSparse outputC = new MatrixSparse(rows, cols);
if (rows != B.rows || cols != B.cols)
throw new IllegalArgumentException("Inconsistent matrix shapes.");
outputC = reshapeOrDeclare(outputC, this, rows, cols);

addOp(this, B, outputC);

return outputC;
}

public static MatrixSparse reshapeOrDeclare(MatrixSparse target, MatrixSparse reference, int rows, int cols) {
if (target == null)
return reference.create(rows, cols);
else if (target.rows != rows || target.cols != cols)
target.reshape(rows, cols);
return target;
}

/**
* Performs matrix addition:<br>
* C = A + B
*
* @param A Matrix
* @param B Matrix
* @param C Output matrix.
*/
private static void addOp(MatrixSparse A, MatrixSparse B, MatrixSparse C) {
double[] x = adjust(new DGrowArray(), A.rows);
int[] w = adjust(new IGrowArray(), A.rows, A.rows);

C.indicesSorted = false;
C.nz_length = 0;

for (int col = 0; col < A.cols; col++) {
C.col_idx[col] = C.nz_length;

multiplyAddColA(A, col, 1, C, col + 1, x, w);
multiplyAddColA(B, col, 1, C, col + 1, x, w);

// take the values in the dense vector 'x' and put them into 'C'
int idxC0 = C.col_idx[col];
int idxC1 = C.col_idx[col + 1];

for (int i = idxC0; i < idxC1; i++) {
C.nz_values[i] = x[C.nz_rows[i]];
}
}
C.col_idx[A.cols] = C.nz_length;
}

/**
* Performs matrix addition:<br>
* C = this - B
*
* @param B Matrix
* @return this - B
*/
public MatrixSparse subtract(MatrixSparse B) {
MatrixSparse outputC = new MatrixSparse(rows, cols);
if (rows != B.rows || cols != B.cols)
throw new IllegalArgumentException("Inconsistent matrix shapes.");
outputC = reshapeOrDeclare(outputC, this, rows, cols);

subtractOp(this, B, outputC);

return outputC;
}

/**
* Performs A - B
* @param A Matrix
* @param B Matrix
* @param C Output matrix.
*/
private static void subtractOp(MatrixSparse A, MatrixSparse B, MatrixSparse C) {
double[] x = adjust(new DGrowArray(), A.rows);
int[] w = adjust(new IGrowArray(), A.rows, A.rows);

C.indicesSorted = false;
C.nz_length = 0;

for (int col = 0; col < A.cols; col++) {
C.col_idx[col] = C.nz_length;

multiplyAddColA(A, col, 1, C, col + 1, x, w);
multiplyAddColA(B, col, -1, C, col + 1, x, w);

// take the values in the dense vector 'x' and put them into 'C'
int idxC0 = C.col_idx[col];
int idxC1 = C.col_idx[col + 1];

for (int i = idxC0; i < idxC1; i++) {
C.nz_values[i] = x[C.nz_rows[i]];
}
}
C.col_idx[A.cols] = C.nz_length;
}

/**
* Performs matrix multiplication. C = this * B
*
* @param B Matrix
* @return this * B
*/
public MatrixSparse multiply(MatrixSparse B) {
MatrixSparse C = new MatrixSparse(this.rows, B.cols);
multiplyOp(this, B, C);
return C;
}

/**
* Performs matrix multiplication. C = A*B
* Performs matrix multiplication. C = this * scalar
*
* @param scalar value
* @return this * scalar
*/
public MatrixSparse multiply(double scalar) {
MatrixSparse C = new MatrixSparse(rows, cols);
for (int col = 0; col < cols; col++) {
int start = col_idx[col];
int end = col_idx[col + 1];
for (int idx = start; idx < end; idx++) {
int row = nz_rows[idx];
double value = nz_values[idx];
C.unsafeSet(row, col, value * scalar);
}
}
return C;
}

/**
* Performs matrix multiplication. C = A * B
*
* @param A Matrix
* @param B Matrix
* @param C Storage for results. Array size is increased if needed.
*/
public static void mult(MatrixSparse A, MatrixSparse B, MatrixSparse C) {
private static void multiplyOp(MatrixSparse A, MatrixSparse B, MatrixSparse C) {

double[] x = adjust(new DGrowArray(), A.rows);
int[] w = adjust(new IGrowArray(), A.rows, A.rows);
Expand All @@ -527,7 +659,7 @@ public static void mult(MatrixSparse A, MatrixSparse B, MatrixSparse C) {
int rowB = B.nz_rows[bi];
double valB = B.nz_values[bi]; // B(k,j) k=rowB j=colB

multAddColA(A, rowB, valB, C, colB + 1, x, w);
multiplyAddColA(A, rowB, valB, C, colB + 1, x, w);
}

// take the values in the dense vector 'x' and put them into 'C'
Expand All @@ -547,10 +679,10 @@ public static void mult(MatrixSparse A, MatrixSparse B, MatrixSparse C) {
*
* <p>NOTE: This is the same as cs_scatter() in csparse.</p>
*/
public static void multAddColA(MatrixSparse A, int colA,
double alpha,
MatrixSparse C, int mark,
double[] x, int[] w) {
private static void multiplyAddColA(MatrixSparse A, int colA,
double alpha,
MatrixSparse C, int mark,
double[] x, int[] w) {
int idxA0 = A.col_idx[colA];
int idxA1 = A.col_idx[colA + 1];

Expand Down Expand Up @@ -579,7 +711,7 @@ public static void multAddColA(MatrixSparse A, int colA,
*/
public MatrixSparse transpose() {
MatrixSparse A_t = reshapeOrDeclare(null, cols, rows, nz_length);
transposeOp(this, A_t, null);
transposeOp(this, A_t);
return A_t;
}

Expand All @@ -591,8 +723,8 @@ public static MatrixSparse reshapeOrDeclare(MatrixSparse target, int rows, int c
return target;
}

private static void transposeOp(MatrixSparse A, MatrixSparse C, IGrowArray gw) {
int[] work = adjust(gw, A.cols, A.rows);
private static void transposeOp(MatrixSparse A, MatrixSparse C) {
int[] work = adjust((IGrowArray) null, A.cols, A.rows);
C.reshape(A.cols, A.rows, A.nz_length);

// compute the histogram for each row in 'a'
Expand All @@ -619,22 +751,17 @@ private static void transposeOp(MatrixSparse A, MatrixSparse C, IGrowArray gw) {
}
}

public static void main(String[] args) {
double[][] matrix = {
{1, 0, 0},
{0, 5, 0},
{0, 0, 10},
};

MatrixSparse sparseCSC = MatrixSparse.from2DArray(matrix, ConstantsETK.DOUBLE_EPS);
for (int col = 0; col < sparseCSC.cols; col++) {
int start = sparseCSC.col_idx[col];
int end = sparseCSC.col_idx[col + 1];
public MatrixDense toDense() {
MatrixDense matrixDense = new MatrixDense(rows, cols);
for (int col = 0; col < cols; col++) {
int start = col_idx[col];
int end = col_idx[col + 1];
for (int idx = start; idx < end; idx++) {
int row = sparseCSC.nz_rows[idx];
double value = sparseCSC.nz_values[idx];
System.out.println("Non-zero element at (" + row + ", " + col + "): " + value);
int row = nz_rows[idx];
double value = nz_values[idx];
matrixDense.unsafeSet(row, col, value);
}
}
return matrixDense;
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,87 @@ public void testUnsafeGet() {
assertEquals(100, sparseCSC.unsafeGet(0, 0), 1e-12);
}

@Test
public void testAdd() {
double[][] matrixA = {
{1, 4, 7},
{2, 5, 8},
{3, 6, 10},
};

double[][] matrixB = {
{5, 9, 13},
{3, 6, 7},
{1, 8, 10},
};

MatrixSparse A = MatrixSparse.from2DArray(matrixA, ConstantsETK.DOUBLE_EPS);
MatrixSparse B = MatrixSparse.from2DArray(matrixB, ConstantsETK.DOUBLE_EPS);
MatrixDense C = A.add(B).toDense();

double[] expected = {6, 13, 20, 5, 11, 15, 4, 14, 20};
assertArrayEquals(expected, C.getArray(), 1e-12);
}

@Test
public void testSubtract() {
double[][] matrixA = {
{1, 4, 7},
{2, 5, 8},
{3, 6, 10},
};

double[][] matrixB = {
{5, 9, 13},
{3, 6, 7},
{1, 8, 10},
};

MatrixSparse A = MatrixSparse.from2DArray(matrixA, ConstantsETK.DOUBLE_EPS);
MatrixSparse B = MatrixSparse.from2DArray(matrixB, ConstantsETK.DOUBLE_EPS);
MatrixDense C = A.subtract(B).toDense();

double[] expected = {-4, -5, -6, -1, -1, 1, 2, -2, 0};
assertArrayEquals(expected, C.getArray(), 1e-12);
}

@Test
public void testMultiply() {
double[][] matrixA = {
{1, 4, 7},
{2, 5, 8},
{3, 6, 10},
};

double[][] matrixB = {
{5, 9, 13},
{3, 6, 7},
{1, 8, 10},
};

MatrixSparse A = MatrixSparse.from2DArray(matrixA, ConstantsETK.DOUBLE_EPS);
MatrixSparse B = MatrixSparse.from2DArray(matrixB, ConstantsETK.DOUBLE_EPS);
MatrixDense C = A.multiply(B).toDense();

double[] expected = {24, 89, 111, 33, 112, 141, 43, 143, 181};
assertArrayEquals(expected, C.getArray(), 1e-12);
}

@Test
public void testMultiplyScalar() {
double[][] matrixA = {
{1, 4, 7},
{2, 5, 8},
{3, 6, 10},
};

MatrixSparse A = MatrixSparse.from2DArray(matrixA, ConstantsETK.DOUBLE_EPS);
MatrixDense C = A.multiply(2).toDense();

double[] expected = {2, 8, 14, 4, 10, 16, 6, 12, 20};
assertArrayEquals(expected, C.getArray(), 1e-12);
}

@Test
public void testLUDecomposition() {
double[][] matrix = {
Expand Down

0 comments on commit 3f5a2d6

Please sign in to comment.