Sparse sub access

This package allows to retrieve and assign values of sparse matrix in one shot.

Description

The indexes must be provided in (row,column) form. Functions in this package are programmed with care so no overflow of linear index is used with large-size matrix. I try to design an engine of automatic selection algorithms so as to provide a good run time performance across various cases. I'm not sure it's worth the effort, but my experience shows that it is easy to mess with the performance when non-appropriate method is using in sparse manipulation. This package does not warranty for best performance, but it just has a modest pretension to prevent user to use bad algorithm. For a better customize solution, a calibration is required on the specific
computer by running this command: spcalib('auto')

• Acknowledgements: This submission has inspired the following:; Robust Sparse data types
• MATLAB release: MATLAB 7.7 (R2008b)

[Anthony Lamb]

Very good!!!! It works well for building global stiffness matrices for finite element computations. While building my stiffness matrices I currently do the following
      val = getspvalmex(A_global, Ic, Jr);
      A_global = setspvalmex(A_global, Ic, Jr, val+X);

Bruno...Is there any other way your code can be used to get this done??

[Anthony Lamb]

Sorry about that.... I think that the following code does the trick
 
A_global = setspvalmex(A_global, Ic, Jr, X,@plus);

[Matt J]

Bruno - When I do

 mex -v -O setspvalmex.c

I get a compilation error:

Error setspvalmex.c: 626 compiler error in _kids--Bad rule number 0
 
  C:\PROGRA~1\MATLAB\R2009B\BIN\MEX.PL: Error: Compile of 'setspvalmex.c' failed.

Any idea what I'm doing wrong?

[Bruno Luong]

Hello Matt,

After investigation, it turns out the "_kids--Bad rule" message is caused by a bug of the LCC compiler.

A workaround is to replace the or "|" operator by plus "+" operator in the cmex file setspvalmex.c, lines #626, 634 and 709, then recompile.

Thanks for report the bug. I'll submit a new version.

Bruno

[Matt J]

Yep, that workaround fixed it. Anyway, I like it!!

[James Tursa]

Bruno,
FYI, my experience is that lcc has bugs with the unsigned long long type ... I often get a "bad rule" compile error. I have gotten into the habit of using signed long long instead and writing extra code to handle situations where I need to shift into or out of the sign bit. A real hassle just so the lcc compiler that ships with MATLAB will work with my code.

[Bruno Luong]

Thanks for the tip James.

[H]

Bruno, Thank you very much for this package. Also thank you for all your help in the forums.

This works appr. 30 times faster than using linear indexing to zero sparse matrix columns that I used before!

[H]

This seems to have one problem though:

Before I used A(:,indexes)=0, it just took a lot of time.

No I tried

[I J]=find(A);
idxofzero = ismember(J,indexes);
A=setsparse(A,I(idxofzero),J(idxofzero),0);

And I get "out of memory error". I tried it several times with both methods and got the same result. Am I missing something here?

My sparse matrix is about 55,000 x 55,000 and it has 4,060,200 nnz elements. The indexes hold about 1/100 of the columns.

[Bruno Luong]

Hello H,

I can't really help with memory problem. The suggestion I can give is to try on your side to increase the free largest memory block from other places. Thanks for the kind words.

Bruno

[shuangou ren]

 

[shuangou ren]

good job, two times faster for finite element assemblage.

[LeFlaux]

Hey Bruno!

In my FE-Code the global matrices are assembled as usual in a loop over the elements, eg:

A(idx, idx) = A(idx, idx) + locA;

where idx defines a sub-matrix of A (eg idx = [3,1,4,12,7,10] and locA is 6x6). I translated this code to Your functions as follows using row and columns indexing:

[cidxs, ridxs] = meshgrid(idxs, idxs);
A = setsparse(A, ridxs(:), cidxs(:), locA(:), @plus);

Unfortunately, this is (little) more expensive than the original version. Do You have any suggestions how to shorten things above a little?

Thank you in anticipation for your support.