Solving [K]{x}={b} efficiently, when [K] is symmetric!

Hi,

I have a fully populated symmetric matrix [K] (like those found in the finite element method).. but to save space for big [K], only the upper triangular is stored (packed) DIRECTLY either in a rectangular matrix OR in a long single column vector, (nearly a factor of 2.0 savings in space this way).

The question: Is there a matlab m file for the upper triangular back substitution matrix used in a Gauss elimination LDU solve for solving linear equations like [K]{x}={b}

Thank you,

Answer 1

[PA00]

In general FEM matrix’s are symmetric but not fully populated, also if you write your matrix with attention to the FEM mesh you can most of the times make the [K] matrix to be a banded matrix.

Since your matrix is fully populated it is not advantageous to use sparse matrix methods and matlab does (like Bruno said) not have symmetric matrix storage capabilities.

You can always write a function to solve your particular matrix, there are some efficient FORTRAN subroutines in manuals that do that for very particular types of matrices and you can rewrite them to matlab however I suggest you write your own. I have written my own to solve FEM symmetric and banded matrices where I store only the upper diagonals into a rectangular matrix and then solve it using Gaussian elimination code with a small twist to only reach that banded part. This was done to use less memory storing the matrix and use less processing time to solve the system.

Anyways, for your specific question you can use “linsolve(A,B,opts)” where in opts you state the type of A matrix you are inserting like a positive definite and upper triangular matrix (see the help file to see how it works).

This function works nice however has a bit of a downside, it does not accept matrix’s stored in a sparse system. It would be great if matlab made a specific matrix store system for symmetric matrix (and banded also if it was not much to ask) for memory savings and also specific solving code for faster solving.

Answer 2

[Bruno Luong]

Short answer no.

Long answer:

MATLAB does not have a mechanism of storing symmetric matrix (saving by a factor of 2 is probably not worth the effort).

Long answer: the CHOL command looks only at the upper half of your matrix, but as stated above you have to provide the full matrix, so I can't see how it helps you, without even look at the requirement of definite positiveness of the matrix.

BTW, FEM matrix is usually sparse. Now that is something worth the effort.

But you should be aware that coding a half of fully-populated matrix in sparse will take even more memory.

To summarize: just don't bother if you want to save half of fully populated matrix.

Bruno

[Liang]

This Q/A was dated 2011. Any added comment for 2024? Thanks.

[Walter Roberson]

Not much has changed since 2011.

You could potentially call directly into a BLAS routine; https://www.mathworks.com/help/matlab/matlab_external/calling-lapack-and-blas-functions-from-mex-files.html

[Bruno Luong]

And I predict it won't change in the next 20 years.

[Steven Lord]

I don't know if it would save time, but the iterative linear equation solvers do accept function handles that accept a vector and compute certain functions involving the coefficient matrix and that vector. If you can compute those certain functions without explicitly creating the coefficient matrix, this might allow you to solve larger problems that wouldn't fit in memory. For example, for the gmres function:

You can optionally specify the coefficient matrix as a function handle instead of a matrix. The function handle returns matrix-vector products instead of forming the entire coefficient matrix, making the calculation more efficient.

To use a function handle, use the function signature function y = afun(x). Parameterizing Functions explains how to provide additional parameters to the function afun, if necessary. The function call afun(x) must return the value of A*x.

If you know how to compute A*x without explicitly forming A, as in the "Using Function Handle Instead of Numeric Matrix" example on the gmres documentation page, you might want to try one of these solvers.

[Bruno Luong]

@Steven Lord is there any matrix algebra function for tall arrays ? I don't think but you know MATLAB features inside out could prove me wrong.

Few other PDE simulation sofware store the matrix on disc to leverage the dimension of the problem that can be solved.

Also missing in MATLAB is some sort of "sparse tall array", some sort of sparse matrix that can be loaded by section into RAM.

II have no idea how many users outthere are interested on such extension, probably very small though.

[Liang]

Thanks again. Please comment.

Among the many different aspects of FEM solution, the two major solution issues are in the original 2011 question, i.e., banded/sparse and symmetry of K.

After reading your recent responses, I looked at the matlab "\" or mdivide documentation. From the diagram, it appears to me that, in order for matlab to detect if a sparse K is symmetric, a complete K must be formed before "\" is used. So, your answer 13 years ago does not need any update.

What can be done differently, potentially, is to use combination of amd (or other similar) and ldl (amd and ldl accept upper/lower sparse matrix but I don't know if the codes actually form the complete A+A') to compute L and then exercise "\". ldl should not need to generate A+A' to compute L. I know very little about amd.

Iteration solution should be a very effective way to take advantage of the band/sparse and symmetry of K. However I have been using "\" for linear equation solution. Don;t have any user experience on iteration solution.