program LinearEquations
! solving the matrix equation A*x=b using LAPACK
  implicit none
  double precision :: M2mat(2,2), T2vec(2), a2vec(2), tmpp(2,2), tmpp2(2)
  integer :: info, pivot2(2)
  character :: msg*100

  M2mat = reshape( (/3.0, 2.0, 2.0, 6.0/),(/2,2/))
  T2vec = (/5.0, 8.0/)
  tmpp = M2mat
  tmpp2 = T2vec
  call dgesv (2, 1, tmpp, 2, pivot2, tmpp2, 2, info)
  if (info /= 0) then
    write (*, '("Error in call to POSV, info=",i4)') info
  end if
  a2vec = tmpp2

  write(*,*) "result of test 2x2 matrix solve is a2vec = ",a2vec(1),a2vec(2)

  tmpp2 = matmul(M2mat,a2vec)

  write(*,*) "check: M2mat multiplying solution yields ",tmpp2(1),tmpp2(2)

end program