Full matrix inversion for triangular and

[John W. Eaton]

On  6-Dec-2006, David Bateman wrote:

| The attached patch adds probing of the matrix type and special casing of
| positive definite and triangular matrices to the Matrix and
| ComplexMatrix classes inverse methods and adjusts the Finv and xpow
| functions appropriately. As an example consider
| 
| N=1000;
| A = 10*eye(N) + randn(N,N);
| t0 = cputime(); BA = inv(A); t(1) = cputime() - t0; ## Normal case
| P = A*A';
| t0 = cputime(); BP = inv(P); t(2) = cputime() - t0; ## PD case
| U = triu(A);
| t0 = cputime(); BU = inv(U); t(3) = cputime() - t0; ## Upper Triangular
| L = tril(A);
| t0 = cputime(); BL = inv(L); t(4) = cputime() - t0; ## Lower Triangular
| t
| 
| which returns
| 
|    7.7648   7.7618   3.8974   4.0144
| 
| for 2.9.9, and for 2.9.9+ with the attached patch returns
| 
|    7.8098   3.7634   1.3958   1.4018
| 
| for significant gains in speed. Note that this essentially makes cholinv
| redundant, though it appears that the same is not so for chol2inv. Consider
| 
| t0 = cputime(); BP = cholinv(P); cputime() - t0
| ans =  3.7494
| 
| R = chol(P); t0 = cputime(); BP = chol2inv(R); cputime() - t0
| ans =  2.6386
| t0 = cputime(); BP = inv(R); BP = BP * BP'; cputime() - t0
| ans =  5.3012
| 
| Now to check correctness, I did
| 
| N = 100;
| n = 10;
| for i=1:n
|   A = 10*eye(N) + randn(N,N);
|   P = A*A';
|   assert(A*inv(A), eye(N), 1e-10);
| endfor
| 
| for i=1:n
|   A = 10*eye(N) + randn(N,N);
|   U = triu(A);
|   assert(U*inv(U), eye(N), 1e-10);
| endfor
| 
| for i=1:n
|   A = 10*eye(N) + randn(N,N);
|   L = tril(A);
|   assert(L*inv(L), eye(N), 1e-10);
| endfor
| 
| which ran through without assert signaling any failures. This patch
| excludes the toeplitz solver code I previously sent

OK, please check in these changes.

Thanks!

jwe