Problem implementing a QR factorization I'm trying to write a Fortran subroutine to compute a QR factorization using the Householder method. To test my routine, I compute the factorization of the following matrix:
$$
A =
 \begin{pmatrix}
  12 & -51 & 4 \\
  6 & 167 & -68  \\
  -4 & 24 & -41 
 \end{pmatrix},
$$
which, if done correctly, will reduce to the following upper triangular matrix:
$$
R =
 \begin{pmatrix}
  14 & 21 & -14 \\
  0 & 175 & -70  \\
  0 & 0 & 35 
 \end{pmatrix}.
$$
However, the matrix I actually get is:
$$
R =
 \begin{pmatrix}
  -14 & -21 & 14 \\
  0 & -175 & 70  \\
  -0 & 0 & 35 
 \end{pmatrix},
$$
which looks almost correct, except for some strange sign changes. I've been staring at my subroutine all day trying to see where these sign changes are being introduced, but I can't identify the problem. 
My algorithm is as follows:
$
for \:\: k \:=\: 1\: to\: n
$
$
\qquad x(k:m) = A(k:m,k)
$ 
$
\qquad v(k:m) = \mathtt{sign}(x(k))||x(k:m)||_{2}e1 + x(k:m)
$
$
\qquad v(k:m) = v(k:m)/||v(k:m)||_{2}
$
$
\qquad A(k:m,k:n) = A(k:m,k:n) - 2vv^{\top}A(k:m,k:n)
$
To calculate the factor $2vv^{\top}A(k:m,k:n),$ I made another subroutine called outer_product to compute the outer product of $v$ with itself, i.e. $vv^{\top}$, and then matrix multiply the result into my submatrix $A(k:m,k:n)$. However, I'm not sure if this is legitimate - I suspect herein lies the problem.
I would really appreciate it if someone could glance at my code to see if there is any obvious reason for the incorrect sign changes: 
integer, parameter :: dp = selected_real_kind(15)

integer, intent(in) :: m, n
real(dp), dimension(m,n), intent(inout) :: A
real(dp), dimension(m,m), intent(out) :: Q

integer :: k
real(dp) :: two_norm
real(dp), dimension(m) :: x, e1
real(dp), dimension(m,n) :: v
real(dp), dimension(m,m) :: outprod_vv

v = 0.0_dp

do k=1,m
    Q(k,k) = 1.0_dp
end do

!Householder triangularization
do k=1,n

    e1(k) = 1.0_dp

    x(k:m) = A(k:m,k)
    v(k:m,k) = sign( sqrt(dot_product(x(k:m),x(k:m))), x(k) )* &
        e1(k:m) + x(k:m)

    v(k:m,k) = v(k:m,k)/(sqrt(dot_product(v(k:m,k),v(k:m,k))))
    call outer_product(v(k:m,k), m-k+1, outprod_vv(k:m,k:m))

    A(k:m,k:n) = A(k:m,k:n) - &
        2.0_dp*matmul(outprod_vv(k:m,k:m), A(k:m,k:n)) 

    !Form Q implicitly    
    Q(k:m,k:m) = Q(k:m,k:m) - 2.0_dp* &
        matmul(outprod_vv(k:m,k:m), Q(k:m,k:m))

end do

Q = transpose(Q)

 A: A Householder reflection $Q:=I-2vv^T$ which transforms a vector $x$ to a distinct vector $y$ (such that $\|x\|_2=\|y\|_2$) can be constructed by setting $v:=(x-y)/\|x-y\|_2$. We have
$$
Qx=\left[1-\frac{2(x-y)^Tx}{(x-y)^T(x-y)}\right]x+\frac{2(x-y)^Tx}{(x-y)^T(x-y)}y.
$$
Using the above assumptions, it is straightforward to verify that $Qx=y$ (the coefficients by $x$ and by $y$ are, respectively, zero and one).
So, to zero out all the components of $x$, $y:=\sigma e_1$, where $|\sigma|=\|x\|_2$. Since one usually wants to avoid severe cancelation in computing $v$, $\sigma$ is chosen so that the sign of the first component of $x$ is the same as the sign of $\sigma$.
Note that, in $Qx=\sigma e_1$, $\sigma$ is negative if the first component of $x$ is positive so that explains why HHQR creates the upper triangular matrix with some negative diagonal entries.
If you would like to have a positive diagonal in $R$, consider multiplying by $-Q=2vv^T-I$ instead of by $Q$.
Since you showed the code, I would have some remarks:


*

*You do not need a special vector for $e_1$ since adding its multiple to a vector does nothing but modifies one entry.

*To compute $QB=B-2vv^TB$, it is better to implement it as $QB=B-2v(v^TB)$ instead of $QB=B-2(vv^T)B$. I believe you can see why if, e.g., $m\gg n$.

*Although in HHQR, it is better to do not store $Q$ explicitly (think again about the case $m\gg n$; not mentioning numerical reasons now), your $Q$ is certainly incorrect. You should update all the columns of $Q$ instead of $k$ to $m$.

*Since by construction of $Q$, you know what should be the content of $A(k:m,k)$ (its first component is up to sign equal to the norm of $v$ and the rest is zero). So computing it explicitly by applying $Q$ is, firstly, waste of time, and, secondly, instead of having zeros in rows $k+1,\ldots,m$, you are very likely to get some small junk due to roundoff in the strictly lower triangular part of $R$. 
For inspiration, you can have a look on my answer here; you might be intereseted there in particular in the implementation of the orthogonalization loop.
I would also suggest to have a look on Chapter 5 here.
A: A quick glance tells me your routine is fine.  But more importantly, since QR factorization is not unique, the $R$ your routine outputs is perfectly fine-- it just means that the $Q$ has absorbed a negative.  Also, for comparison, the Householder routine in Matlab gives the same $R$ as your Fortran routine.
If you implement the modified Gram-Schmidt (MGS) routine, you will obtain the $R$ you originally sought, since MGS sets the diagonal entries to positive values.
In practice, this difference in sign doesn't make a difference, since one usually cares more about the product $QR$ than about the individual factors.
