wo80 / csparse.net Goto Github PK
View Code? Open in Web Editor NEWA concise library for solving sparse linear systems with direct methods.
License: GNU Lesser General Public License v2.1
A concise library for solving sparse linear systems with direct methods.
License: GNU Lesser General Public License v2.1
Inside Dulmage–Mendelsohn decomposition file there is a section for finding Strongly Connected Components (SCC). Can this code be used to find SCC of a given sparse matrix? how to interpret result?
i need to go through all edges and identify only edges that do connect to different components. eg. in below image:
need to identify that bc, bf, ef, cg and hg edges are connection between different components.
Is it possible that your code be used for such purpose?
---- edit
All i need is to know how to assign same number to nodes in a single SCC?
Thanks
Hi there.
How can someone use CSparse.NET in a commercial closed source software?
Does it need re-licensing? or pay for new license?
Thanks
Hi, first of all thanks for a great library, you are really saving my thesis =). I was wondering is there some parts of unmanaged code in CSparse or it is a direct translation to C# of Tim Davis`s original code? If there are native code parts, can you describe them a bit more in detail (which part exactly, and is it made in native Fortran C/C++) Thanks again.
Hi,
Thank you very much for publishing this library. Really great to have an open source sparse solver in .NET.
I've got just one question.
I encountered an error with feeding a rectangular matrix defined in MathNet.Numerics to CSparse for QR decomposition.
I followed the tutorial from this page (https://github.com/wo80/CSparse.NET/wiki/Math.NET-Numerics-and-CSparse) but didn't manage to figure out the cause. It seems to work fine for square matrices.
Have anyone made it work?
If it is supported, could you elaborate the procedure?
Best,
Kenryo
I require an inverse of a sparse matrix. The LU implementation is very suitable to do this operation, but it lacks the final step.
In fact the inverse is solving matrix M with identity matrix as right-hand side (matrix instead of rhz vector)
Hello,
Here's a test case
[Test]
public void TestSparseCholesky()
{
double[,] A = new double[,]
{
{1.0, 0.0, 4.0, 0.0 , 0.2 },
{0.0, 1.0, 6.0, 0.0 , 0.0 },
{4.0, 6.0, 6.0, 3.0 , 3.0 },
{0.0, 0.0, 3.0, 0.5 , 0.0 },
{0.2, 0.0, 3.0, 0.0 , 0.5 }
};
var x = new double[] { 1, 2, 3, 4, 5 };
var y = new double[5];
var ccs = Converter.ToCompressedColumnStorage(Converter.FromDenseArray(A));
var chol = new CSparse.Double.Factorization.SparseCholesky(ccs, ColumnOrdering.MinimumDegreeAtPlusA);
chol.Solve(x, y);
}
Hi,
You have any existing code that be able to use Intel's MKL or CUDA for cholesky factorization of a CompressedColumnStorage ?
I have created a list of triplets. But I don’t know if I can convert them into a CSparse matrix. Besides, I very like your work but where can I get some documentation or tutorials?
Hi,
Is there code for LDLt decomposition of symmetric matrices with the library?
Hi!
I have several issues, most of which stem from the first one here (let SparseMatrix
be the one from Math.NET and CSparseMatrix
the one from Csparse.NET):
CSparseMatrix
yields incorrect results.a
var a = (SparseMatrix)SparseMatrix.Build.SparseOfRowArrays([[1, 2], [3, 4], [5, 6]]);
var aStorage = (SparseCompressedRowMatrixStorage<double>)a.Storage;
var aSparse = new CSparseMatrix(a.RowCount, a.ColumnCount)
{
// Based on https://github.com/wo80/CSparse.NET/wiki/Math.NET-Numerics-and-CSparse.
RowIndices = aStorage.ColumnIndices,
ColumnPointers = aStorage.RowPointers,
Values = aStorage.Values
};
var aDash = aSparse.Transpose();
aDash
to be
1 3 5
2 4 6
1 2 0
3 4 0
%%MatrixMarket matrix coordinate real general
2 3 4
1 1 1
2 1 3
1 2 2
2 2 4
MatrixMarketWriter.WriteMatrix
).aSparse
is incorrectly dumped as
%%MatrixMarket matrix coordinate real general
3 2 4
1 1 1
2 1 2
1 2 3
2 2 4
Values
property is [1,2,3,4,5,6]
(with 3 rows and 2 columns).A = B*C
equivalent to var A = B.ParallelMultiply(C)
or var A = C.ParallelMultiply(B)
?var aTilde = aDash.ParallelMultiply(aSparse);
I get
5 11
11 25
var aTilde = aSparse.ParallelMultiply(aDash);
I get
10 14 0
14 20 0
0 0 0
aDash*aSparse
and aSparse*aDash
should be
10 14
14 20
5 11 0
11 25 0
0 0 0
aDash
and aSparse
that were computed by the library (i.e. the unexpected ones).Solve
rs.var solver = SparseQR.Create(aTilde, ColumnOrdering.Natural);
Solve
(with ColumnOrdering.Natural
) yields the wrong result.I want to numerically solve the regularized l_2-minimization problem|| Ax -b ||_2 + λ|| x ||_2
, where A
is some noise operator, b
a noisy signal, x
the denoised signal that I want to find, and λ > 0
is an arbitrary regularization parameter.
´x´ and ´b´ have the same number of entries.
This is equivalent to an l_2-minimization problem of the form || Cx_c - y ||_2
, with C in |R^mxn
being a sparse matrix dependent on λ
with m > n
.
According to the documentation of SparseQR.Solve
,
If m >= n a least-squares problem (min |Ax-b|) is solved.
which is exactly what I want there.
Unfortunately, this returns wrong values for my actual problem that I want to solve, where m
is about 8000 and n
about 4000.
So I heavily simplified the problem, which yields the results above.
When multiplying two sparse matrices that are all zeros, the following exception occurs:
Unhandled Exception:
System.NullReferenceException: Object reference not set to an instance of an object
at CSparse.Double.SparseMatrix.Multiply (CSparse.Storage.CompressedColumnStorage`1[T] other) [0x000ae] in <441d951ef3a340c28331a195ffb93d84>:0
at EmptyMultiply.Main (System.String[] args) [0x0001c] in <265f17d61719442385620412a8f37a5e>:0
[ERROR] FATAL UNHANDLED EXCEPTION: System.NullReferenceException: Object reference not set to an instance of an object
at CSparse.Double.SparseMatrix.Multiply (CSparse.Storage.CompressedColumnStorage`1[T] other) [0x000ae] in <441d951ef3a340c28331a195ffb93d84>:0
at EmptyMultiply.Main (System.String[] args) [0x0001c] in <265f17d61719442385620412a8f37a5e>:0
This is due to the fact that Values
isn't allocated for empty matrices: https://github.com/wo80/CSparse.NET/blob/master/CSparse/Storage/CompressedColumnStorage.cs#L57
result
is an empty matrix if both inputs are empty: https://github.com/wo80/CSparse.NET/blob/master/CSparse/Double/SparseMatrix.cs#L393
However, Resize
doesn't check if Values
is null
: https://github.com/wo80/CSparse.NET/blob/master/CSparse/Storage/CompressedColumnStorage.cs#L762
Sample program:
using System;
using CSparse;
using CSparse.Double;
public class EmptyMultiply
{
public static void Main(string[] args)
{
var a = SparseMatrix.OfColumnMajor(3, 2, new double[6]);
var b = SparseMatrix.OfColumnMajor(2, 4, new double[8]);
a.Multiply(b);
}
}
Hi,
It's me again. Sorry for interrupting.
I'm working on a Finite element system based on Unity. I used CSparse.NET and Math.NET. I found that the bottleneck is matrix processing and computing, such as Multiply(), Get_Item() and CompressedColumnStorage.OfIndexed(). It means that matrix editing and computing is much slower than solving the linear system.
The Multiply() and Get_Item() I used is Math.NET. But how can I implement them in CSparse? I want to make it as fast as my Matlab program.
Best wishes
Hello,
I'm writing a Finite element framework. And I need to remove several rows and columns from a CSC sparse matrix A
. In Matlab, I can use slicing operations like A([2, 3, 4], [6, 7, 8])
to select several rows and columns for computing. How can I implement this using CSparse?
Or can I access a value from the CSC sparse matrix by inputting (i, j)
?
Cheers
Hello,
Is there any straight way to compute a "Reduced row echelon form" of a n*m
matrix in CSparse.Net?
Definition of Reduced row echelon form in wikipedia.
Thanks
To start: Awesome project. Would AVX vectorization with System.Numerics.Vector be possible? Is the original CSparse making use of SIMD?
In SparseMatrix.Scatter()
move the null
check out of the for loop:
internal override int Scatter(int j, double beta, int[] w, double[] x, int mark,
CompressedColumnStorage<double> mat, int nz)
{
int i, p;
if (w == null || mat == null) return -1; // check inputs
var cj = mat.RowIndices;
var ap = ColumnPointers;
var ai = RowIndices;
var ax = Values;
if (x != null)
{
for (p = ap[j]; p < ap[j + 1]; p++)
{
i = ai[p]; // A(i,j) is nonzero
if (w[i] < mark)
{
w[i] = mark; // i is new entry in column j
x[i] = beta * ax[p]; // x(i) = beta*A(i,j)
cj[nz++] = i; // add i to pattern of C(:,j)
}
else
{
x[i] += beta * ax[p]; // i exists in C(:,j) already
}
}
}
else
{
for (p = ap[j]; p < ap[j + 1]; p++)
{
i = ai[p]; // A(i,j) is nonzero
if (w[i] < mark)
{
w[i] = mark; // i is new entry in column j
cj[nz++] = i; // add i to pattern of C(:,j)
}
}
}
return nz;
}
Is it currently foreseen in the library that a value in the SparseMatrix will be changed after first initialization?
I can see we are able to retrieve the value of sparseMatrix via .At(int i, int j). But I don't see a way to change the value.
The use case is that I need to increase the value of a diagonal cell with a penalty value after some sparsematrix addition and multiplication has taken place.
Firstly, thanks for your wonderful effort.
My question is if it possible to include matrix multiplication in CSparse.NET.
I currently have a loop using CSparse.NET to solve equations and Mathnet.Numerics to do multiplication. So i am having to convert back and forth between the the different formats of MathNet Vector and CSparse.Net double().
I don't know if there will be any performance improvements. This is my code.
Dim KeGlobalFree As SparseMatrix = MathNet.Numerics.LinearAlgebra.Double.SparseMatrix.Build.Sparse(FreeDOFs.Count, FreeDOFs.Count)
Dim KgGlobalFree As SparseMatrix = MathNet.Numerics.LinearAlgebra.Double.SparseMatrix.Build.Sparse(FreeDOFs.Count, FreeDOFs.Count)
........Missing Code
Dim EigenVectorEstimate = MathNet.Numerics.LinearAlgebra.Double.SparseVector.Build.Dense(KgGlobalFree.RowCount, 1)
Dim Errors As Double = 10 ' = MathNet.Numerics.LinearAlgebra.Double.SparseVector.Build.Dense(KgGlobalFree.RowCount, 1)
Dim EigenValueEstimate As Double = 1.0
Dim y2 = CSparse.Double.Vector.Create(KeGlobalFree.ColumnCount, 0)
While Errors >= 10 ^ -8
Dim y1 = -KgGlobalFree * EigenVectorEstimate
Dim storage = DirectCast(KeGlobalFree.Storage, MathNet.Numerics.LinearAlgebra.Storage.SparseCompressedRowMatrixStorage(Of Double)) ' Get CSR storage.
Dim A = New CSparse.Double.SparseMatrix(KeGlobalFree.ColumnCount, KeGlobalFree.ColumnCount) With {.ColumnPointers = storage.RowPointers, .RowIndices = storage.ColumnIndices, .Values = storage.Values} ' Create CSparse matrix and Assign storage arrays.
Dim SC = CSparse.Double.Factorization.SparseCholesky.Create(A, CSparse.ColumnOrdering.MinimumDegreeAtPlusA) ' Apply column ordering to A to reduce fill-in.
SC.Solve(y1.ToArray, y2)
Dim y2vec = MathNet.Numerics.LinearAlgebra.Double.SparseVector.Build.DenseOfArray(y2)
Dim norm = y2vec * -KgGlobalFree * y2vec
EigenVectorEstimate = y2vec.Divide(Math.Sqrt(norm))
Dim NewEigenValueEstimate = EigenVectorEstimate * KeGlobalFree * EigenVectorEstimate
Errors = Math.Abs((NewEigenValueEstimate - EigenValueEstimate) / EigenValueEstimate)
EigenValueEstimate = NewEigenValueEstimate
End While
Hi,
Thanks for great code. After that i've had problem with a particular matrix and cholesky solver.
For that matrix cholesky solver gave me an array of double.NaN
. i am not sure if my matrix was positive definite, but i was expecting that if it is not SPD, then i get an matrix should be positive definite
error. is it true?
Thanks
Greetings,
why adding a matrix B to an existing matrix A requires matrix C,
I think the default behavior should be adding B to A and return A itself.
/// <summary>
/// Adds two matrices, A = alpha*A + beta*B, where A is current instance.
/// </summary>
public abstract void Add(T alpha, T beta, CompressedColumnStorage<T> other);
CSparse.NET/CSparse/Storage/CompressedColumnStorage.cs
Lines 406 to 419 in 1f1357a
Hello,
For making the CSparse.NET CLSCompliant, these minor changes are needed:
Matrix.rowCount
and Matrix.colCount
accesors from protected
into internal
Converter.FromDenseArray<T>(T[][] array)
or FromDenseArray<T>(T[,] array)
in order to have different names.Do you agree to do these changes? and then from other .net languages like VB.Net, CSparse.NET is accessible
although we already have 4 good solvers, i think we still need UMFPACK, which has the best performance for large square non-symmetric matrices,
e.g. solving a 1 000 000 x 1 000 000 sparse matrix with ~900 000 non-zero elements using Cholesky with A+A` column ordering takes about 3 minutes and about 2 gb of ram for just the factorization.
but with umfpack it takes only ~50 seconds
Hello,
I'm trying to use your library to solve a large sparse matrix. Since I'm experimenting a performance issue, I would like to know if there is a fix or workaround to resolve it.
Particularly, I have noticed that almost all the CPU time is consumed by method GraphHelper.DepthFirstSearch (see the attachment for details).
I'm wondering if there is an optimization problem or perhaps the algorithm is already optimized and there is nothing to do...?
Thank you very much,
Mac
When I create a coo matrix using pre-defined arrays, the created coo matrix doesn't have the number of non-zeros. Therefore, when I convert the COO matrix into a CSC matrix, all values are lost. Because the converting code will check the number of non-zeros. If the number of non-zeros is 0, the converting would be failed.
For example, ik, jk, sk are pre-defined arrays. The created COO matrix doesn't have the number of non-zeros. Then, matrix A has no any values.
var coo = new CoordinateStorage<double>(2 * (nelx + 1) * (nely + 1), 2 * (nelx + 1) * (nely + 1), ik, jk, sK);
var A = (SparseMatrix)SparseMatrix.OfIndexed(coo);
I can get the value using matrix.At(row, col)
, but I can't find a way to set the value.
The user should be able to configure the degree of parallelism in the parallel matrix multiplication algorithms by passing ParallelOptions to the Parallel.For
calls. This could either be done by adding a method parameter or by adding some configuartion options to the library.
The latter is done in MathNet (see Control.cs) and I'd prefer to do it the same way.
Hello! Thank you sincerely for providing a nice sparse matrix library.
I would like to apply Span and Memory to function arguments provided by CSparse.NET.
For example, SparseCholesky.Solve(ReadOnlySpan input, Span result) for SparseCholesky.Solve(double[] input, double[] result) Is it possible to add an overload on the CSparse.NET side so that I can do this?
Hi!
Would be great for debugging if you published those data with your NuGet package as described here.
Without these, I was not able to generate any usable debugging information even when decompiling.
Hi @wo80 ,
Would it be possible to provide public access for the Lower Matrix L in the SparseCholesky Factorization?
I am trying to implement the below method, so that would be of great help.
Regards,
Anthony
How can I create a matrix from multiple diagonals, just like the behavior of scipy.sparse.diags?
scipy.sparse.diags([a,b,c], [-1,0,1], ...)
Thank you for your help.
Hi, do you know any solver for nonlinear eq. system with enhanced methods like arc-length or work control specially for Finite Element Analysis of solids and structures with C#?
Thanks
Hi,
Using CoordinatedStorage like this:
using Coord = CSparse.Storage.CoordinateStorage<double>;
var buf = new Coord(10, 10, 0);
buf.At(2,3,1.0);
will throw exception of IndexOutOfRangeException
, but this will not:
using Coord = CSparse.Storage.CoordinateStorage<double>;
var buf = new Coord(10, 10, 1);
buf.At(2,3,1.0);
I tried to use the Francis algorithm to obtain eigenvalues of a sparse matrix, but it always diverged to NaNs. So, I started by checking if the QR factorization works correctly. I used the following code:
The results seem to be wrong, because the bottom-right element in R seems to be negated compared to results in Octave, and, most importantly, multiplying Q * R does not give A:
A:
7 2
2 1
Q:
14.2801098892805 0
2 0.824163383692134
R:
-7.28010988928052 -2.19776902317902
0 -0.412081691846067
Q * R:
-103.960769224964 -31.3843831622532
-14.560219778561 -4.73516068786748
Octave code and results for comparison:
A = [ 7 2; 2 1 ]
A =
7 2
2 1
[Q, R] = qr(A)
Q =
-0.96152 -0.27472
-0.27472 0.96152
R =
-7.28011 -2.19777
0.00000 0.41208
Q * R
ans =
7.0000 2.0000
2.0000 1.0000
Hello,
First, thanks a lot for this library, I'm glad that there is finally a good native C# open-source sparse matrix library! I have a beginner's question regarding the solution of underdetermined systems.
I am currently doint a full C# implementation of the Heat Method by Keenan Crane et al., which requires to solve a linear system Ax=y
, where A
is a square matrix, which is not full-rank (A
is actually the Laplacian matrix of a graph, so the vector (1,1,1,1,1,1,1...,1)
is in the null space of A
). I would like to solve the system in a least square sense, and I have read that QR factorization is well-suited for this purpose.
I applied the same code as the sparseQR example, but with a laplacian matrix. The result of Solve is a vector proportional (1,1,1,1,1,1,1...,1)
with a very large norm (each component of the vector has a value superior to 1 million, while the y
vector has components inferior to 10.
Is there a specific trick to solve this in the least square sense with the SparseQR Factorization?
Best,
Romain
Hi
The lowest .NET version that is supported by the csparse nuget package is net 4.6.1 by now. Is it possible to also include older ones, like net45 or even fossile ones like net40?
Thanks
Currently method for get or set member at i,j
of a Double.DenseMatrix
is the SparseMatrix.At
method,
Also possible to have indexing access for the instance:
/// <summary>
/// Gets or sets the member at i,j
/// </summary>
/// <param name="i">the row</param>
/// <param name="j">the column</param>
/// <returns>member at location i,j</returns>
public double this[int i, int j]
{
get
{
return this.At(i, j);
}
set
{
this.At(i, j, value);
}
}
Final result is member at i,j
can be accessed in C# like this:
metrix[i,j] = value
Hi,
How are you?
I think it would be good if library does not throw exception if someone wants to create a 0x0 matrix. Logically it does not makes any sense to create a 0x0 matrix, but say in fem problems it is possible in situations that a 0x0 matrix should be created.
Thanks
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.