spFindElement

[Stuart Thorncraft]

Hi,

I have been playing around with Sparse 1.4. A couple of questions / comments:

* There doesn't appear to be any support for basic matrix operations like multiply and add - is this right?

* Is it possible to create non-square sparse matrices?

* I wrote some code (see below) to iterate over a sparse matrix to investigate the elements contained in the matrix.  I noticed that the spFindElement(...) function seems to return the element residing in a given memory location rather than the position within the matrix.  This is inconsistent with how the function spGetElement(...) works.  Is there another way I can get the elements in the matrix without having to create a new entry in the sparse matrix? Or is my use of spFindElement(...) completely wrong.  The output of the program is summarised below.

Thanks very much for your assistance

Stuart

-------------OUTPUT--------------

Size of matrix = 3 x 3.

Matrix before factorization:
             3         1         2

   3         1       ...       ...
            2j
   1       ...       ...       0.3
                                1j
   2       ...         4       ...
                    0.5j

Largest element in matrix = 4.5.
Smallest element in matrix = 1.3.

Largest pivot element = 3.
Smallest pivot element = 3.

Density = 33.33%.

( 1, 1) = 1.000  |  2.000j
( 2, 3) = 0.300  |  1.000j
( 3, 2) = 4.000  |  0.500j

-------------------------

----------CODE---------

#include <stdio.h>
#include "spMatrix.h"

void main()
{
    spMatrix pM;
    spElement * pE;
    spError err;
    pM = spCreate(3,1,&err);

    pE = spGetElement(pM,3,3);
    spADD_COMPLEX_ELEMENT(pE,1.0,2.0);

    pE = spGetElement(pM,1,2);
    spADD_COMPLEX_ELEMENT(pE,0.3,1.0);

    pE = spGetElement(pM,2,1);
    spADD_COMPLEX_ELEMENT(pE,4.0,0.5);

    spPrint(pM,1,1,1);

    for (int i=1; i<=3; ++i)
    {
        for (int j=1; j<=3; ++j)
        {
            pE = spFindElement(pM,i,j);
            if (pE != 0)
                printf("(%2d,%2d) = %4.3f  |  %4.3fj\n", i, j, *pE, *(pE+1));
        }
    }
}

[Ken Kundert]

Stuart,
1. You are correct, there is no support for basic matrix operations like multiply and add. Just matrix factorization and back-solve, matrix-vector multiplication, and limited support for transposing.

2. There is no support for non-square matrices.

3. I do not understand your comment about the inconsistency between spGetElement() and spFindElement(). They seem quite consistent to me. They take the same arguments and provide the same return value. The only difference is that one creates the element if it does not exist, and the other does not.

The function spFindElement() was added on to the program as a result of user request. I personally have never found a use for it. Be aware that it is not coded very efficiently. Each call will result in either a row or column search.

-Ken

[Stuart Thorncraft]

Hi, thanks for the responses. 

In relation to the third point:

The code I provided puts the following into a sparse matrix:

1+2j at position (3,3)
0.3+1.0j at position (1,2) &
4+0.5j at position (2,1)

Then when I execute the spFindElement command on element (1,1) I get a non-null spElement back, indicating that something must be in position (1,1) but there isn't.  If instead, I was to execute the spGetElement command for position (1,1), I would end up putting a 0+0j into the sparse matrix at position (1,1). 

The results returned by spFindElement seem to depend on the order in which you place elements into a Sparse matrix, while spGetElement doesn't.  Is this intended?

[Ken Kundert]

Stuart,
    I think you have found a bug. As I mentioned last time, spFindElement() has only recently been made available for external use; it was originally an internal utility function. In its roll as an internal utility, it used only internal row and column numbers. When it was opened up, I don't believe I added the code to properly translate the external row and column numbers to internal row and column numbers before I performed the look up.

So now my choice is to either remove the function or fix it. Do you need it?

-Ken

[Stuart Thorncraft]

Ken,
Presently I am looking at writing some helper functions to add/multiply the sparse matrices. Given that I am looking at using the code in a situation where performance is not super-critical, I was using the spFindElement in a fairly dumb algorithm.  So fixing its implementation is not really critical.

What is the best way to iterate over the complex elements of a sparse matrix?  I guess I should take a look at the code behind spPrint(...)

Thanks
Stuart

[Ken Kundert]

spPrint() is probably not the best function to use as an example. Try spMultiply() and spMultiplyTransposed(). One shows how step through the matrix by rows, the other by columns.

-Ken

[Smeuuh]

Hi,
Not sure anyone will actually read this message, but well :)
Thank you for spFindElement, I actually wrote it myself in sparse 1.3, for debugging purposes.
spFindElement doesn't translate indices, which means its behavior is inconsistent with spGetElement.
The fix is quite simple, just add the beginning of spGetElement to the beginning of spFindElement. It's quick, and it works fine here.
Thank you a lot for the library, it gives us a big speedup over our naïve dense solver in our simulator (http://edacious.hypertriton.com)
Antoine Levitt