how to apply a block of additional householder reflectors?

[Azua_Garcia__Giovann]

Hello,

I have the R part matrix of an A=QR decomposition e.g.
R =
1.1032 1.1223 1.1583
0 0.1705 -0.3269
0 0 0.0945
Now I add a row of ones to A and therefore to the old R:
Rnew =
1.1032 1.1223 1.1583
0 0.1705 -0.3269
0 0 0.0945
1 1 1
I can quickly compute the elementary reflectors (and taus) that can annihilate the ones below e.g.
V =
-2.5923 -26.7224 3.3802
tau =
0.2591 0.0028 0.1610
Now I would like to apply this compact form set of Householder reflectors to a vector b i.e. to compute c=Q'*b as part of a solve process and don't know how to do this. Basically I need to:

1) c=Q'*b this I do using LAPACKE_dormqr
2) c=H3*H2*H1*c this I don't know how to do given that I have V and tau.

The closest match I found was "dlarzb" but I see no tau as part of its input ... it also requires a T parameter "Array, DIMENSION (ldt,k). The triangular k-by-k matrix T in the representation of the block reflector" which I am not sure what it means. Is there such a function?

For example building and applying the reflectors in matlab would be:

% build first Householder reflector
v1 = [1 0 0 V(1)]';
H1 = eye(4) - tau(1)*(v1*v1')

% build second Householder reflector
v2 = [0 1 0 V(2)]';
H2 = eye(4) - tau(2)*(v2*v2')

% build third Householder reflector
v3 = [0 0 1 V(3)]';
H3 = eye(4) - tau(3)*(v3*v3')

c = H3*H2*H1*[[Q; zeros(1, 3)] I(:,4)]'*b

Thanks in advance,
Best regards,
Giovanni

[uhrm]

For updating a QR factorization, what one usually uses are Givens rotations. Your example would be updated like (Matlab pseudo-code):

R = [1.1032, 1.1223,  1.1583;
     0,      0.1705, -0.3269;
     0,      0     ,  0.0945];

A = [R; ones(1,3)];

G1 = planerot(A([1,4],1));
A([1,4],:) = G1*A([1,4],:)

G2 = planerot(A([2,4],2));
A([2,4],2:end) = G2*A([2,4],2:end)

G3 = planerot(A([3,4],3));
A([3,4],3:end) = G3*A([3,4],3:end)

In Matlab, see also "qrinsert" and "qrdelete".

Doing the same with MKL, you find the functions in BLAS. Creating the Givens rotation is called "rotg"; for applying it to two vectors you use "rot".

Hope this helps.

[Azua_Garcia__Giovann]

Hello uhrm,

Thank you for the hint :)

Best regards,
Giovanni