title 'IMLEQN: Solution of Systems of Equations';
  *--  ------- ---------------------------------;
proc iml;
   reset print log fuzz fw=5;
   %include iml(matlib);
*-- [1.] 3 EQUATIONS IN 3 UNKNOWNS (Consistent);
A= {1 1 -4, 1 -2 1, 1 1 1};
b= {2, 1, 0};
xx = t('X1' : 'X3');
print A '*' xx '=' b;
* they are consistent, since r(A) = r(A b);
print (r(A)) (r(A || b));
*-- A solution exists: use inv() or solve();
x = inv(A) * b;
x = solve(A,b);
print A ' * ' x '=' (A * x) '=' b;
*-- Echelon form of (A || b) shows solution;
r = echelon(A || b);

*-- [2.] 4 EQUATIONS IN 3 UNKNOWNS (Consistent);
A= {1 1 -4, 1 -2 1, 1 1 1, 2 -1 2};
b= {2, 1, 0, 1};
* they are consistent, since r(A) = r(A b);
print (r(A)) (r(A || b));
*-- A solution exists, but more equations than unknowns,
    so use ginv;
x = ginv(A) * b;
print A ' * ' x '=' (A * x) '=' b;
*-- Echelon form of (A || b) shows solution;
r = echelon(A || b);

*-- [3.] 3 EQUATIONS IN 3 UNKNOWNS (Inconsistent);
A={1 3 1, 1 -2 -2, 2 1 -1};
b={2,3,6};
*-- inconsistent, since r(A) < r(A b);
print (r(A)) (r(A || b));
r = echelon(A || b);
*-- inv() fails, since A is singular;
x = inv(A) * b;
* g-inverse provides approx. solution;
x = ginv(A) * b;
print A '*' x ' = ' (A * x);
*-- change rhs to give consistent set;
b = {2,3,5};
r = r(A || b);
r = echelon(A || b);
*-- since r(A) < ncol(a) solution is not unique;
*   but ginv finds a solution;
x = ginv(A) * b;
print A '*' x ' = ' (A * x) '=' b;
*-- Equivalent way of producing a solution;
A11 = A[{1 2},{1 2}];
A12 = A[{1 2}, 3];
b1  = b[{1 2},];
x2  = -1.1;
* solve for 2 independent unknowns;
x1 = inv(A11) * b1 - inv(A11) * A12 * x2;
quit;