48/50 Linear System Solver

[Han]

This was a question that I stumbled upon while working on the solver an equation library program. Solve for x in Ax=b if:

Let A = [[12.25, 1.25, 4.5],[1.25,12.25,4.5],[4.5,4.5,3]]
Let b = [[-166.2],[157.2],[-3]]

Note that A is singular. On the HP48, with flag -22 SET (infinite -> MAXREAL), I get:

b A / -> x=[[-25.0333...],[4.3666...],[30]]
b A LSQ -> x=[[-14.972727...],[14.42727...],[-.1818...]]

The b A LSQ method gives a solution with smaller norm than the b A / method. I am able to compute the same solution using SVD decomposition on A and applying the pseudo-inverse.

Code:

<<
A SVD
OBJ-> DROP DROP @ get nonzero eigenvalues
{} + + INV @ turn into a list and compute reciprocals
0 + EVAL ->V3 { 3 3 } DIAG-> @reconstruct inverse 3x3 diagonal matrix
SWAP TRN SWAP * SWAP TRN * @ compute the pseudo inverse
b *
>>

However, I am not sure how the HP48 got the answer for b A / ... Comparing the norms of the residuals, b A / gives a smaller norm for Ax-b, but clearly the norm of the solution is larger than that from b A LSQ. So what algorithm does b A / apply? (Would really like to not have to break out Jazz and decompile some ROM routines just to figure out the algorithm...)

EDIT: Additionally, RSD according the the advanced user's guide is supposed to compute Ax-b. However, I seem to get different answers when I manually compute Ax-b as compared to using RSD. Is this the case for anyone else? What am I overlooking?