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

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

call qrdec(A,R)
q=a
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,r)
call qrsolve(A,r,v,x)
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,r)
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