gnucap-devel
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Query regarding bsmatrix


From: al davis
Subject: Re: Query regarding bsmatrix
Date: Thu, 20 Jul 2023 16:11:16 -0400


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.




reply via email to

[Prev in Thread] Current Thread [Next in Thread]