#% FILE antonNM.mpl

interface(echo=2);
# Example 2 - AR6 Section 9.2 p442; AR7 Section 9.9 p512
# CAUTION. Anton makes U unit upper triangular;
# below we make L unit lower triangular.
 with(linalg):
# Find the LU decomposition of the matrix
 A:= matrix([ [ 2,  6,  2],
               [-3, -8,  0],
               [ 4,  9,  2] ]):
# Set the diagonal entries of  L  equal to 1 and  U  equal to A
 L:= diag(1,1,1):                               U:= A:
#
# Now reduce  U  to row-echelon form, storing the INVERSES of the row
# operations used in the matrix  L



 L[2,1]:= U[2,1]/U[1,1]:               U:= addrow( U,1,2,-L[2,1]):
  L[3,1]:= U[3,1]/U[1,1]:               U:= addrow( U,1,3,-L[3,1]):
  L[3,2]:= U[3,2]/U[2,2]:               U:= addrow( U,2,3,-L[3,2]):
  print( L,U );

# Verify that    A = L*U
 equal( A,evalm(L&*U) );



# Exercise 18 - AR6 Section 9.2 p449; AR7 Section 9.9 p519
# An example of the  PA = LU  decomposition of a matrix
 with(linalg):
 A:= matrix([ [3,-1, 0],
               [3,-1, 1],
               [0, 2, 1] ]):
# The following row operations will reduce  A  to Echelon form:
L:= diag(1,1,1):                               U:= swaprow(A,2,3):
L[3,1]:= U[3,1]/U[1,1]:               U:= addrow( U,1,3,-L[3,1]);
# Alternatively, perform the above row interchanges on the identity matrix
 L:= diag(1,1,1): P23:= swaprow(L,2,3);
# Now perform an LU decomposition on PA, as in LU Example above
 PA:= evalm( P23&*A );
# Looking at the result on the screen, and working the steps
# one gets to L
   L:=matrix([[1,0,0],[0,1,0],[1,0,1]]):
   equal(A,P23 &* L &* U);
   equal( PA,evalm(L&*U) );


 A:= array( symmetric,1..3,1..3 ):
 with(linalg):
 equal( A,transpose(A) );

interface(echo=1):

