Re: Query regarding bsmatrix
al davis via Gnucap-devel <[email protected]> Thu, 20 Jul 2023 16:11:16 -0400
| Newsgroups | gmane.comp.gnu.gnucap.devel |
|---|---|
| Message-ID | <20230720161116.5f67bf71@s> |
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.