program assgin7_1
!
use denseop
integer :: j
real(8)::a(3,3)=reshape( (/36,30,18,30,41,23,18,23,14/),(/3,3/) ),&
b(3)=(/150,181,106/),x(3)
real(8)::v(3,3),d(3),l(3,3),i(3,3),c(3,3),dd(3,3),det

call printmat(a,'A matrix')

l=fchol(a)
v=matmul(l,transpose(l))
call printmat(v,'a) After cholesky decomposition')

call eigen(a,d,v)
call printmat(v,'b) Eigenvalues')
dd=diag(d)
!dd=0
!do j=1,3
! dd(j,j)=d(j) ! dd=diag(d) does not work with r4
!enddo
call printmat(dd,'Eigenvectors')
c=matmul(v,matmul(dd,transpose(v)))
call printmat(c,'VDV')
i=matmul(v,transpose(v))
call printmat(i,'VV')
det=fdet_s(a)
print*,'Determinant',det

call solve_s(A,b,x)
print*,'d) solutions:',x
c=finverse_s(a)
x=matmul(c,b)
print*,' solutions:',x

i=matmul(a,c)

call printmat(c,'A inverse')
call printmat(i,'I=AA-1')

end program assgin7_1