Only sorted access to entries in HSS form

[bonevbs]

Hello again, sorry if this is a trivial question - I wasn't sure where to ask this. I ran across the issue that individual entries in HSS form A(I,j) can only be accessed if the index sets i and j are sorted. So I had to resort to a wrapper in the form of

function out = sorted_access(Afun, i, j)
  [~,is] = sort(i); [~,js] = sort(j);
  out(is,js) = Afun(i(is), j(js));
end

Edit: After diving into subs ref I do see now how unsorted indices can be an issue. I'll leave the snippet here in case it helps someone else.

I am sorry to spam like this but I'd like to also ask how to get both generators for off diagonal blocks. In the current form calling

[U,V] = hssA.hss_offdiag('upper')

causes the error Output argument "varargout{2}" (and maybe others) not assigned during call to "hss/subsref". which seems to be a matlab related issue. It seems Matlab only supports one output with the current implementation. Any hints on how I could resolve this?

[massei]

Hi Boris! No worries for the spam :), the discussion is useful. For subsref with ordered indices: we wanted subsref to return a sub matrix in the hss format and if the indices are not sorted can be subtle. if you just want the output in the dense format: we overload the function "get" (basically, we copied your snipplet) so that you can call get(A, I, J) to get A(I, J) even I and J are not ordered. The second issue is caused by the fact that hss_offdiag is in the private repository. If you add a file offdiag.m in the @hss repository containing

function [U, V] = offdiag(H, ul) [U, V] = hss_offdiag(H, ul); end

then you can just call offdiag(A, 'upper') for your purposes.

[bonevbs]

Hi Stefano, thanks a lot for your quick reply, it's very helpful!

In fact one other thing I noticed is that there seems to be an issue when accessing submatrices of HSS matrices with non-standard blocking. I guess non-standard blocking is still experimental?

So, when I run full(hssA(I,J)) and hssA is such a matrix I get a dimension mismatch. If you are interested I can replicate the matrix and post it here.

[massei]

Yes, if you can post your example that would be useful.

[bonevbs]

schur.mat.zip Ok I just saved them into a file. Running

load('schur.mat')
full(S(I,J))

should give you the error Error using * Incorrect dimensions for matrix multiplication..

Btw, I checked out the latest commit. It seems you forgot to add get? I am not sure, at least running get(S,I,J) gives me Too many input arguments.

[massei]

Thank you for spotting this! There was a mistake in the transpose for hss matrices with non standard partitioning. Now it should be fixed.

[bonevbs]

Thanks a lot!! :) And glad I can help.. in fact I might have found some more I'm afraid: I want to invert the topmost HSS block and inv(hssA.A11) should work right in my understanding. But even if I am trying to invert the whole matrix with non-standard clustering it gives me an error when I do inv(hA). I have attached the matrix. inv.mat.zip

[massei]

The problem with using hssA.A11 directly is that there are some fileds (for instance topnode) that have to be updated before applying some hss operations. If you want to do such operation you need to clean the exactred submatrix with a function like the following

function H = clean_hss(H) H.topnode = 1; H = clean_ric(H); end

function H = clean_ric(H) if isempty(H.A11) && isempty(H.A22) H.leafnode = 1; H.Wl = []; H.Wr = []; H.Rl = []; H.Rr = []; H.ml = []; H.nl = []; H.mr = []; H.nr = []; else H.leafnode = 0; H.U = []; H.V = []; if ~isempty(H.A11) H.A11 = clean_ric(H.A11); end if ~isempty(H.A22) H.A22 = clean_ric(H.A22); end end end

[bonevbs]

I see! I wasn't aware that there is more to fix than topside = 1. This should cover the things I am trying to do right now, so thanks for all the helpful answers! :)

I don't know whether you did - you might want to have a look at the matrix that I attached anyway. I was trying to invert the entire matrix, not just the top block and it didn't work for some reason. Again it's a matrix with non-standard blocks, but maybe I am doing something wrong?

[massei]

Hi Boris, the problem should be fixed, it was affecting backslash for HSS with non standard partitioning. Let us know if it works now.

[bonevbs]

Hi Stefano, thanks a lot I will have a look and let you know! Boris