[Feature request] Add Matrix Multiplication about csparse.net HOT 8 CLOSED

1. Actually, you shouldn't use operator overloads at all. Using MathNet, you are tempted to write code like -KgGlobalFree * EigenVectorEstimate, but this hides some important internals from the user. For example, the unary minus operator for sparse matrices will create a copy of the data. Doing this inside a loop is a bad idea, if you care about performance (lots of unnecessary memory allocations). The same argument holds for all operator overloads (matrix-vector etc.).

2. Since KeGlobalFree isn't modified inside the loop, you should create the factorization outside the loop.

3. Why don't you use CSparse all the way? The following C# code wasn't tested. It's just to get an idea:

using CSparse.Double;
using CSparse.Double.Factorization;
using CSparse.Storage;

void Example()
{
    int N = FreeDOFs.Count;

    // Use coordinate storage to assemble the matrices.
    var cKeGlobalFree = new CoordinateStorage<double>(N, N, 2 * N);
    var cKgGlobalFree = new CoordinateStorage<double>(N, N, 2 * N);
    
    // Convert to sparse matrix.
    var KeGlobalFree = CSparse.Converter.ToCompressedColumnStorage(cKeGlobalFree);
    var KgGlobalFree = CSparse.Converter.ToCompressedColumnStorage(cKgGlobalFree);

    double Errors = 10;
    double EigenValueEstimate = 1;

    var EigenVectorEstimate = Vector.Create(N, 1);

    var y1 = Vector.Create(N, 0);
    var y2 = Vector.Create(N, 0);
    
    var SC = SparseCholesky.Create(KeGlobalFree, CSparse.ColumnOrdering.MinimumDegreeAtPlusA);

    while (Errors >= 10e-8)
    {
        // y1 = -KgGlobalFree * EigenVectorEstimate
        KgGlobalFree.Multiply(-1, EigenVectorEstimate, 0, y1);

        // EigenVectorEstimate = KeGlobalFree \ y1
        SC.Solve(y1, EigenVectorEstimate);

        // y2 = -KgGlobalFree * EigenVectorEstimate
        KgGlobalFree.Multiply(-1, EigenVectorEstimate, 0, y2);

        // norm = EigenVectorEstimate * -KgGlobalFree * EigenVectorEstimate
        double norm = Vector.DotProduct(EigenVectorEstimate, y2);

        // EigenVectorEstimate = EigenVectorEstimate / sqrt(norm)
        Vector.Scale(Math.Sqrt(norm), EigenVectorEstimate);

        // y2 = KeGlobalFree * EigenVectorEstimate
        KeGlobalFree.Multiply(1, EigenVectorEstimate, 0, y2);

        // NewEigenValueEstimate = EigenVectorEstimate * KeGlobalFree * EigenVectorEstimate
        double NewEigenValueEstimate = Vector.DotProduct(y2, EigenVectorEstimate);

        Errors = Math.Abs((NewEigenValueEstimate - EigenValueEstimate) / EigenValueEstimate);

        EigenValueEstimate = NewEigenValueEstimate;
    }
}

from csparse.net.

Are you multiplying a sparse matrix with a sparse vector, or sparse matrix with a dense vector?
I think dense vector multiplication is already in the library.

from csparse.net.

I am multiplying sparse matrices with dense vectors.
How do i use it? Do i just use the multiplication operator like A * B? where A is double(,) and B is double()

Thanks.

from csparse.net.

I didn't know that a minus sign creates a copy of the data.
Thanks for your help! I will implement your suggestions.

from csparse.net.

wo80, I would like to thank you again! Your algorithm is a lot faster, about 70 times faster for a 700x700 sparse matrix! I still haven't utilized your idea of using CoordinateStorage because when i assign cKeGlobalFree it is always in an additive manner (cKeGlobalFree (i,j)=cKeGlobalFree (i,j)+x) and i couldn't figure out how to do this with CoordinateStorage.

For the benefit of someone who might be looking at this later, the following line of code should be changed:
From: Vector.Scale(Math.Sqrt(norm), EigenVectorEstimate);
To: Vector.Scale(1.0/Math.Sqrt(norm), EigenVectorEstimate);

I realize i am going off topic here, but do you know of any fast C# EigenProblem Solvers?
The EigenValue solver in Math.Net is really slow as it finds all eigenvalues. In reality, someone might be interested in the first couple EigenValues. The code above uses the inverse iteration method and is for solving the buckling problem. There are faster, more complicated methods out there such as the subspace iteration method. Any suggestions would be appreciated.

from csparse.net.

The CoordinateStorage will automatically add up all values you pass in, so

// Create a 2x2 matrix.
var C = new CoordinateStorage<double>(2, 2, 4);

// Add some values.
C.At(0, 0, 1.0);
C.At(1, 1, 1.0);
C.At(1, 1, 0.5); // Again position (1 ,1)

// Convert to sparse matrix.
var A = CSparse.Converter.ToCompressedColumnStorage(C);

will produce the matrix A with entries

[ 1.0  0.0 ]
[ 0.0  1.5 ]

Unfortunately, I don't know any sparse Eigenproblem solvers for .NET. I've used ARPACK / ARPACK++, calling the native code from C# using P/Invoke.

from csparse.net.

Thanks!
As you might have already noticed, i am not a programmer by training. Do you have any example using ARPACK in C# ?

I looked at http://sweb.cityu.edu.hk/kincau/blog/files/0e9a2b646554191e9ddc831b697cf04d-1.html and tried it in visual studio but i am left with 3 errors!

from csparse.net.

Even if the C# code would compile, you'd still need to build the native ARPACK DLL, which is rather involved.

I've used a lot of native solvers from C# and I'm planning to release all the code here on Github (including ARPACK), but at the moment, I don't have much time to work on it.

from csparse.net.