program main use matrix implicit none integer, parameter :: n = 3, m = 3 real :: A(1:n,1:m), Q(1:n,1:m), R(1:m,1:m) real :: inv(n,n) real :: b(n) ! Define the matrix: A = 1.0 A(2,2) = 3.0 A(3,3) = 5.0 b = 1.0 b(2) = 4.0 call MatOut('A =', A) call VecOut('b =', b) print *, '---------------------' ! Do QR-decomposition: call qrdec(A, Q, R) ! Check: call MatOut('Q =', Q) call MatOut('R =', R) call MatOut('QR =', matmul(Q,R)) call MatOut('Q^T*Q =', matmul(transpose(Q),Q)) ! Backsubstitution: call VecOut('Ax=b => x =', backsub(R,b)) ! Determinant: print *, '---------------------' print '("abs(det(A)) = ",f8.3)', determinant(R) ! Invers: print *, '---------------------' inv = inverse(Q, R) call MatOut('A^-1 =', inv) call MatOut('A^-1*A =', matmul(A,inv)) end program main