DOK Sparse Matrix Extensions
In .Net, extension methods are fantastic ways of adding new features to types that can share the same behavior (static methods) without necessarily changing the types' source codes. Extension methods can be maintained in separate code files. IDOKSparseMatrixExtensions provides extension methods to sparse matrices that implement IDOKSparseMatrix.
IDOKSparseMatrixExtensions includes some basic sparse matrix operations:
Here's the IDOKSparseMatrixExtensions code in full.
Coding matrix operations can seem complex. However, note that coordinate-wise operations with the sparse matrix are actually the simplest to code for. Coding operations for the compressed matrix formats is where it can actually be much, much tougher.
IDOKSparseMatrixExtensions includes some basic sparse matrix operations:
- Transpose() transposes the sparse matrix to a new sparse matrix.
- MatrixSum(), MatrixDifference() and MatrixAddOrSubrtact() perform addition/subtraction operation on the sparse matrix and save the results to a new sparse matrix. Observe in code that it uses 2 loops: The first loop simply copies the sparse matrix to the result sparse matrix. The 2nd loop operates on the input sparse matrix, adding to or updating the result matrix as necessary.
- MatrixScalarProduct() performs matrix-scalar multiplication and saves the results to a new sparse matrix.
- MatrixScalarProductInPlace() performs matrix-scalar multiplication and saves the results in place. Note that the code is simply a call to MatrixScalarProduct() passing the sparse matrix also as the result sparse matrix.
- MatrixVectorProduct() performs matrix-vector multiplication and returns the result as an IList<P>.
- MatrixProduct() performs matrix multiplication and saves the results to a new sparse matrix. It applies the dot product method. Observe in code the loop in loop construct, using the sparse matrix for rows and the input sparse matrix for columns, in order to build the result sparse matrix.
MDAS.cs
namespace Mendz.Library.Matrices
{
public enum MDAS
{
Multiply,
Divide,
Add,
Subtract
}
}
There are three compression methods in IDOKSparseMatrixExtensions: ToCVS(), ToCRS and ToCCS().Here's the IDOKSparseMatrixExtensions code in full.
IDOKSparseMatrixExtensions.cs
using System;
using System.Collections.Generic;
using System.Threading.Tasks;
namespace Mendz.Library.Matrices
{
public static class IDOKSparseMatrixExtensions
{
public static void Transpose<K, T>(
this IDOKSparseMatrix<K, T> sparseMatrix,
IDOKSparseMatrix<K, T> result)
{
if (!result.Size.Equals(sparseMatrix.Size))
{
throw new InvalidOperationException(
"Result matrix size must match matrix size.");
}
bool isLinearIndexed = (typeof(K) == typeof(int));
Action<KeyValuePair<K, T>> transpose = (entry) =>
{
dynamic key;
if (isLinearIndexed)
{
key = MatrixCoordinates.TransposeLinearIndex(
result.Size, (dynamic)entry.Key, result.LinearIndexMode);
}
else
{
key = MatrixCoordinates.TransposeCoordinates((dynamic)entry.Key);
}
result.SetEntry(key, entry.Value);
};
if (result.IsConcurrent)
{
Parallel.ForEach(sparseMatrix.Matrix, (entry) =>
{
transpose(entry);
});
}
else
{
foreach (var entry in sparseMatrix.Matrix)
{
transpose(entry);
}
}
}
public static void MatrixSum<K, T, A, S>(
this IDOKSparseMatrix<K, T> sparseMatrix1,
IDOKSparseMatrix<K, A> sparseMatrix2,
IDOKSparseMatrix<K, S> result)
{
MatrixAddOrSubtract(sparseMatrix1, sparseMatrix2, result, MDAS.Add);
}
public static void MatrixDifference<K, T, A, D>(
this IDOKSparseMatrix<K, T> sparseMatrix1,
IDOKSparseMatrix<K, A> sparseMatrix2,
IDOKSparseMatrix<K, D> result)
{
MatrixAddOrSubtract(sparseMatrix1, sparseMatrix2, result, MDAS.Subtract);
}
private static void MatrixAddOrSubtract<K, T, A, R>(
IDOKSparseMatrix<K, T> sparseMatrix1,
IDOKSparseMatrix<K, A> sparseMatrix2,
IDOKSparseMatrix<K, R> result, MDAS mdas)
{
if (!sparseMatrix1.Size.Equals(sparseMatrix2.Size))
{
throw new InvalidOperationException(
"Matrices must have the same size.");
}
if (!result.Size.Equals(sparseMatrix1.Size)
|| !result.Size.Equals(sparseMatrix2.Size))
{
throw new InvalidOperationException(
"Result matrix size must match matrices size.");
}
Action<KeyValuePair<K, A>> operate = (entry) =>
{
var key = entry.Key;
dynamic value;
if (result.Matrix.ContainsKey(key))
{
value = result.GetEntry(key);
if (mdas == MDAS.Subtract)
{
value -= entry.Value;
}
else
{
value += entry.Value;
}
if (value == default(R))
{
result.Matrix.Remove(key);
}
else
{
result.Matrix[key] = value;
}
}
else
{
value = entry.Value;
result.SetEntry(key, value);
}
};
if (result.IsConcurrent)
{
Parallel.ForEach(sparseMatrix1.Matrix, (entry) =>
{
result.SetEntry(entry.Key, (dynamic)entry.Value);
});
Parallel.ForEach(sparseMatrix2.Matrix, (entry) =>
{
operate(entry);
});
}
else
{
foreach (var entry in sparseMatrix1.Matrix)
{
result.SetEntry(entry.Key, (dynamic)entry.Value);
};
foreach (var entry in sparseMatrix2.Matrix)
{
operate(entry);
}
}
}
public static void MatrixScalarProduct<K, T, S, P>(
this IDOKSparseMatrix<K, T> sparseMatrix, S scalar,
IDOKSparseMatrix<K, P> result)
{
if (!result.Size.Equals(sparseMatrix.Size))
{
throw new InvalidOperationException(
"Result matrix size must match matrix size.");
}
if (scalar.Equals(default(S)))
{
result.Matrix.Clear();
}
else
{
Action<KeyValuePair<K, T>> multiply = (entry) =>
{
result.SetEntry(entry.Key, entry.Value * (dynamic)scalar);
};
if (result.IsConcurrent)
{
Parallel.ForEach(sparseMatrix.Matrix, (entry) =>
{
multiply(entry);
});
}
else
{
foreach (var entry in sparseMatrix.Matrix)
{
multiply(entry);
}
}
}
}
public static void MatrixScalarProductInPlace<K, T>(
this IDOKSparseMatrix<K, T> sparseMatrix, T scalar)
{
MatrixScalarProduct(sparseMatrix, scalar, sparseMatrix);
}
public static IList<P> MatrixVectorProduct<K, T, V, P>(
this IDOKSparseMatrix<K, T> sparseMatrix, IList<V> vector)
{
(int rows, int columns) size = sparseMatrix.Size;
int rows = size.rows;
P[] product = new P[rows];
(int row, int column) coordinates;
bool isLinearIndexed = (typeof(K) == typeof(int));
foreach (var item in sparseMatrix.Matrix)
{
dynamic key = item.Key;
if (isLinearIndexed)
{
coordinates = MatrixCoordinates.ToCoordinates(
size, key, sparseMatrix.LinearIndexMode);
}
else
{
coordinates = key;
}
product[coordinates.row] +=
item.Value * (dynamic)vector[coordinates.column];
}
return product;
}
public static void MatrixProduct<K, T, M, P>(
this IDOKSparseMatrix<K, T> sparseMatrix1,
IDOKSparseMatrix<K, M> sparseMatrix2,
IDOKSparseMatrix<K, P> result)
{
if (sparseMatrix1.Size.rows != sparseMatrix2.Size.columns)
{
throw new InvalidOperationException(
"Input matrix number of columns must match matrix number of rows.");
}
if (result.Size.rows != sparseMatrix1.Size.rows
|| result.Size.columns != sparseMatrix2.Size.columns)
{
throw new InvalidOperationException(
"Result matrix size m x n must match
matrix rows m and input matrix columns n.");
}
bool isLinearIndexed = (typeof(K) == typeof(int));
Func<(int rows, int columns), K, MatrixLinearIndexMode,
(int row, int column)> convert = (size, k, linearIndexMode) =>
{
(int row, int column) c;
if (isLinearIndexed)
{
c = MatrixCoordinates.ToCoordinates(
size, (dynamic)k, linearIndexMode);
}
else
{
c = (dynamic)k;
}
return c;
};
Action<(int row, int column), T,
(int row, int column), M> multiply = (c1, v1, c2, v2) =>
{
(int row, int column) coordinates = (c1.row, c2.column);
dynamic key;
if (isLinearIndexed)
{
key = MatrixCoordinates.ToLinearIndex(
result.Size, coordinates, result.LinearIndexMode);
}
else
{
key = coordinates;
}
P value =
result[coordinates.row, coordinates.column] + (v1 * (dynamic)v2);
if (value.Equals(default(P)))
{
result.Matrix.Remove(key);
}
else
{
result.SetEntry(key, value);
}
};
if (result.IsConcurrent)
{
Parallel.ForEach(sparseMatrix1.Matrix, (entry1) =>
{
(int row, int column) c1 = convert(
sparseMatrix1.Size, entry1.Key, sparseMatrix1.LinearIndexMode);
T v1 = entry1.Value;
Parallel.ForEach(sparseMatrix2.Matrix, (entry2) =>
{
multiply(c1, v1,
convert(
sparseMatrix2.Size, entry2.Key,
sparseMatrix2.LinearIndexMode),
entry2.Value);
});
});
}
else
{
foreach (var entry1 in sparseMatrix1.Matrix)
{
(int row, int column) c1 = convert(
sparseMatrix1.Size, entry1.Key, sparseMatrix1.LinearIndexMode);
T v1 = entry1.Value;
foreach (var entry2 in sparseMatrix2.Matrix)
{
multiply(c1, v1,
convert(
sparseMatrix2.Size, entry2.Key,
sparseMatrix2.LinearIndexMode),
entry2.Value);
}
}
}
}
public static CVS<T> ToCVS<K, T>(
this IDOKSparseMatrix<K, T> sparseMatrix,
MatrixLinearIndexMode linearIndexMode = MatrixLinearIndexMode.RowMajorOrder)
{
CVS<T> cvs = new CVS<T>(linearIndexMode, sparseMatrix.Size);
cvs.Compress(sparseMatrix);
return cvs;
}
public static CRS<T> ToCRS<K, T>(this IDOKSparseMatrix<K, T> sparseMatrix)
{
CRS<T> crs = new CRS<T>(sparseMatrix.Size);
crs.Compress(sparseMatrix);
return crs;
}
public static CCS<T> ToCCS<K, T>(this IDOKSparseMatrix<K, T> sparseMatrix)
{
CCS<T> ccs = new CCS<T>(sparseMatrix.Size);
ccs.Compress(sparseMatrix);
return ccs;
}
}
}
Quite a lot, isn't it? The purpose of the codes "(typeof(K) == typeof(int))" is to check for the sparse matrix keys' type. If it's an int, then it is consumed as a linear index. Otherwise, it's consumed as coordinates. The purpose of the codes "if (result.IsConcurrent)" for example is to check if the result sparse matrix can be used in a parallel operation. Concurrent sparse matrix types are expected to be thread-safe.Coding matrix operations can seem complex. However, note that coordinate-wise operations with the sparse matrix are actually the simplest to code for. Coding operations for the compressed matrix formats is where it can actually be much, much tougher.
Comments
Post a Comment