Mar 9, 2011 at 4:46 PM
Edited Mar 9, 2011 at 4:57 PM

I am working on a group of classes that handle Symmetric dense matrices. I need some input on structure/architecture/design/whatever_this_is_called to make it fit nicely.
Symmetric matrices provide an immediate memory advantage. A (square) dense matrices has n^2 values whereas a symmetric dense needs n*(n+1)/2 (roughly half+) Also, we save a lot on calculations (eg when creating it we only need to set half+ the values etc).
It would be very advantageous to be able to specify when a matrix is symmetric. By convention, I am saving the upper triangle of the matrix rowwise (the partial row from the diagonal element and to the right, then the next partial row etc).
A symmetric matrix can be saved with 2 major ways:
1) Vector form. Since we know exactly how many elements are saved on each row, it is easy to calculate where an element [i,j] is located on the vector and access it.
2) Jagged Array form(eg double[][]). Each subarray contains exactly as many elements as needed. Eg the first one contains n elements, the second one contains n1 elements etc.
I am strongly in favor of the second form because its a lot clearer and maintainable. However, I noticed that even for DenseMatrices MathNet uses double[]. I would think that having to perform (column * RowCount) + row to access every element would be less
efficient than double[,] (even though [,] would make that sort of calculation internally) but I am sure this has been researched (I think Iridium used to use [,] or [][]). Please clarify which form is preferred.
Next issue: For the DenseMatrix class of Double/Single/Complex/Complex32 we got:
DenseMatrix : Matrix : Matrix<T>
Symmetric is a property of the matrix, so an interface ISymmetric seems natural. But since I cannot provide a body to the methods, I thought of creating a new abstract class next to Matrix. Currently I am thinking of:
SymmetricDenseMatrix : SymmetricMatrix : Matrix<T>
for each one of Double/Single/Complex/Complex32. This is would help in creating a: SymmetricSparseMatrix : SymmetricMatrix : Matrix<T>
How does this seem? I also need to take into account that the matrix is always square. Do I need a SquareMatrix<T>?
All input is greatly appreciated :).


Mar 9, 2011 at 6:06 PM
Edited Mar 9, 2011 at 6:07 PM

A symmetric matrix would be great addition.
>However, I noticed that even for DenseMatrices MathNet uses double[]. I would think that having to perform (column * RowCount) + row
>to access every element would be less efficient than double[,]
We are trying to keep the matrix storage compatible with native BLAS/LAPACK matrix storage schemes. For dense matrices, the matrices are stored as a single dimensional array with LAPACK requiring the data to be stored in a columnmajor format (rowmajor
might be an option now with the new C API for LAPACK, but I haven't checked and I don't know how well that is supported with different LAPACK implementations). I believe [,] is stored in contiguous memory in .NET so it can passed as a single
dimensional array to the native libraries, but it will be in rowmajor form. Also in earlier versions of the .NET runtime, [,] based matrices were exceptionally slow. I don't know if the newer versions of the runtime have improved performance. It might be
time to redo our benchmark comparing [],[][], and [,] ([] were the fastest and [,] were the slowest in the benchmarks we did years ago for dnAnalytics even with the [] indexing).
As for the dense symmetric matrix storage, I think we need to use packed symmetric storage in a 1D array to be compatible with the native libraries. I'll double check tomorrow, as well see if the CSB sparse storage scheme is compatible.
> Do I need a SquareMatrix<T>?
I wouldn't add one.



If we want to take advantage of native BLAS/LAPACK routines for the symmetric matrix, we should use packed storage. See this link for a description: http://software.intel.com/sites/products/documentation/hpc/mkl/webhelp/appendices/mkl_appB_MA.html
We are planning on using native routines for the sparse matrix, so we need to select a storage scheme that is compatible with MKL (and hopefully other providers). Here is a list of the MKL sparse storage formats: http://software.intel.com/sites/products/documentation/hpc/mkl/webhelp/appendices/mkl_appA_SMSF.html
We'd have to add a second sparse matrix type if we wanted to support CSB storage. User's could then select either based on their needs.



Ok, I will change the symmetric storage to the scheme used in MKL.



Are there methods in the LinearAlgebraProvider(s) for symmetric matrices?



No. We are focusing on dense matrices for the first release. Fell free to modify the LinearAlgebraProvider interface to handle them.

