N   = 10;
A   = rand(N,N);
AtA = A'*A;

cond(A)

(cond(A))^2

cond(AtA)

## Input independent vectors
v1 = [1.; 0.; 0.];
v2 = [1.; 1.; 1.];
v3 = [1.; 1.; 0.];
A = [v1 v2 v3]

# VECTOR v1.
# Normalize v1:
a1 = v1/norm(v1)   ## book calls these a_i but q_i is better

# VECTOR v2.
# Remove the span of a1 from v2:
a2 = v2 - a1*a1'*v2

# normalize a2:
a2 /= norm(a2)

# VECTOR v3.
# Remove the span of a1 & a2 from v3: 
a3 = v3 - a1*a1'*v3 - a2*a2'*v3

# normalize a3:
a3 /= norm(a3)

## Orthogonal basis matrix formed by Gram-Schmidt:
Q = [a1 a2 a3]

Q'*Q

## Upper triangle matrix (encoded in our operations above)
R = Q'*A



?(qr)

?(qrfact)



M = 5;
N = 2;
A = rand(M,N)

QR1 = qrfact(A, Val{false})  ## QR with pivot=false (for demo purposes)

Q1 = QR1[:Q]   ## "Fat" QR with null-space basis, too

R1 = QR1[:R]

## Note: Julia permits dimension mismatch here with QRCompactWYQ
Q1 * R1  

A



(Q2, R2) = qr(A)  ## Default: unpivoted "thin QR"

Q2  ## thin basis

R2

Q2*R2

A

b = ones(M)/sqrt(M); ## unit vector
x1 = QR1 \ b ## QR1 from qrfact (fat Q)

x2 = R2 \ (Q2'*b) ## QR2 from qr (thin Q, no pivoting)

## Residual:
norm(A*x1-b), norm(A*x2-b)