wo80 / csparse.net Goto Github PK

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);

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

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.

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

To start: Awesome project. Would AVX vectorization with System.Numerics.Vector be possible? Is the original CSparse making use of SIMD?

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,
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);

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

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;
}

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?

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

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

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

I can get the value using matrix.At(row, col), but I can't find a way to set the value.

Hello,
For making the CSparse.NET CLSCompliant, these minor changes are needed:

• Convert Matrix.rowCount and Matrix.colCount accesors from protected into internal
• Change name of either 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

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, 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

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.

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); 

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 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,
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

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

Hi,
Is there code for LDLt decomposition of symmetric matrices with the library?

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

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

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!

Feature request

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!

Issue

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):

• Transposing a CSparseMatrix yields incorrect results.
  E.g. when transposing sparse matrix a

  var a = (SparseMatrix)SparseMatrix.Build.SparseOfRowArrays([[1, 2], [3, 4], [5, 6]]);

  using

  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();

  one would expected aDash to be

  1 3 5
  2 4 6
  

  in matrix notation.
  However, I instead get

  1 2 0
  3 4 0
  

  or

  %%MatrixMarket matrix coordinate real general
  2 3 4
  1 1 1
  2 1 3
  1 2 2
  2 2 4
  

  in Matrix Market format (obtain via 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
  

  although its Values property is [1,2,3,4,5,6] (with 3 rows and 2 columns).
• Multiplying matrices is unclear to me.
  Is supposed A = B*C equivalent to var A = B.ParallelMultiply(C) or var A = C.ParallelMultiply(B)?
  When I execute var aTilde = aDash.ParallelMultiply(aSparse); I get

  5 11
  11 25
  

  and with var aTilde = aSparse.ParallelMultiply(aDash); I get

  10 14 0
  14 20 0
  0 0 0
  

  .
  Mathematically, aDash*aSparse and aSparse*aDash should be

  10 14
  14 20
  

  and

  5 11 0
  11 25 0
  0 0 0
  

  respectively.
  Here, I used the values for aDash and aSparse that were computed by the library (i.e. the unexpected ones).
• These issues of course propagate to or also happen in the Solvers.
  E.g. var solver = SparseQR.Create(aTilde, ColumnOrdering.Natural);
  yields incorrect factorization matrices and applying Solve (with ColumnOrdering.Natural) yields the wrong result.

My original mathematical problem

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.

Versions

• .NET SDK: 8.0.202
• CSparse 3.8.1
• MathNet.Numerics 5.0.0
• MathNet.Numerics.MKL.Win-x64 3.0.0

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 yvector 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 @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

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

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:

Program.zip

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,

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

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,
You have any existing code that be able to use Intel's MKL or CUDA for cholesky factorization of a CompressedColumnStorage ?