What is the largest allowable size of matrices and vectors to pass to mldivide?

I am using mldivide to solve a sparse matrix problem within a finite element code. The matrices are defined as sparse from the begining by

A = sparse(nGdof,nGdof)

For a system of size around 16000x16000, the code works well. For larger matrices I get strange results. What brings me to the question, is there a limit of size on mldivide? or more generally, for sparse matrices in general?

[Sean de Wolski]

What do you mean by "strange" results? Please clarify; are your results: wrong values, wrong sizes, not being attained etc?

[Andrew Newell]

It would also help to know how many nonzero elements there are.

[gabriel villalobos]

Thanks for the interest!

This is a finite element code, so I solve a matrix equation of the sort:

U=stiffness(activeDof,activeDof)\force(activeDof);

where activeDof is a vector of positions in the matrix, stiffness is sparse. U is used to produce a new force vector, that I display in the plot. The strange result is that this new force vector is zero for most of the entries.

Which is why I asked, ?is there a limit size on mldivide?

thanks!

Answer 1

[Walter Roberson]

The size limit is only due to available memory.

What is the condition number of the system? If your system is not well conditioned, you could be hitting underflows that propagate as zeros.

[Sean de Wolski]

Something very plausible if the boundary conditions aren't applied or aren't applied properly to make the system non-singular...

[gabriel villalobos]

When the matrix is non invertible mldivide returns an error that stops the code, so this is not the case either!

[gabriel villalobos]

Thanks both for answering. But in that case matlab would stop... this is not the case here.

[Sean de Wolski]

Not necessarily true, it can be invertible at the precision of which you're working...

[Walter Roberson]

It is not clear to me that "no underflows" would be the same as "non invertible".

[gabriel villalobos]

Thanks Sean and Walter.

I am applying this many times, the matrix changes each time, and thus the condition changes also. But, for instance now the last condition number I printed was:

7.084516338167707e+27

[Walter Roberson]

In http://www.mathworks.com/matlabcentral/newsreader/view_thread/18017 Jeffery J. Leader wrote,

>> As a rule of thumb, if the condition number is roughly 10^k, expect to lose about k significant digits when solving Ax=b.

If that estimate holds for your circumstances, then you are losing more digits (about 27) than exist in your data (about 16).

The estimate Jeffrey Leader provided agrees with something John d'Errico posted several years later, indicating that a "too large" condition number was roughly 1/eps as that would be about 5E15

[gabriel villalobos]

Sean and Walter, thanks a lot.

Answer 2

[Sean de Wolski]

No, there is no limit. Does your stiffness matrix have boundary conditions applied?

EDIT example for clarification

K = [1 2 3;4 5 6;7 8 9];
Kbc = [1 0 0;0 5 0;0 0 1]; %DOFs 1,3 fixed

[gabriel villalobos]

The boundary conditions are given in the form of the activeDof vector. This is, activeDof only has the number of the DOF of the system which are free to move.

[Sean de Wolski]

And thus your system is unstable!!!

If everything can move and you apply a force, will it ever stop? Infinite displacements!! (well in theory)

[Sean de Wolski]

Instead, apply the BCs to your stiffness matrix before hand by zeroing out all rows and columns of it that correspond to fixed DOF and setting a 1 to the diagonal of those DOFs.

[gabriel villalobos]

No, the system is not unstable. For M*x=u, with M a matrix, x & u vectors, there is solution if M is invertible.

[Sean de Wolski]

"gabriel villalobos about 3 hours ago

The boundary conditions are given in the form of the activeDof vector. This is, activeDof only has the number of the DOF of the system which are free to move"

Thus extracting only that parts that are free to move leave you system not stable.