- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

Link Copied

2 Replies

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Hello uhrm,

Thank you for the hint :)

Best regards,

Giovanni

Thank you for the hint :)

Best regards,

Giovanni

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page