[50G] Question

[Gilles]

(03-24-2015 06:41 AM)Tugdual Wrote:  

(03-23-2015 12:39 PM)Gilles Wrote:  Another way is to use the SOLVE command. I will post the source if I found some time...

You tickle my curiosity. Do you mean something like building and solving a system of equations by developing LU=A?

yes...

EDIT : this v2 allows parameters in the input matrix

Code:

LU2
«
 DUP SIZE HEAD → d     @ Empile A. Ordre de la matrice : d
 « d DUP
  « → I J 
   « CASE 
     'I==J' THEN 1 END
     'I<J'  THEN 0 END
     "l" I + J + S~N 
     END
   » 
  » LCXM               @ Initialse L

  d DUP
  « → I J 
   « CASE 
    'I>J' THEN 0 END
     "u" I + J + S~N  
     END
   » 
  » LCXM              @ Initialise U
  
  → L U
  «
   L U * LVAR UNROT
   AXL SWAP AXL 2. «=» DOLIST
   «EVAL» DOSUBS AXL  SWAP SOLVE 
   IF DUP {} == THEN DROP "No LU factorisation" DOERR END
   HEAD DUP
   
   L SWAP SUBST        @ renvoie L
   SWAP U SWAP SUBST   @ renvoie U
  » 
 » 
»

The code is from HPUserEdit
Exact mode, and flag -86 set (or 256 ATTACH)
To understand how it works, the debug mode is perfect.

Code:

 Examples :

[[ 2 -1 0 ]
 [ -1 2 -1 ]
 [ 0 -1 2 ]]
->

[[ 1 0 0 ]
 [ '-1/2' 1 0 ]
 [ 0 '-2/3' 1 ]]

[[ 2 -1 0 ]
 [ 0 '3/2' -1 ]
 [ 0 0 '4/3' ]]

--------------

[[ 1 5 8 0 ]
 [ 0 2 6 9 ]
 [ 1 5 11 7 ]
 [ 0 2 6 13 ]]

->

[[ 1 0 0 0 ]
 [ 0 1 0 0 ]
 [ 1 0 1 0 ]
 [ 0 1 0 1 ]]

[[ 1 5 8 0 ]
 [ 0 2 6 9 ]
 [ 0 0 3 7 ]
 [ 0 0 0 4 ]]

-------------
[[ 4 0 1 ]
 [ 2 1 0 ]
 [ 2 2 3 ]]

->

[[ 1 0 0 ]
 [ '1/2' 1 0 ]
 [ '1/2' 2 1 ]]

[[ 4 0 1 ]
 [ 0 1 '-1/2' ]
 [ 0 0 '7/2' ]]

--------------
[[4  4  3  ]
[3  -10  8 ]
[1  -4  'z'  ]]

->

[[ 1 0 0 ]
 [ '3/4' 1 0 ]
 [ '1/4' '5/13' 1 ]]

[[ 4 4 3 ]
 [ 0 -13 '23/4' ]
 [ 0 0 '(26*z-77)/26' ]]