--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolbox/comp.m	Thu May 07 18:36:24 2015 +0200
@@ -0,0 +1,41 @@
+function C = comp(A, k)
+%COMP    Comparison matrices.
+%        COMP(A) is DIAG(B) - TRIL(B,-1) - TRIU(B,1), where B = ABS(A).
+%        COMP(A, 1) is A with each diagonal element replaced by its
+%        absolute value, and each off-diagonal element replaced by minus
+%        the absolute value of the largest element in absolute value in
+%        its row.  However, if A is triangular COMP(A, 1) is too.
+%        COMP(A, 0) is the same as COMP(A).
+%        COMP(A) is often denoted by M(A) in the literature.
+
+%        Reference (e.g.):
+%        N.J. Higham, A survey of condition number estimation for
+%        triangular matrices, SIAM Review, 29 (1987), pp. 575-596.
+
+if nargin == 1, k = 0; end
+[m, n] = size(A);
+p = min(m, n);
+
+if k == 0
+
+% This code uses less temporary storage than the `high level' definition above.
+   C = -abs(A);
+   for j=1:p
+     C(j,j) = abs(A(j,j));
+   end
+
+elseif k == 1
+
+   C = A';
+   for j=1:p
+       C(k,k) = 0;
+   end
+   mx = max(abs(C));
+   C = -mx'*ones(1,n);
+   for j=1:p
+       C(j,j) = abs(A(j,j));
+   end
+   if all( A == tril(A) ), C = tril(C); end
+   if all( A == triu(A) ), C = triu(C); end
+
+end