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

keep working on PLN facotrize

Browse Source
This commit is contained in:
Youmu Liu
2025-12-24 18:09:16 +08:00
parent b60ce56f1c
commit 3b3a90a2af
12 changed files with 428 additions and 305 deletions
Download Patch File Download Diff File Expand all files Collapse all files
StgSharp.Mathematics/Numeric/MatrixCompute.Factorize.cs
+54 -2
View File
@@ -38,17 +38,69 @@ namespace StgSharp.Mathematics.Numeric
public static void FactorizePLU(
Matrix<float> matrix,
Span<float> p,
Span<float> 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.LowerTriangle) {
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 < matrix.KernelColumnLength; i++) { }
// Naive Doolittle-style PLU without row swaps (for illustration).
for (int i = 0; i < n; i++)
{
pivot[i] = i; // identity pivot
lower[i, i] = 1.0f; // unit diagonal for L
}
for (int k = 0; k < n; k++)
{
// 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[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 demonstration PLU.");
}
lower[i, k] = (matrix.GetElementUnsafe(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/MatrixCompute.cs
+9
View File
@@ -200,6 +200,10 @@ namespace StgSharp.Mathematics.Numeric
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;
@@ -246,6 +250,11 @@ namespace StgSharp.Mathematics.Numeric
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)
{