HP-71 matrix methods (LU decomp?)



#2

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 ]
[ 3 4 ]

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 ]
[ 0 0 ] [ 1 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:


PA = LU => A = P-1LU

det(A) = det(P-1) * det(L) * det(U) = (-1 or 1) * 1 * det(U)


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


A-1 = (P-1LU)-1 = U-1L-1P

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.


For the matrix

A = [ 1  2 ]
[ 3 4 ]

the following can be extracted from the HP-15C calculation of determinant or system solution:

P = [ 0  1 ]  L = [ 1             0 ]  U = [ 3               4 ]
[ 1 0 ] [ 0.3333333333 1 ] [ 0 0.6666666668 ]

The "LU" command of the HP-48G returns the following:

P = [ 0  1 ]  L = [ 3               0 ]  U = [ 1  1.33333333333 ]
[ 1 0 ] [ 1 0.666666666667 ] [ 0 1 ]

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 = 

exact, HP-48G HP-71B, HP-28C, HP-42S, HP-48S(?) HP-15C
[ -2 1 ] [ -1.99999999998 0.999999999994 ] [ -2 0.9999999999 ]
[ 1.5 -0.5 ] [ 1.49999999999 -0.499999999997 ] [ 1.5 -0.5 ]

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
(3. * .66666666667), as well as the slightly-inexact elements of the inverse.

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 ]  
[ 0 0 ]
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:
The HP48G and its successors have matrix arithmetic that was reworked by Paul McClellan to use 15 digit arithmetic for internal operations, rounding to 12 digits only occurring after everything is done.

The HP71 and HP48S rounded internal operations to 12 digits at intermediate steps in the carrying out of matrix operations.


Further insights, anyone?

-- KS


Edited: 7 Sept 2010, 4:58 a.m. after one or more responses were posted


#3

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.


#4

Quote:
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.

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.


#5

Quote:
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?

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:
I would doubt that any official HP product for the HP-71B would address these minor issues.

Agreed.

#6

Quote:
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.

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

#7

- the HP48G and successors get the inverse right thanks to the 15 digits used throughout.
- the 42S calculates the determinant of a 2x2 matrix directly. If you calculate the determinant of
[[ 1 2 0 ]
[ 3 4 0 ]
[ 0 0 1 ]],
it is no longer exactly 2.
- the 48G series and successors have an additional test: if the matrix contains only integers, the determinant returned will be rounded to the nearest integer.

#8

Werner --

Quote:
- the HP48G and successors get the inverse right thanks to the 15 digits used throughout.

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:
- the 42S calculates the determinant of a 2x2 matrix directly. If you calculate the determinant of
[[ 1 2 0 ]
[ 3 4 0 ]
[ 0 0 1 ]],
it is no longer exactly 2.

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.


Possibly Related Threads…
Thread Author Replies Views Last Post
  AFTER HP-Prime update, Shift+Matrix CRASHES Joseph Ec 3 2,108 12-06-2013, 11:06 AM
Last Post: Joseph Ec
  HP Prime Matrix TERRIBLE bug and question uklo 19 5,357 11-25-2013, 12:10 PM
Last Post: Mic
  HP Prime - Reset methods bluesun08 3 3,338 11-21-2013, 09:49 AM
Last Post: Joe Horn
  HP Prime: editing a matrix Alberto Candel 6 2,423 11-20-2013, 06:26 PM
Last Post: Helge Gabert
  Absolute Value and Matrix BruceTTT 5 2,197 11-11-2013, 11:52 PM
Last Post: Walter B
  How to move lexfiles from PC to 71 w/o HP-IL? Joe Horn 9 3,391 10-18-2013, 03:50 PM
Last Post: Christoph Giesselink
  WP-34S Matrix operations with routine-local registers? Tom Grydeland 1 1,230 09-04-2013, 10:46 AM
Last Post: Marcus von Cube, Germany
  Matrix Characteristic Polynomial - Reloaded. Ángel Martin 12 3,629 08-22-2013, 05:33 PM
Last Post: Thomas Klemm
  Matrix Richard Berler 3 1,497 08-18-2013, 06:24 PM
Last Post: Paul Dale
  Advantage/CCD Matrix Challenge Ángel Martin 1 1,211 08-09-2013, 06:22 PM
Last Post: Thomas Klemm

Forum Jump: