function [AH,BH,CH,DH] = dshift_nare(A,B,C,D,v,w,eta1,eta2)
% [AH,BH,CH,DH]=DSHIFT_NARE(A,B,C,D,v,w,eta1,eta2) applies the double shift 
% technique to the NARE C + XA + DX - XBX = 0
% A, B, C, D: matrix coefficients of the NARE such that  H = [A -B; -C -D]
%             has a double eigenvalue equal to zero
% v: a column vector in the right kernel of H
% w: a column vector such that w' is in the left kernel of H
% eta1, eta2: two scalars 
% AH, BH, CH, DH: matrix coefficients of the new NARE such that  
% [AH -BH; -CH -DH]*v = eta1*v, w'*[AH -BH; -CH -DH] = eta2*w'
H = [A,-B;-C,-D];
n = size(A,1);
m = size(D,1);
p = v/norm(v)^2;
q = w/norm(w)^2;
HH = H + eta1*v*p' + eta2*q*w';
AH = HH(1:n,1:n);
BH = -HH(1:n,n+1:n+m);
CH = -HH(n+1:n+m,1:n);
DH = -HH(n+1:n+m,n+1:n+m);