hssrank: now it computes the maximum rank of the hss blocks

Author: massei

@hss/hodlrrank.m

@@ -0,0 +1,3 @@
+function k = hodlrrank(A)
+k = hss_hodlr_rank(A);
+end

@hss/hssrank.m

@@ -1,3 +1,4 @@
 function k = hssrank(A)
-k = hss_rank(A);
+% Maximum rank of all the HSS block rows at each level
+	k = hss_rank(A);
 end

@hss/private/hss_chol_fact.m

@@ -43,7 +43,7 @@
     k = size(A.U, 2);
     [row, col] = size(A.D);
     if k == row
-        x = zeros(col, size(b,2)); D = A.D; U = A.U; ind = zeros(1,0); cind = [1:col];
+        D = A.D; U = A.U; ind = zeros(1,0); cind = [1:col];
         return
     end
     ind = [1:row-k];  % Indices of computed entries in the solution, at the end

@hss/private/hss_chol_solve.m

@@ -4,7 +4,7 @@
 %	       The procedure is based on the implicit generalized Cholesky factorization
 %	       of A and is taken from the paper:
 %
-%	       Xia, J., Chandrasekaran, S., Gu, M., & Li, X. S. (2010). Fast algorithodlrs for
+%	       Xia, J., Chandrasekaran, S., Gu, M., & Li, X. S. (2010). Fast algorithms for
 %              hierarchically semiseparable matrices. Numerical Linear Algebra with Applications, 17(6), 953-976.
 %
 if A.leafnode == 1

@hss/private/hss_hodlr_rank.m

@@ -0,0 +1,8 @@
+function k = hss_hodlr_rank(A)
+% maximum rank of the off-diagonal blocks
+	if A.leafnode == 1
+		k = 0;
+	else
+		k = max([rank(A.B12), rank(A.B21), hss_hodlr_rank(A.A11), hss_hodlr_rank(A.A22)]);
+	end
+end

@hss/private/hss_rank.m

@@ -1,7 +1,23 @@
 function k = hss_rank(A)
-if A.leafnode == 1
-    k = 0;
-else
-    k = max([rank(A.B12), rank(A.B21), hss_rank(A.A11), hss_rank(A.A22)]);
+	k = hss_rank_rec(A, zeros(size(A.B12, 1), 0), zeros(size(A.B21, 1), 0));
 end
+function k = hss_rank_rec(A, M1, M2)
+% Compute the rank of all the HSS block rows at each level
+	if A.leafnode == 1
+		k = 0;
+	else
+		k1 = rank([A.B12, M1]);
+		k2 = rank([A.B21, M2]);
+		if A.A11.leafnode == 0
+			k3 = hss_rank_rec(A.A11, A.A11.Rl * [A.B12, M1], A.A11.Rr * [A.B12, M1]);
+		else
+			k3 = 0;
+		end
+		if A.A22.leafnode == 0	
+			k4 = hss_rank_rec(A.A22, A.A22.Rl * [A.B21, M2], A.A22.Rr * [A.B21, M2]);
+		else
+			k4 = 0; 		
+		end	
+		k = max([k1, k2, k3, k4]);
+	end	
 end

@hss/spy.m

@@ -16,7 +16,7 @@ function spy(H)
 axis([ 0 n 0 m ]);
 hold on;
 
-spy_draw_block(H, 0, 0, hssrank(H));
+spy_draw_block(H, 0, 0, hodlrrank(H));
 
 hold off;
 
@@ -65,4 +65,4 @@ function spy_draw_block(H, xoffset, yoffset, rk)
         sprintf('%d', thisrk), ...
         'HorizontalAlignment', 'center');
 end
-end
+end

rk_sylv.m

@@ -123,7 +123,7 @@
         check_residual = it >= needed_iterations;
     end
 
-    % check_residual = true;
+     %check_residual = true;
 
     if check_residual
         if poleA ~= inf