program main
use qrgivens
implicit none
real*4,allocatable::A(:,:),B(:,:),Q(:,:),ai(:,:),R(:,:),v(:),x(:)
integer::n=4,m=3,i,j

print*,"QR-DECOMPOSITION BY GIVENS ROTATIONS:"
allocate(A(n,m))
call RANDOM_NUMBER(A)
print*,"random tall matrix A:"; call printm(A);
allocate(B(n,m))
B(:,:)=A(:,:)

call qrdec(A)
allocate(Q(n,m))
call qrgetq(A,Q)
allocate(R(m,m))
call qrgetr(A,R)
print*,"check: Q*R: ==A?"; call printm(MATMUL(Q,R))
print*,"check: Q^T*Q: ==1?"; call printm(MATMUL(transpose(Q),Q))
print*,"check: R: upper triangular?"; call printm(R)
print*,"check: Q^T*A: ==R?"; call printm(MATMUL(transpose(Q),B))

print*,"LINEAR SYSTEM OF EQUATIONS:"
deallocate(A,B,Q,R)
allocate(A(m,m),B(m,m),Q(m,m),R(m,m),v(m),x(m))
call RANDOM_NUMBER(A)
print*,"random square matrix A:"; call printm(A);
B(:,:)=A(:,:)
call RANDOM_NUMBER(v)
print*,"random right-hand-side v:"; call printv(v)

call qrdec(A)
x=qrsolve(A,v)
print*,"solution x to A*x=v:"; call printv(x)
print*,"check: A*x: ==v?";call printv(matmul(B,x))

print*,"MATRIX INVERSE:"
ai=qrinv(A)
print*,"inverse matrix A^-1:";call printm(ai)
print*,"check: A^-1*A: ==1?";call printm(matmul(ai,b))
print*,"check: A*A^-1: ==1?";call printm(matmul(b,ai))

end program main