% discussed pivoting some more, including the difference between partial
% pivoting and scaled partial pivoting, then looked at the construction of an
% LU factorization for some suitably permuted version PA of A, storing
% the multipliers used in the places that supposedly contain zero entries
% because of elimination. The multipliers form the strictly lower-triangular
% part of the unit lower triangular factor L, while the upper triangular matrix
% U comprises the pivot rows used. 
% Also re-iterated that, after the LU factorization
% is found, one solves A?=b in two steps, first solving L? = Pb (by forward
% substitution), getting some solution y, then solving U?=y (by backward
% substitution), getting the solution x.  
% 
% Also pointed out that the textbook erroneously refers to matlab's LU command
% which is, actually, the lu command.
% To make this clear, did
>> help lu

 LU     LU factorization.
    [L,U] = LU(X) stores an upper triangular matrix in U and a
    "psychologically lower triangular matrix" (i.e. a product
    of lower triangular and permutation matrices) in L, so
    that X = L*U.  X must be square.
 
    [L,U,P] = LU(X) returns lower triangular matrix L, upper
    triangular matrix U, and permutation matrix P so that
    P*X = L*U.
 
    LU(X), with one output argument, returns the output from
    LINPACK'S ZGEFA routine.
 
    LU(X,THRESH) controls pivoting in sparse matrices, where
    THRESH is a pivot threshold in [0,1].  Pivoting occurs
    when the diagonal entry in a column has magnitude less than
    THRESH times the magnitude of any sub-diagonal entry in that
    column.  THRESH = 0 forces diagonal pivoting.  THRESH = 1 is
    the default.
 
    See also LUINC, QR, RREF.

>> A = [.01 1
>>      1   1];
>> [l,u,p] = LU(A)
??? [l,u,p] = LU
              |
Undefined function or improper matrix assignment.

>> which LU
C:\MATLABR11\toolbox\matlab\matfun\lu.m
>> which lu
lu is a built-in function.
>> [L,U,P] = lu(A)
L =
    1.0000         0
    0.0100    1.0000
U =
    1.0000    1.0000
         0    0.9900
P =
     0     1
     1     0

% note that, somewhat wastefully, the record of the row permutations used is
% returned as a full matrix, P , when it would have sufficed to return the
% n-vector piv actually generated for which  Pb = b(piv) for any n-vector b.
  
% Next, I was asked what would happen if A is not invertible.
% Here is the A we just used:
>> A =
    0.0100    1.0000
    1.0000    1.0000
% we make it noninvertible by changing its first entry this way:

>> A(1,1) =1
A =
     1     1
     1     1
[L,U,P] = lu(A)
L =
     1     0
     1     1
U =
     1     1
     0     0
P =
     1     0
     0     1
% So, it seems that the factorization proceeds undisturbed, corresponding to the
% the fact that if we cannot find a pivot row for some unknown, then there is
% nothing to eliminate. 
% The difficulty comes when we want to back substitute, i.e., solve U?=y.

% This is also evident when we try to solve A?=b directly:
>> A\[1;1]
Warning: Matrix is singular to working precision.
ans =
   Inf
   Inf
% What happens when we make A almost singular?
>> A(1,1) = 1+eps
A =
    1.0000    1.0000
    1.0000    1.0000
>> A\[1;1]
Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 5.551115e-017.
ans =
     0
     1
% i.e., we get the more complicated message about some rcond being very small.
% this will be discussed in the next lecture.
>> diary off