HP-71 matrix methods (LU decomp?) - Printable Version +- HP Forums (https://archived.hpcalc.org/museumforum) +-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html) +--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html) +--- Thread: HP-71 matrix methods (LU decomp?) (/thread-171592.html) |
HP-71 matrix methods (LU decomp?) - Karl Schneider - 09-05-2010 A recent post noted a bug in the HP-71B (ROM 1BBBB and earlier) in which some leap years were incorrectly identified. J-F Garnier responded that it was a known bug, but did not cite his own web page listing them: http://www.jeffcalc.hp41.eu/emu71/bug71.html At the bottom of the page are listed some inexact results of matrix arithmetic using the HP-71 Math ROM. E.g., the calculated determinant of the matrix
[ 1 2 ] is not quite an integer (-2.00000000001 instead of -2), and the elements of the calculated inverse similarly are not quite exact. Also, the calculated determinant of the singular matrices
[ 0 1 ] and [ 0 0 ] are returned respectively as 0 with a message "WRN: Overflow", and -0 (with no message). I suspect that the HP-71B (1983) and its Math Pac ROM (1984) utilizes linear-algebra methods similar to those of the 1982 HP-15C (with its 10-digit Coconut processor), adapted for the 12-digit Saturn processor. These methods seem to have been refined further for the Saturn-based HP-28C (1986), HP-42S (1988), HP-48S (1990), and HP-48G (1993), which do not exhibit all of these minor flaws. Valentin Albillo has explained that the functionality contained in the HP-71 Math ROM was removed from the base HP-71 ROM in order to make room for other features (i.e., support for the IEEE floating-point standard that was then a proposal and was later adopted). So, time might have been short for incorporating all refinements in the HP-15C to the HP-71's new Saturn processor. For reasons of RAM size and execution speed, the HP-15C performs an LU decomposition of a square matrix A (larger than 1 x 1) prior to calculating determinants, inverses, and system solutions. The Doolittle algorithm is used, producing a lower-triangular unit matrix L having ones on the main diagonal, while the upper-triangular matrix U has non-zero values. The permutation matrix P performs any necessary row transpositions, and is an identity matrix if none are required:
where det (U) is the product of its main-diagonal entries, and det (P-1) is odd for an odd number of row transpositions.
It follows that
On the HP-15C, L and U are composed "in place", and row-swap information is stored in spare bits (no P matrix is stored). L and U can be inverted in place, and the resulting inverses can also be multiplied together in place. In the case of a singular or near-singular input matrix, the HP-15C will also modify elements slightly to allow the necessary LU decomposition without divide-by-zero or overflow. The Saturn-based models use a Crout algorithm, which produces an upper-triangular unit matrix U. Otherwise, it appears that all of these models take the same basic approach as that of the HP-15C. The HP-71B's methods are stated on pp. 86-88 of the Math Pac Owner's Manual.
A = [ 1 2 ] the following can be extracted from the HP-15C calculation of determinant or system solution:
P = [ 0 1 ] L = [ 1 0 ] U = [ 3 4 ] The "LU" command of the HP-48G returns the following:
P = [ 0 1 ] L = [ 3 0 ] U = [ 1 1.33333333333 ]
It is easy to see where these methods may produce non-integer or inexact values in determinants and inverses. However, the HP-15C, HP-42S, and HP-48G did produce 'cleaner' results, as listed below: A-1 = As J-F Garnier pointed out, and as is documented on page 88 of the HP-71B Math Pac Owner's Manual, inverses can be computed more accurately by using the linear-system command MAT X = SYS(A,I) with an identity matrix as I. In this case, an exact inverse is returned by the HP-71B (but not by the HP-42S "SIMQ" command). It seems possible that some algorithmic refinements in the HP-15C did not get incorporated into the HP-71 Math ROM. For example: 1. Not adjusting elements of the singular matrix A might explain the "overflow" messages that likely occurred during LU decomposition.
2. Some inexactitude in the calculations using elements of the Crout-algorithm L and U matrices might explain the determinant of 2.0000000001 The HP-28C, HP-42S, and HP-48G return integer values for the determinants, and in the process do not warn of overflows or return "-0". However, the HP-28C and HP-42S return the very same inexact matrix-inverse values as the HP-71B does, while the HP-48G returns the exact values.
The HP-48G does return "LU error: Overflow" for LU decomposition of the singular matrix [ 0 1 ]apparently becuase no non-zero element is present in the left column. However, a valid Crout LU decomposition can apparently be produced; please see this Forum thread from 2006. A comment by Rodger Rosenbaum in a thread from 2007:
Quote: Further insights, anyone? -- KS
Edited: 7 Sept 2010, 4:58 a.m. after one or more responses were posted
Re: HP-71 matrix methods (LU decomp?) - Egan Ford - 09-05-2010 There is a newer MATH ROM I can send you. More info here: http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv017.cgi?read=112732.
Re: HP-71 matrix methods (LU decomp?) - Karl Schneider - 09-05-2010 Quote: Hello, Egan -- That seemed to be an extremely generous offer, given that I'm using a physical HP-71B with Math ROM version 1A. (I don't know if there was a later version.) I read the archived thread you posted, but still don't quite understand the "MATHROM" product that is recorded on EPROM. Is it a hardware item? Is it designed to plug into the HP-71B, or must it be loaded via HP-IL or other means? Or, is it an emulator for the HP-48GX or HP-49g+/50g? (I have neither.) A comment by Rodger Rosenbaum in the same 2007 thread you cited indicates that further refinement was done for the HP-48G, so I would doubt that any official HP product for the HP-71B would address these minor issues. My main objective in starting this new thread was to explain the HP-71 results listed by J-F Garnier in his bug list, and to re-introduce HP matrix methods and LU decomposition. Thanks, -- KS
Edited: 6 Sept 2010, 4:18 a.m.
Re: HP-71 matrix methods (LU decomp?) - Werner - 09-06-2010 - the HP48G and successors get the inverse right thanks to the 15 digits used throughout. Re: HP-71 matrix methods (LU decomp?) - Egan Ford - 09-06-2010 Quote:I do not have a MATH ROM for the 71B, so I opted to have one burned to an EPROM. The individual that did this for me had a newer unofficial version of the MATH ROM. It is discussed in the aforementioned thread (last post).
I can send you the LEX file. If you have enough RAM and HP-IL you can transfer it into memory. To avoid confusion you may want to remove your MATH ROM first. Or I can send you an HDRIVE1.DAT with MATHROM.LEX for testing with EMU71. Quote:Agreed. Re: HP-71 matrix methods (LU decomp?) - Karl Schneider - 09-09-2010 Quote: Hello again, Egan -- Thank you. I do have a physical HP-71B with 36 kB added RAM and HP-IL. However, I do not have Jean-Francois' PIL-Box HPIL-USB device, HP's HPIL-to-ISA bus device, or Emu71 software. I'd have to obtain some hardware before being able to put MATHROM.LEX to use. I'll send an e-mail before too long.
-- KS
Re: HP-71 matrix methods (LU decomp?) - Karl Schneider - 09-09-2010 Werner --
Quote: Yes. Incidentally, I added that piece of info independently as an edit to my original post, after Egan provided a link to an archived thread.
Quote: That makes sense, as the determinant of a 2 x 2 matrix is more easily calculated in a straightforward manner than via LU decomposition. Incidentally, for a 1 x 1 matrix, the HP-15C still produces an (identical) 1 x 1 result matrix with the determinant, but doesn't identify it as an LU decomposition (e.g., "C--"). All things considered, that is the best approach -- and another example of attention to detail in the HP-15C. -- KS
Edited: 9 Sept 2010, 4:35 a.m.
|