Re: Query regarding bsmatrix

[Jacob M]

Hello Al,
Thanks for getting back to me.  What you're saying makes a great deal of
sense.  A couple of things:

"To go along with this, the matrix is stored by row below the diagonal,
and by column above."

Yes, it's ingenious.  I find it easier to picture this matrix as a stack of
reversed 'L' shapes, lining up along the diagonal.  Since circuit matrices
have a roughly symmetrical sparsity pattern, it's an efficient storage
scheme.  And, as you point out, means that the LU decomposition routines
benefit from very good cache locality.

The one potential problem is that, in some cases, the 'arms' of the L
shapes can become very long.  I don't know how much of a problem that is in
the real world.  Two questions: is there an easy way to print out the
matrices?  And do you think it'd be beneficial to implement a reordering
scheme like reverse Cuthill-McKee?  (There's an open source implementation
in Boost, I believe.) I see that the code currently reverses the node
order, but surely that can't be close to optimal?

"So this variant here is designed so you can make a small number of
updates to the "A" matrix, leaving most of it as it was, then only
recalculate the few spots in LU that need to be recalculated, leaving
most of it as it was."

Just an observation: even if you were doing a full LU decomposition every
time, it still makes sense to maintain separate A and LU matrices as it
means that you can be 'lazy' about stamping only the elements that change
into the matrix.  Performing LU decomposition in place, as Spice, Qucs and,
I think, Xyce do it is not a good strategy, even though it uses half the
memory.

One thing puzzles me, and I can't find an immediate answer in the code:
when does it become preferable to do a full decomposition?  That is to say,
is there a percentage of changed values above which it's more efficient to
start over from scratch?  And at what point does accumulated rounding error
become significant?

Also, comparing the two lu_decomp functions:

 u(bn,mm) /= d(bn,bn);

becomes:

u(bn,mm) = aa.u(bn,mm) / d(bn,bn);

Okay, this makes sense if the diagonal is unchanged, but what if a cell on
the diagonal changes?  Don't the stored values in U need to be rescaled?
(Don't get me wrong, I'm not trying to pick holes, it's just that this
implementation of incremental update to LU seems far simpler than the
literature I've been able to find suggests it ought to be, and I'm trying
to understand why that is.)

"That is what "prop" does.  It keeps track of
what else needs to be recalculated when you make a spot change in A."

Ok, thanks - I think I understand better.  Just to make sure I understand
the logic - if you had a 'pinch point' at which there's JUST a diagonal,
with no arms on the 'L', it would act as a brick wall preventing earlier
calculations from influencing anything further on.  Is that right?

And I'm still trying to understand this:

prop = mm;
if (bn < mm) {
    prop = mm;

It would make sense, I think, if the third line was "prop = bn", but as
things are, prop is going to be set to mm whether or not bn < mm, so the
second assignment isn't actually doing anything (and is doubtless being
removed by the compiler).  If it's fine to remove the second assignment,
it'd make the code a lot easier to follow. :-)  On the other hand, was the
intent to set prop = std::min(mm, bn)?

Regarding licensing of this and certain other files:

"Nobody really knows where it came from, so nobody really has the right to
claim ownership of it, even if this version is a little different from
another."

I appreciate the clarification.  Thank you.

Best wishes,
Jacob

On Thu, Jul 20, 2023 at 9:12 PM al davis <ad211@freeelectron.net> wrote:

>
>
> To understand BSMATRIX you must first understand classic LU
> decomposition for a full matrix.  I will not repeat that here, but
> rather refer you to any numerical analysis text.
>
> The classic LU decomposition has 3 nested loops ...
> .. down the diagonal
> .... by row
> ...... by column
>
> These can be in any order, in theory with identical results and run
> time.  Since Fortran stores by column, for Fortran code, that way works
> best, because down a column means adjacent spots in memory.
>
> In practice, number of math-operations (multiply and add) is identical,
> but the time it takes to find a particular spot in the matrix varies
> somewhat.
>
> This led to accusations that C is inherently slower because the usual
> storage is by row, making the "by column" loop slower because of slower
> indexing and cache misses.  This is easily solved by:
> .. down the diagonal
> .... by column
> ...... by row
>
> The code here uses a variant based on the Crout or Doolittle algorithm
> ,, down the diagonal
> --- if above diagonal
> .... by row
> ...... by column
> --- else (below diagonal)
> .... by column
> ...... by row
> ---
>
> To go along with this, the matrix is stored by row below the diagonal,
> and by column above.  This gives us two advantages ....
> 1. In the innermost loop, memory access is along a vector.  Every spot
> is adjacent to the previous.
> 2. The calculation of a particular spot in L or U is encapsulated.  The
> "Lvalue" in the innermost loop does not move.
>
> Aside from that, it's standard LU decomposition.
>
> In circuit simulation, for large circuits, it is common for parts of
> the circuit to converge faster than other parts.  In multi-rate, even
> more so.
>
> So this variant here is designed so you can make a small number of
> updates to the "A" matrix, leaving most of it as it was, then only
> recalculate the few spots in LU that need to be recalculated, leaving
> most of it as it was.  That is what "prop" does.  It keeps track of
> what else needs to be recalculated when you make a spot change in A.
>
>
> On the licensing .....
>
> That manual is very old and not maintained any more.  It's all in the
> wiki now.  I had forgotten about that statement, but  basically it says
> that those files are not really original anyway, and I don't remember
> where they came from.  That's what "public domain" really is.  Nobody
> really knows where it came from, so nobody really has the right to
> claim ownership of it, even if this version is a little different from
> another.
>
>