Thread Subject: FEM substructuring (static condensation) -> out ot memory; ANSYS can do it

Hi,

I'm coding my own FEM code in Matlab and want to implement
the so-called substructuring technique, which basically is a
static condensation.

However Matlab runs out of memory quiet fast. As soon as I
want to reduce my stiffness matrix by about 5000dof the
calculations fails.

Here's the formula for substructuring, so you can get an
idea of it, if you don't know it.

K*u = F

u= [ua; ub];
F = [Fa; Fb];
K = [Kaa Kab; Kba Kbb]
Kab=Kba';

ua are the degrees of freedom i want to keep, whereas ub are
the ones which should be removed. Therefore I obtain a new
matrix for the left side and a new vector on the right side;

K_new = [Kaa - Kab*Kbb^-1*Kba]
F_new = [Kaa - Kab*Kbb^-1*Fb]

Now of course I'm not computing the full inverse of Kbb, but

Kaa-Kab/Kbb*Kba

I'm not quiet sure why Matlab runs out of memory when Kbb
reaches about n=5000 (as I mentioned above). Either it's the
computation of Kab/Kbb or maybe as Kab/Kbb most probably
becomes a full matrix it just needs to much memory.

The only thing I know for sure is that ANSYS can calculate
that kind of substructuring without any problems. It even
manages to remove well over 100 000 degrees of freedom
without any major computing needs (takes only about 5 min).

So there must be a major difference in the way Matlab (or
myself) is calculating K_new compared to ANSYS.

Is it the way matrix are being stored. As far as I know
ANSYS uses Harwell-Boeing (not quiet sure). This maybe is
more memory efficient then the way Matlab is storing sparse
matrixes, but that still wouldn't explain why the
calculation fails.

Any ideas or reason on how I would have to change my code?

Thanks a lot!!!

"A B" <gitsnedbutzi@hotmail.com> wrote in message

...
> K*u = F
>
> u= [ua; ub];
> F = [Fa; Fb];
> K = [Kaa Kab; Kba Kbb]
> Kab=Kba';
>
> ua are the degrees of freedom i want to keep, whereas ub are
> the ones which should be removed. Therefore I obtain a new
> matrix for the left side and a new vector on the right side;
>
> K_new = [Kaa - Kab*Kbb^-1*Kba]
> F_new = [Kaa - Kab*Kbb^-1*Fb]
>
> Now of course I'm not computing the full inverse of Kbb, but
>
> Kaa-Kab/Kbb*Kba

what are the dimensions of these matrices? Kbb is dimension
5000?

What's happening here is that you're computing your own
Schur complement. It might be too much to represent this
explicitly. Perhaps ANSYS does this implicitly, and
updates/downdates the factorization of K accordingly.

MATLAB uses the same compressed-sparse column form that the
Harwell-Boeing format uses (mostly, except for different
methods used inside LU, Chol, backslash, etc, that aren't
visible to the user - M or mex or otherwise).

> what are the dimensions of these matrices? Kbb is dimension
> 5000?
>
> What's happening here is that you're computing your own
> Schur complement. It might be too much to represent this
> explicitly. Perhaps ANSYS does this implicitly, and
> updates/downdates the factorization of K accordingly.
>
> MATLAB uses the same compressed-sparse column form that the
> Harwell-Boeing format uses (mostly, except for different
> methods used inside LU, Chol, backslash, etc, that aren't
> visible to the user - M or mex or otherwise).
>

Hi Tim,

Kbb could be larger than 5000, but 5000 is approximately the
dimension where I get problems with MATLAB running out of
memory.

How would I calculate the Schur Complement implicitly? Is
there a special algorithm I could look up, or is an implicit
calculation of the Schur complement included in Matlab?

Thanks

Btw.

Nice book you wrote about sparse linear systems. I have to
admit I got lost pretty soon, as I'm only a mechanical
engineer and never learnt graphs, but it helped me a lot to
speed up my code for creating large sparse matrices.

"A B" <gitsnedbutzi@hotmail.com> wrote in message ...
> Hi Tim,
>
> Kbb could be larger than 5000, but 5000 is approximately the
> dimension where I get problems with MATLAB running out of
> memory.
>
> How would I calculate the Schur Complement implicitly? Is
> there a special algorithm I could look up, or is an implicit
> calculation of the Schur complement included in Matlab?
>
> Thanks
>
> Btw.
>
> Nice book you wrote about sparse linear systems. I have to
> admit I got lost pretty soon, as I'm only a mechanical
> engineer and never learnt graphs, but it helped me a lot to
> speed up my code for creating large sparse matrices.

Thanks - glad the book helped, even if just a bit. Yes,
it's got a lot of combinatorial stuff in it - more of that
than linear algebra. That's the nature of sparse matrix
computations. They're not so much linear algebra as graph
theoretic algorithms with reals that go along for the ride.

I'm not sure off-hand how you'd deal with the Schur
complement implicitly. The Sherman-Morrison update might be
used - I'm not sure.

How bad is it to let backslash handle the whole K? The
"static condensation" is exploited internally, in either a
multifrontal (LU) or supernodal (chol) way.

Consider this graph game. Take a symmetric matrix. Let
there be one node per row/col, and let there be an edge
(i,j) if A(i,j) is nonzero.

Now, pick a node, any node. Add edges between all of its
neighbors (they now form a clique), then eliminate the node
you picked (and its incident edges). Continue until the
graph is gone.

Every edge that ever appeared ... gives you the nonzero
pattern of L for the Cholesky factorization of L*L'=P*A*P'.
 The P gives you the order in which you eliminated the nodes.

Now, for a FEM matrix, your graph starts out as a collection
of cliques. Static condensation consists of eliminating
nodes whose neighbors already form a clique, so no new edges
are added.

Those cliques give LU, or Chol, a way of exploiting dense
matrix kernels (LU and CHOL both use the level-3 BLAS). If
you pre-eliminate via static condensation, your cliques are
smaller ... but there's less chance for using the BLAS.

So LU and CHOL naturally do "static condensation", if given
the right ordering.

See what happens if you do a CHOL with the nodes to
statically condense ordered first. You can do this with
AMD, or if you download CAMD from the File Exchange (in
SuiteSparse), you can do it more easily with that (give CAMD
a set of constraints to tell it how to order the nodes, with
the statically condensed nodes ordered first, and the rest
2nd, so that you have 2 constraint sets).

Then try CHOL with the resulting permutation as the
fill-reducing ordering.

The lower right corner of L will be the factor of your
statically-condensed matrix.

> Thanks - glad the book helped, even if just a bit. Yes,
> it's got a lot of combinatorial stuff in it - more of
that
> than linear algebra. That's the nature of sparse matrix
> computations. They're not so much linear algebra as
graph
> theoretic algorithms with reals that go along for the
ride.
>
> I'm not sure off-hand how you'd deal with the Schur
> complement implicitly. The Sherman-Morrison update
might be
> used - I'm not sure.
>
> How bad is it to let backslash handle the whole K? The
> "static condensation" is exploited internally, in either
a
> multifrontal (LU) or supernodal (chol) way.
>
> Consider this graph game. Take a symmetric matrix. Let
> there be one node per row/col, and let there be an edge
> (i,j) if A(i,j) is nonzero.
>
> Now, pick a node, any node. Add edges between all of its
> neighbors (they now form a clique), then eliminate the
node
> you picked (and its incident edges). Continue until the
> graph is gone.
>
> Every edge that ever appeared ... gives you the nonzero
> pattern of L for the Cholesky factorization of
L*L'=P*A*P'.
> The P gives you the order in which you eliminated the
nodes.
>
> Now, for a FEM matrix, your graph starts out as a
collection
> of cliques. Static condensation consists of eliminating
> nodes whose neighbors already form a clique, so no new
edges
> are added.
>
> Those cliques give LU, or Chol, a way of exploiting dense
> matrix kernels (LU and CHOL both use the level-3 BLAS).
If
> you pre-eliminate via static condensation, your cliques
are
> smaller ... but there's less chance for using the BLAS.
>
> So LU and CHOL naturally do "static condensation", if
given
> the right ordering.
>
> See what happens if you do a CHOL with the nodes to
> statically condense ordered first. You can do this with
> AMD, or if you download CAMD from the File Exchange (in
> SuiteSparse), you can do it more easily with that (give
CAMD
> a set of constraints to tell it how to order the nodes,
with
> the statically condensed nodes ordered first, and the
rest
> 2nd, so that you have 2 constraint sets).
>
> Then try CHOL with the resulting permutation as the
> fill-reducing ordering.
>
> The lower right corner of L will be the factor of your
> statically-condensed matrix.
>


Thanks Tim for your hints, I'll try the ordering of the
nodes and then use CHOL.

Can you suggest any introductory book on graphs and/or the
combinatorial stuff which would help me make your book
more understandable? Maybe there is some book for the
introduction to the analysis of sparse algorithms?

"A B" <gitsnedbutzi@hotmail.com> wrote in message ...

> Thanks Tim for your hints, I'll try the ordering of the
> nodes and then use CHOL.
>
> Can you suggest any introductory book on graphs and/or the
> combinatorial stuff which would help me make your book
> more understandable? Maybe there is some book for the
> introduction to the analysis of sparse algorithms?

Try using CAMD (in SuiteSparse on the File Exchange). It's
like AMD, except that you can give it a set of ordering
constraints. You can tell it to order all the "statically
condensed" nodes first, then order all the rest (2
constraint sets). Then within the sets, the approximate
minimum degree ordering will work its magic, to reduce
fill-in. Then use the resulting ordering to permute the
matrix for CHOL, or LU.

Off-hand, I'm not sure what book I could suggest, since it
would depend on what other background you have. I give a
*really* short review of the graph theory in my book, and I
cite some others there. I think Cormen, Leieserson, and
Rivest (Intro to Algorithms, I think the title is), MIT
Press, is a good book. That might be too advanced,
depending on your background.

The way it is done in SDT www.sdtools.com/sdt

is to use a factored matrix object for K_bb^-1 (very
critical if you want to have 100 000 DOFs). Then, you solve
for Kbb^-1*Kba by piece by block. Then you project the
matrices.

You have to realize that 1e5*5000*8/1023^3= 3.7 GB
If you want that to happen on a 32 bit machine you need to
go out of core (like SDT and ANSYS do).