youmuliu/stgsharp
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

swiss ok, needs more benchmark

Browse Source
This commit is contained in:
Youmu Liu
2026-03-06 18:22:30 +08:00
parent 83997f64d0
commit f6e897ed06
70 changed files with 1889 additions and 3371 deletions
Download Patch File Download Diff File Expand all files Collapse all files
StgSharp.Mathematics/Numeric/Linear/Exception.cs
+65
View File
@@ -0,0 +1,65 @@
//-----------------------------------------------------------------------
// -----------------------------------------------------------------------
// file="Exception"
// Project: StgSharp
// AuthorGroup: Nitload
// Copyright (c) Nitload. All rights reserved.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//
// -----------------------------------------------------------------------
// -----------------------------------------------------------------------
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;
namespace StgSharp.Mathematics.Numeric
{
public class MatrixDimensionMismatchException<T> : Exception where T : unmanaged, INumber<T>
{
public MatrixDimensionMismatchException(
Matrix<T> right,
Matrix<T> ans
)
: base($"""
Matrix dimensions must agree but current:
right({right?.ColumnLength ?? 0}×{right?.RowLength ?? 0}),
ans({ans?.ColumnLength ?? 0}×{ans?.RowLength ?? 0})
"""
) { }
public MatrixDimensionMismatchException(
Matrix<T> left,
Matrix<T> right,
Matrix<T> ans
)
: base($"""
Matrix dimensions must agree but current:
left({left?.ColumnLength ?? 0}×{left?.RowLength ?? 0}),
right({right?.ColumnLength ?? 0}×{right?.RowLength ?? 0}),
ans({ans?.ColumnLength ?? 0}×{ans?.RowLength ?? 0})
"""
) { }
}
}
StgSharp.Mathematics/Numeric/Linear/Matrix.cs
+62 -41
View File
@@ -2,8 +2,8 @@
// -----------------------------------------------------------------------
// file="Matrix"
// Project: StgSharp
// AuthorGroup: Nitload Space
// Copyright (c) Nitload Space. All rights reserved.
// AuthorGroup: Nitload
// Copyright (c) Nitload. All rights reserved.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
@@ -29,17 +29,17 @@ using StgSharp.Mathematics;
using StgSharp.Mathematics.Numeric;
using System;
using System.Numerics;
using System.Runtime.CompilerServices;
// using hlsfHandle = global::StgSharp.HighPerformance.Memory.HybridLayerSegregatedFitAllocationHandle;
using System.Runtime.CompilerServices;
// using hlsfHandle = global::StgSharp.HighPerformance.Memory.HybridLayerSegregatedFitAllocationHandle;
namespace StgSharp.Mathematics.Numeric
{
public unsafe partial class Matrix<T> where T: unmanaged, INumber<T>
public unsafe partial class Matrix<T> where T : unmanaged, INumber<T>
{
private protected static readonly int t_count = 16 / sizeof(T);
private protected static readonly int t_count = 16 / sizeof(T);
private protected static readonly int t_size = sizeof(T);
private readonly MatrixKernelBoarder _boarder;
private readonly MatrixSize _size;
@@ -48,8 +48,8 @@ namespace StgSharp.Mathematics.Numeric
private protected Matrix(
MatrixKernelBoarder boarder,
MatrixSize size,
MatrixStorage<T> buffer
MatrixStorage<T> buffer
)
{
_boarder = boarder;
@@ -57,7 +57,10 @@ namespace StgSharp.Mathematics.Numeric
_buffer = buffer;
}
public T this[int x, int y]
public T this[
int x,
int y
]
{
get
{
@@ -78,40 +81,43 @@ namespace StgSharp.Mathematics.Numeric
long offset = GetElementOffset(x, y);
if (offset < 0) {
return;
}
}
((T*)Buffer)[offset] = value;
}
}
}
/// <summary>
/// Matlab style indexer to generate a new matrix
/// </summary>
/// </summary>
/// <param name="column">
///
/// </param>
/// </param>
/// <param name="row">
///
/// </param>
/// </param>
/// <returns>
///
/// </returns>
/// </returns>
/// <exception cref="NotImplementedException">
///
/// </exception>
public unsafe Matrix<T> this[ReadOnlySpan<MatrixIndex> column, ReadOnlySpan<MatrixIndex> row]
/// </exception>
public unsafe Matrix<T> this[
ReadOnlySpan<MatrixIndex> column,
ReadOnlySpan<MatrixIndex> row
]
{
get { throw new NotImplementedException(); }
set { throw new NotImplementedException(); }
}
}
/// <summary>
/// Length of the matrix in columns.
/// </summary>
public int ColumnLength => _size.ElementColumnLength;
/// </summary>
public int ColumnLength => _size.ElementColumnLength;
/// <summary>
/// Length of the matrix in rows.
/// </summary>
/// </summary>
public int RowLength => _size.ElementRowLength;
public MatrixLayout Layout { get; init; }
@@ -127,49 +133,64 @@ namespace StgSharp.Mathematics.Numeric
protected internal int KernelRowLength => _size.KernelRowLength;
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public unsafe ref MatrixKernel<T> GetMatrixKernel(int x, int y)
public unsafe ref MatrixKernel<T> GetMatrixKernel(
int x,
int y
)
{
return ref _boarder.IsInBuffer(x, y, in _size) ?
ref _buffer[_boarder.GetIndexOffset(x, y, in _size)] :
ref Unsafe.NullRef<MatrixKernel<T>>();
}
internal long GetElementOffset(int x, int y)
internal long GetElementOffset(
int x,
int y
)
{
int
k_c = y / 4,
c = y % 4,
k_r = x / 4,
int
k_c = y / 4,
c = y % 4,
k_r = x / 4,
r = x % 4;
if (!_boarder.IsInBuffer(k_r, k_c, in _size)) {
return -1;
}
long kernelOffset = _boarder.GetIndexOffset(k_r, k_c, in _size);
long kernelOffset = _boarder.GetIndexOffset(k_r, k_c, in _size);
return kernelOffset * 16 + 4 * r + c;
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal unsafe ref T GetElementUnsafe(int x, int y)
internal unsafe ref T GetElementUnsafe(
int x,
int y
)
{
// Caller guarantees: indices in range, layout is dense square.
int k_c = y / 4;
int k_r = x / 4;
int c = y % 4;
int r = x % 4;
int r = x % 4;
ref MatrixKernel<T> kernel = ref _buffer[_boarder.GetIndexOffset(k_r, k_c, in _size)];
return ref kernel[c, r];
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal unsafe MatrixKernel<T>* GetKernelAddressUnsafe(int x, int y)
internal unsafe MatrixKernel<T>* GetKernelAddressUnsafe(
int x,
int y
)
{
return Buffer + _boarder.GetIndexOffset(x, y, in _size);
}
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal unsafe ref MatrixKernel<T> GetKernelUnsafe(int x, int y)
internal unsafe ref MatrixKernel<T> GetKernelUnsafe(
int x,
int y
)
{
// Caller guarantees: indices in range, layout is dense square.
int k_c = y / 4;
@@ -180,4 +201,4 @@ namespace StgSharp.Mathematics.Numeric
}
}
}
}
StgSharp.Mathematics/Numeric/Linear/MatrixCompute.Factorize.cs.old
+160
View File
@@ -0,0 +1,160 @@
//-----------------------------------------------------------------------
// -----------------------------------------------------------------------
// file="MatrixCompute.Factorize"
// Project: StgSharp
// AuthorGroup: Nitload
// Copyright (c) Nitload. All rights reserved.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//
// -----------------------------------------------------------------------
// -----------------------------------------------------------------------
using StgSharp.Internal;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace StgSharp.Mathematics.Numeric
{
public static partial class MatrixCompute
{
public static unsafe void FactorizePLU(
Matrix<float> matrix,
Span<int> pivot,
Matrix<float> lower,
Matrix<float> upper
)
{
if (matrix.ColumnLength != matrix.RowLength) {
throw new InvalidOperationException("Matrix must be square.");
}
if (matrix.Layout != MatrixLayout.DenseRectangle) {
throw new NotSupportedException("Only dense square matrices are supported.");
}
int n = matrix.RowLength;
if (pivot.Length < n) {
throw new ArgumentException("Pivot span length is insufficient.", nameof(pivot));
}
if (lower.RowLength != n ||
lower.ColumnLength != n ||
lower.Layout != MatrixLayout.UpperTriangle) {
throw new ArgumentException("Lower matrix must be lower-triangular with same dimensions as input.", nameof(lower));
}
if (upper.RowLength != n ||
upper.ColumnLength != n ||
upper.Layout != MatrixLayout.UpperTriangle) {
throw new ArgumentException("Upper matrix must be upper-triangular with same dimensions as input.", nameof(upper));
}
int size = matrix.KernelColumnLength * matrix.KernelRowLength;
if (size < 256)
{
for (int i = 0; i < n; i++)
{
pivot[i] = i; // identity permutation
}
// Aggressive: always call native pivot to fetch four-column candidates
ulong* pivots = stackalloc ulong[4];
for (int kBlock = 0; kBlock < n; kBlock += 4)
{
NumericalModule.GlobalContext
.mat_sk[(int)MatrixElementType.F32].pivot(matrix.Buffer,
(ulong)(kBlock >> 2),
(ulong)(matrix.KernelRowLength - (kBlock >> 2)),
(ulong)pivots);
int colsInBlock = Math.Min(4, n - kBlock);
for (int offset = 0; offset < colsInBlock; offset++)
{
int k = kBlock + offset;
int pivotRow = k;
float maxAbs;
int candidateRow = (int)pivots[offset];
int idx = -1;
for (int r = k; r < n; r++)
{
if (pivot[r] == candidateRow)
{
idx = r;
break;
}
}
if (idx != -1)
{
pivotRow = idx;
maxAbs = MathF.Abs(matrix.GetElementUnsafe(pivot[idx], k));
} else
{
maxAbs = MathF.Abs(matrix.GetElementUnsafe(pivot[k], k));
for (int r = k + 1; r < n; r++)
{
float val = MathF.Abs(matrix.GetElementUnsafe(pivot[r], k));
if (val > maxAbs)
{
maxAbs = val;
pivotRow = r;
}
}
}
if (maxAbs == 0f) {
throw new InvalidOperationException("Zero pivot encountered in PLU.");
}
if (pivotRow != k) {
(pivot[k], pivot[pivotRow]) = (pivot[pivotRow], pivot[k]);
}
lower[k, k] = 1.0f; // unit diagonal for L
// Compute U[k, j] for j >= k
for (int j = k; j < n; j++)
{
float sum = 0f;
for (int s = 0; s < k; s++) {
sum += lower[k, s] * upper[s, j];
}
upper[k, j] = matrix.GetElementUnsafe(pivot[k], j) - sum;
}
// Compute L[i, k] for i > k
for (int i = k + 1; i < n; i++)
{
float sum = 0f;
for (int s = 0; s < k; s++) {
sum += lower[i, s] * upper[s, k];
}
float u = upper[k, k];
if (u == 0f) {
throw new InvalidOperationException("Zero pivot encountered in PLU.");
}
lower[i, k] = (matrix.GetElementUnsafe(pivot[i], k) - sum) / u;
}
}
}
// TODO: For larger kernels, consider calling native pivot helper to extract 4 pivots at once.
} else if (size < 1024) { } else { }
}
}
}
StgSharp.Mathematics/Numeric/Linear/MatrixCompute.Float.cs.old
+314
View File
@@ -0,0 +1,314 @@
//-----------------------------------------------------------------------
// -----------------------------------------------------------------------
// file="MatrixCompute.Float"
// Project: StgSharp
// AuthorGroup: Nitload
// Copyright (c) Nitload. All rights reserved.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//
// -----------------------------------------------------------------------
// -----------------------------------------------------------------------
using StgSharp.Collections;
using StgSharp.Internal;
using StgSharp.Mathematics;
using StgSharp.Mathematics.Internal;
using System;
using System.Collections.Generic;
using System.ComponentModel.DataAnnotations.Schema;
using System.Linq;
using System.Numerics;
using System.Security.AccessControl;
using System.Text;
using System.Threading.Tasks;
namespace StgSharp.Mathematics.Numeric
{
public static partial class MatrixCompute
{
public static unsafe void Add(
Matrix<float> left,
Matrix<float> right,
Matrix<float> ans
)
{
MatrixElementType f32 = MatrixElementType.F32;
if (left.ColumnLength != right.ColumnLength ||
left.RowLength != right.RowLength ||
left.ColumnLength != ans.ColumnLength ||
left.RowLength != ans.RowLength) {
throw new IndexOutOfRangeException("Matrix dimensions must agree.");
}
long count = (long)left.KernelColumnLength * left.KernelRowLength;
if (count <= 256)
{
MatrixKernel<float>* leftPtr = left.Buffer;
MatrixKernel<float>* rightPtr = right.Buffer;
MatrixKernel<float>* ansPtr = ans.Buffer;
NumericalModule.GlobalContext.mat_mk[f32].buffer_add(leftPtr, rightPtr, ansPtr,
count);
return;
} else if (count <= 1024)
{
MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
MatrixKernel<float>* leftPtr = left.Buffer;
MatrixKernel<float>* rightPtr = right.Buffer;
MatrixKernel<float>* ansPtr = ans.Buffer;
NumericalModule.GlobalContext.mat_mk[f32].buffer_add(leftPtr, rightPtr, ansPtr,
count);
MatrixParallel.ReturnParallelResources(ref pool);
return;
}
MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
package->Mat1 = (MatrixKernel*)left.Buffer;
package->Mat2 = (MatrixKernel*)right.Buffer;
package->Mat3 = (MatrixKernel*)ans.Buffer;
package->ElementSize = sizeof(float);
package->PrimCount = ans.KernelColumnLength;
package->SecCount = ans.KernelRowLength;
package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f32].buffer_add;
MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
taskGroup.LaunchParallel(SleepMode.DeepSleep);
taskGroup.Dispose();
return;
}
public static unsafe void Fill(
Matrix<float> ans,
float value
)
{
MatrixElementType f32 = MatrixElementType.F32;
long count = (long)ans.KernelColumnLength * ans.KernelRowLength;
ScalarPacket* scalar = MatrixParallelFactory.CreateScalarPacket();
scalar->Data<float>(0) = value;
if (count <= 256)
{
NumericalModule.GlobalContext.mat_sk[f32].buffer_fill(ans.Buffer, scalar, count);
return;
} else if (count <= 1024)
{
MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
NumericalModule.GlobalContext.mat_sk[f32].buffer_fill(ans.Buffer, scalar, count);
MatrixParallel.ReturnParallelResources(ref pool);
return;
}
MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
package->Answer = (nint)ans.Buffer;
package->ElementSize = sizeof(float);
package->PrimCount = ans.KernelColumnLength;
package->SecCount = ans.KernelRowLength;
package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.ANS_SCALAR_PARAM);
package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f32].buffer_fill;
package->Scalar = scalar;
MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
taskGroup.LaunchParallel(SleepMode.DeepSleep);
MatrixParallelFactory.Release(scalar);
taskGroup.Dispose();
}
public static unsafe void Mul(
Matrix<float> left,
Matrix<float> right,
Matrix<float> ans
)
{
MatrixElementType f32 = MatrixElementType.F32;
// For matrix multiplication (left: MxK) * (right: KxN) => ans: MxN
// ColumnLength = number of columns (N), RowLength = number of rows (M)
if (left.ColumnLength != right.RowLength ||
ans.ColumnLength != right.ColumnLength ||
ans.RowLength != left.RowLength)
{
throw new IndexOutOfRangeException("Matrix dimensions must agree.");
}
if (left.Layout != MatrixLayout.DenseRectangle ||
right.Layout != MatrixLayout.DenseRectangle ||
ans.Layout != MatrixLayout.DenseRectangle) {
throw new NotSupportedException("Only dense rectangle matrices are supported.");
}
long count = (long)left.KernelColumnLength * left.KernelRowLength;
if (count <= 1024)
{
MatrixParallelThread[] pool = null;
if (count > 256) {
pool = MatrixParallel.LeadParallel(1);
}
M128* enumeration = stackalloc M128[1];
enumeration->Member<int>(2) = left.KernelColumnLength;
enumeration->Member<int>(3) = ans.KernelColumnLength;
for (int i = 0; i < ans.KernelRowLength; i++)
{
for (int j = 0; j < ans.KernelColumnLength; j++)
{
enumeration->Member<int>(0) = i;
enumeration->Member<int>(1) = j;
NumericalModule.GlobalContext
.mat_sk[f32].kernel_tile_fma((MatrixKernel*)left.Buffer,
(MatrixKernel*)right.Buffer,
(MatrixKernel*)ans.Buffer,
enumeration);
}
}
if (count > 256) {
MatrixParallel.ReturnParallelResources(ref pool!);
}
return;
}
MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
package->Left = (nint)left.Buffer;
package->Right = (nint)right.Buffer;
package->Answer = (nint)ans.Buffer;
package->AnsPrimOffset = ans.KernelRowLength;
package->AnsSecOffset = ans.KernelColumnLength;
package->LeftPrimOffset = left.KernelColumnLength;
package->PrimCount = ans.KernelColumnLength;
package->SecCount = ans.KernelRowLength;
package->ElementSize = (int)MatrixElementType.F32;
package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.SPECIAL, MatrixOperationSpecial.SEQ_FMA);
MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
taskGroup.LaunchParallel(SleepMode.DeepSleep);
taskGroup.Dispose();
return;
}
public static unsafe void Mul(
Matrix<float> left,
float right,
Matrix<float> ans
)
{
MatrixElementType f32 = MatrixElementType.F32;
if (left.ColumnLength != ans.ColumnLength || left.RowLength != ans.RowLength) {
throw new IndexOutOfRangeException("Matrix dimensions must agree.");
}
if (left.Layout != MatrixLayout.DenseRectangle ||
ans.Layout != MatrixLayout.DenseRectangle) {
throw new NotSupportedException("Only dense rectangle matrices are supported.");
}
long count = (long)left.KernelColumnLength * left.KernelRowLength;
ScalarPacket* scalar = MatrixParallelFactory.CreateScalarPacket();
scalar->Data<float>(0) = right;
if (count <= 256)
{
MatrixKernel<float>* leftPtr = left.Buffer;
MatrixKernel<float>* ansPtr = ans.Buffer;
NumericalModule.GlobalContext.mat_sk[f32].buffer_scalar_mul(leftPtr, ansPtr, scalar,
count);
return;
} else if (count <= 1024)
{
MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
MatrixKernel<float>* leftPtr = left.Buffer;
MatrixKernel<float>* ansPtr = ans.Buffer;
NumericalModule.GlobalContext.mat_sk[f32].buffer_scalar_mul(leftPtr, ansPtr, scalar,
count);
MatrixParallel.ReturnParallelResources(ref pool);
return;
}
MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
package->Right = (nint)left.Buffer;
package->Answer = (nint)ans.Buffer;
package->ElementSize = sizeof(float);
package->PrimCount = left.KernelColumnLength;
package->SecCount = left.KernelRowLength;
package->Scalar = scalar;
package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.RIGHT_ANS_SCALAR_PARAM);
package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f32].buffer_scalar_mul;
MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
taskGroup.LaunchParallel(SleepMode.DeepSleep);
taskGroup.Dispose();
MatrixParallelFactory.Release(scalar);
return;
}
public static unsafe void Sub(
Matrix<float> left,
Matrix<float> right,
Matrix<float> ans
)
{
MatrixElementType f32 = MatrixElementType.F32;
if (left.ColumnLength != right.ColumnLength ||
left.RowLength != right.RowLength ||
left.ColumnLength != ans.ColumnLength ||
left.RowLength != ans.RowLength) {
throw new IndexOutOfRangeException("Matrix dimensions must agree.");
}
if (left.Layout != MatrixLayout.DenseRectangle ||
right.Layout != MatrixLayout.DenseRectangle ||
ans.Layout != MatrixLayout.DenseRectangle) {
throw new NotSupportedException("Only dense rectangle matrices are supported.");
}
long count = (long)left.KernelColumnLength * left.KernelRowLength;
if (count <= 256)
{
MatrixKernel<float>* leftPtr = left.Buffer;
MatrixKernel<float>* rightPtr = right.Buffer;
MatrixKernel<float>* ansPtr = ans.Buffer;
NumericalModule.GlobalContext.mat_sk[f32].buffer_sub(leftPtr, rightPtr, ansPtr,
count);
return;
} else if (count <= 1024)
{
MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
MatrixKernel<float>* leftPtr = left.Buffer;
MatrixKernel<float>* rightPtr = right.Buffer;
MatrixKernel<float>* ansPtr = ans.Buffer;
NumericalModule.GlobalContext.mat_sk[f32].buffer_sub(leftPtr, rightPtr, ansPtr,
count);
MatrixParallel.ReturnParallelResources(ref pool);
return;
}
MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
package->Left = (nint)left.Buffer;
package->Right = (nint)right.Buffer;
package->Answer = (nint)ans.Buffer;
package->ElementSize = sizeof(float);
package->PrimCount = left.KernelColumnLength;
package->SecCount = left.KernelRowLength;
package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f32].buffer_sub;
MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
taskGroup.LaunchParallel(SleepMode.DeepSleep);
taskGroup.Dispose();
return;
}
}
}
StgSharp.Mathematics/Numeric/Linear/MatrixCompute.cs
+55
View File
@@ -0,0 +1,55 @@
//-----------------------------------------------------------------------
// -----------------------------------------------------------------------
// file="MatrixCompute"
// Project: StgSharp
// AuthorGroup: Nitload Space
// Copyright (c) Nitload Space. All rights reserved.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//
// -----------------------------------------------------------------------
// -----------------------------------------------------------------------
using StgSharp.Internal;
using StgSharp.Mathematics;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Text;
using System.Threading.Tasks;
namespace StgSharp.Mathematics.Numeric
{
public partial class MatrixCompute
{
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static unsafe MatrixKernel* GetKernelAddressUnsafe<T>(
MatrixKernel<T>* buffer,
int colLength,
int x,
int y)
where T: unmanaged, INumber<T>
{
return (MatrixKernel*)(buffer + (colLength * x + y));
}
}
}
StgSharp.Mathematics/Numeric/Linear/MatrixCompute.double.cs.old
+320
View File
@@ -0,0 +1,320 @@
//-----------------------------------------------------------------------
// -----------------------------------------------------------------------
// file="MatrixCompute.double"
// Project: StgSharp
// AuthorGroup: Nitload
// Copyright (c) Nitload. All rights reserved.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//
// -----------------------------------------------------------------------
// -----------------------------------------------------------------------
using StgSharp.Internal;
using StgSharp.Mathematics;
using StgSharp.Mathematics.Internal;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using unsafe MkProc = delegate* unmanaged[Cdecl]<StgSharp.Mathematics.Internal.MatrixParallelTask, void>;
namespace StgSharp.Mathematics.Numeric
{
public partial class MatrixCompute
{
public static unsafe void Add(
Matrix<double> left,
Matrix<double> right,
Matrix<double> ans
)
{
MatrixElementType f64 = MatrixElementType.F64;
if (left.ColumnLength != right.ColumnLength ||
left.RowLength != right.RowLength ||
left.ColumnLength != ans.ColumnLength ||
left.RowLength != ans.RowLength) {
throw new MatrixDimensionMismatchException<double>(
left, right, ans
);
}
long count = (long)left.KernelColumnLength * left.KernelRowLength;
if (count <= 256)
{
MatrixKernel<double>* leftPtr = left.Buffer;
MatrixKernel<double>* rightPtr = right.Buffer;
MatrixKernel<double>* ansPtr = ans.Buffer;
MkProc add = (MkProc)NumericalModule.GlobalContext.
mat_mk[f64.IntrinsicNode].
Get(MatrixIntrinsicHandle.Add);
for (int i = 0; i < count; i++) {
add(new MatrixParallelTask
{
Mat1 = (MatrixKernel*)(leftPtr + i),
Mat2 = (MatrixKernel*)(rightPtr + i),
Mat3 = (MatrixKernel*)(ansPtr + i),
ElementType = MatrixElementType.F64,
});
}
// .buffer_add(leftPtr, rightPtr, ansPtr, count);
return;
} else if (count <= 1024)
{
MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
MatrixKernel<double>* leftPtr = left.Buffer;
MatrixKernel<double>* rightPtr = right.Buffer;
MatrixKernel<double>* ansPtr = ans.Buffer;
NumericalModule.GlobalContext.mat_sk[f64].buffer_add(leftPtr, rightPtr, ansPtr,
count);
MatrixParallel.ReturnParallelResources(ref pool);
return;
}
MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
package->Left = (nint)left.Buffer;
package->Right = (nint)right.Buffer;
package->Answer = (nint)ans.Buffer;
package->ElementSize = sizeof(double);
package->PrimCount = ans.KernelColumnLength;
package->SecCount = ans.KernelRowLength;
package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f64].buffer_add;
MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
taskGroup.LaunchParallel(SleepMode.DeepSleep);
taskGroup.Dispose();
return;
}
public static unsafe void Fill(
Matrix<double> ans,
double value
)
{
int f64 = (int)MatrixElementType.F64;
long count = (long)ans.KernelColumnLength * ans.KernelRowLength;
ScalarPacket* scalar = MatrixParallelFactory.CreateScalarPacket();
scalar->Data<double>(0) = value;
if (count <= 256)
{
NumericalModule.GlobalContext.mat_sk[f64].buffer_fill(ans.Buffer, scalar, count);
return;
} else if (count <= 1024)
{
MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
NumericalModule.GlobalContext.mat_sk[f64].buffer_fill(ans.Buffer, scalar, count);
MatrixParallel.ReturnParallelResources(ref pool);
return;
}
MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
package->Answer = (nint)ans.Buffer;
package->ElementSize = sizeof(double);
package->PrimCount = ans.KernelColumnLength;
package->SecCount = ans.KernelRowLength;
package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.ANS_SCALAR_PARAM);
package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f64].buffer_fill;
package->Scalar = scalar;
MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
taskGroup.LaunchParallel(SleepMode.DeepSleep);
MatrixParallelFactory.Release(scalar);
taskGroup.Dispose();
}
public static unsafe void Mul(
Matrix<double> left,
Matrix<double> right,
Matrix<double> ans
)
{
int f64 = (int)MatrixElementType.F64;
// For matrix multiplication (left: MxK) * (right: KxN) => ans: MxN
// ColumnLength = number of columns (N), RowLength = number of rows (M)
if (left.ColumnLength != right.RowLength ||
ans.ColumnLength != right.ColumnLength ||
ans.RowLength != left.RowLength)
{
throw new IndexOutOfRangeException("Matrix dimensions must agree.");
}
long count = (long)left.KernelColumnLength * left.KernelRowLength;
if (count <= 1024)
{
MatrixParallelThread[] pool = null;
if (count > 256) {
pool = MatrixParallel.LeadParallel(1);
}
M128* enumeration = stackalloc M128[1];
enumeration->Member<int>(2) = left.KernelColumnLength;
enumeration->Member<int>(3) = ans.KernelColumnLength;
for (int i = 0; i < ans.KernelRowLength; i++)
{
for (int j = 0; j < ans.KernelColumnLength; j++)
{
enumeration->Member<int>(0) = i;
enumeration->Member<int>(1) = j;
NumericalModule.GlobalContext
.mat_sk[f64].kernel_tile_fma((MatrixKernel*)left.Buffer,
(MatrixKernel*)right.Buffer,
(MatrixKernel*)ans.Buffer,
enumeration);
}
}
if (count > 256) {
MatrixParallel.ReturnParallelResources(ref pool!);
}
return;
}
MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
package->Left = (nint)left.Buffer;
package->Right = (nint)right.Buffer;
package->Answer = (nint)ans.Buffer;
package->AnsPrimOffset = ans.KernelRowLength;
package->AnsSecOffset = ans.KernelColumnLength;
package->LeftPrimOffset = left.KernelColumnLength;
package->PrimCount = ans.KernelColumnLength;
package->SecCount = ans.KernelRowLength;
package->ElementSize = (int)MatrixElementType.F64;
package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.SPECIAL, MatrixOperationSpecial.SEQ_FMA);
MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
taskGroup.LaunchParallel(SleepMode.DeepSleep);
taskGroup.Dispose();
return;
}
public static unsafe void Mul(
Matrix<double> left,
double right,
Matrix<double> ans
)
{
int f64 = (int)MatrixElementType.F64;
if (left.ColumnLength != ans.ColumnLength || left.RowLength != ans.RowLength) {
throw new IndexOutOfRangeException("Matrix dimensions must agree.");
}
if (left.Layout != MatrixLayout.DenseRectangle ||
ans.Layout != MatrixLayout.DenseRectangle) {
throw new NotSupportedException("Only dense rectangle matrices are supported.");
}
long count = (long)left.KernelColumnLength * left.KernelRowLength;
ScalarPacket* scalar = MatrixParallelFactory.CreateScalarPacket();
scalar->Data<double>(0) = right;
if (count <= 256)
{
MatrixKernel<double>* leftPtr = left.Buffer;
MatrixKernel<double>* ansPtr = ans.Buffer;
NumericalModule.GlobalContext.mat_sk[f64].buffer_scalar_mul(leftPtr, ansPtr, scalar,
count);
return;
} else if (count <= 1024)
{
MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
MatrixKernel<double>* leftPtr = left.Buffer;
MatrixKernel<double>* ansPtr = ans.Buffer;
NumericalModule.GlobalContext.mat_sk[f64].buffer_scalar_mul(leftPtr, ansPtr, scalar,
count);
MatrixParallel.ReturnParallelResources(ref pool);
return;
}
MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
package->Right = (nint)left.Buffer;
package->Answer = (nint)ans.Buffer;
package->ElementSize = sizeof(double);
package->PrimCount = left.KernelColumnLength;
package->SecCount = left.KernelRowLength;
package->Scalar = scalar;
package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.RIGHT_ANS_SCALAR_PARAM);
package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f64].buffer_scalar_mul;
MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
taskGroup.LaunchParallel(SleepMode.DeepSleep);
taskGroup.Dispose();
MatrixParallelFactory.Release(scalar);
return;
}
public static unsafe void Sub(
Matrix<double> left,
Matrix<double> right,
Matrix<double> ans
)
{
int f64 = (int)MatrixElementType.F64;
if (left.ColumnLength != right.ColumnLength ||
left.RowLength != right.RowLength ||
left.ColumnLength != ans.ColumnLength ||
left.RowLength != ans.RowLength) {
throw new IndexOutOfRangeException("Matrix dimensions must agree.");
}
if (left.Layout != MatrixLayout.DenseRectangle ||
right.Layout != MatrixLayout.DenseRectangle ||
ans.Layout != MatrixLayout.DenseRectangle) {
throw new NotSupportedException("Only dense rectangle matrices are supported.");
}
long count = (long)left.KernelColumnLength * left.KernelRowLength;
if (count <= 256)
{
MatrixKernel<double>* leftPtr = left.Buffer;
MatrixKernel<double>* rightPtr = right.Buffer;
MatrixKernel<double>* ansPtr = ans.Buffer;
NumericalModule.GlobalContext.mat_sk[f64].buffer_sub(leftPtr, rightPtr, ansPtr,
count);
return;
} else if (count <= 1024)
{
MatrixParallelThread[] pool = MatrixParallel.LeadParallel(1);
MatrixKernel<double>* leftPtr = left.Buffer;
MatrixKernel<double>* rightPtr = right.Buffer;
MatrixKernel<double>* ansPtr = ans.Buffer;
NumericalModule.GlobalContext.mat_sk[f64].buffer_sub(leftPtr, rightPtr, ansPtr,
count);
MatrixParallel.ReturnParallelResources(ref pool);
return;
}
MatrixParallelTask* package = MatrixParallelFactory.CreateBaseTask();
package->Left = (nint)left.Buffer;
package->Right = (nint)right.Buffer;
package->Answer = (nint)ans.Buffer;
package->ElementSize = sizeof(double);
package->PrimCount = left.KernelColumnLength;
package->SecCount = left.KernelRowLength;
package->ComputeMode = new MatrixOpMode(MatrixIndexStyle.BUFFER_INDEX, MatrixOperationParam.LEFT_RIGHT_ANS_PARAM);
package->ComputeHandle = (IntPtr)NumericalModule.GlobalContext.mat_sk[f64].buffer_sub;
MatrixParallelWrap taskGroup = MatrixParallel.ScheduleTask(package);
taskGroup.LaunchParallel(SleepMode.DeepSleep);
taskGroup.Dispose();
return;
}
}
}