This project has moved and is read-only. For the latest updates, please go here.

Iterative solver for Complex matrix

Dec 20, 2010 at 10:30 AM

Hello,

I tried to use SparseMatrix of Complex type with an iterative solver and get an OutOfMemory exception. I haven't looked at your code, but I guess the iterative solvers shouldn't call any ToArray methods.

My matrix has a size of approx. 31500x31500, so making it an array is probably no good idea :-)

Here's the exception

   at MathNet.Numerics.LinearAlgebra.Generic.Matrix`1.ToArray()
   at MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners.IncompleteLU.Initialize(Matrix`1 matrix)
   at MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative.BiCgStab.Solve(Matrix`1 matrix, Vector`1 input, Vector`1 result)
   at MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative.BiCgStab.Solve(Matrix`1 matrix, Vector`1 vector)
   

Since I use sparse matrices a lot:

Is there a way to iterate just the non-zero entries (getting column, row and value in each iteration)?

Is there a way, to check the symmetry of the sparse matrix, without having to iterate all rows and columns?

 

Thanks,

Chris

Dec 20, 2010 at 2:36 PM
This discussion has been copied to a work item. Click here to go to the work item and continue the discussion.
Dec 20, 2010 at 2:41 PM
This discussion has been copied to a work item. Click here to go to the work item and continue the discussion.
Dec 20, 2010 at 2:53 PM
This discussion has been copied to a work item. Click here to go to the work item and continue the discussion.
Dec 20, 2010 at 2:55 PM

Hi Chris, 

>My matrix has a size of approx. 31500x31500, so making it an array is probably no good idea :-)
That was implemented poorly. I'll rework it later this week.

>Is there a way to iterate just the non-zero entries (getting column, row and value in each iteration)?
>Is there a way, to check the symmetry of the sparse matrix, without having to iterate all rows and columns?
I've created work items for these. I'll get to them in a couple weeks.

For the first, you could use the RowEnumerator that returns a row, vector KeyValuePair (should that be changed to a tuple?). Then use the IndexedEnumerator on the vector - returns a index, value pair. So indirectly you'd be able to iterate over all the non-zero elements. The IndexedEnumerator would work better.

Regards,
Marcus

Dec 20, 2010 at 6:58 PM
Edited Dec 20, 2010 at 6:59 PM

Thanks again.

The following code will plot the sparsity pattern to a bitmap:

 

private Bitmap Spy(SparseMatrix m)
{
    // Size of the bitmap
    int size = 500;

    // Create bitmap
    Bitmap bmp = new Bitmap(size, size);

    // Fill background
    Graphics g = Graphics.FromImage(bmp);
    g.Clear(Color.White);
    g.Dispose();

    // Size of the matrix
    int N = m.ColumnCount;

    var rows = m.RowEnumerator();

    foreach (var row in rows)
    {
        var nonzeros = row.Value.GetIndexedEnumerator();

        foreach (var col in nonzeros)
        {
            bmp.SetPixel((col.Key * size) / N, (row.Key * size) / N, Color.Black);
        }
    }

    return bmp;
}

 

Unfortunately, the performance isn't quite what I expected. My 31500x31500 matrix has ~220000 nonzeros. The loops take about 150 seconds on my pc. That's pretty much!

Is there a way to speed things up?

Regards,
Chris

Dec 21, 2010 at 7:22 AM
This discussion has been copied to a work item. Click here to go to the work item and continue the discussion.
Dec 21, 2010 at 7:26 AM

The problem seems to be that we haven't optimized the Row/ColumnEnumerators for sparse matrices. I've added a work item. I'll implement the IndexedEnumerator first since it will perform better than an optimized version of RowEnumerator.

Regards,
Marcus 

Dec 25, 2010 at 4:37 PM
wo80 wrote:

Hello,

I tried to use SparseMatrix of Complex type with an iterative solver and get an OutOfMemory exception. I haven't looked at your code, but I guess the iterative solvers shouldn't call any ToArray methods.

My matrix has a size of approx. 31500x31500, so making it an array is probably no good idea :-)

Here's the exception

 

   at MathNet.Numerics.LinearAlgebra.Generic.Matrix`1.ToArray()
   at MathNet.Numerics.LinearAlgebra.Complex.Solvers.Preconditioners.IncompleteLU.Initialize(Matrix`1 matrix)
   at MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative.BiCgStab.Solve(Matrix`1 matrix, Vector`1 input, Vector`1 result)
   at MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative.BiCgStab.Solve(Matrix`1 matrix, Vector`1 vector)

Fixed in commit 029c2bbc44e4.