program sgetrf_example
!***********************************************************************
!*                                                                     *
!* Kadin Tseng                                                         *
!* Scientific Computing and Visualization                              *
!* Boston University                                                   *
!* 1998 (modified 10/5/2006)                                           *
!*                                                                     *
!* This program solves                                                 *
!*                                                                     *
!*  3x - 2y + 2z = 10                                                  *
!*   x + 2y - 3z = -1                                                  *
!*  4x +  y + 2z =  3                                                  *
!*                                                                     *
!* using the LAPACK routines sgemv, sgetrf and sgetrs.                 *
!*                                                                     *
!***********************************************************************
   implicit none
   integer :: info, n=3   ! n is square matrix size
   integer, dimension(:), allocatable :: ipiv
   real,    dimension(:), allocatable :: x, b   ! solution and RHS
   real,  dimension(:,:), allocatable :: a      ! LHS matrix a

! Demonstrates LAPACK usage with Ax = b
   allocate(a(n,n), ipiv(n), x(n), b(n))

   x = (/ 2., -3., -1. /)  ! solution x
   a = reshape( (/ 3.,1.,4.,-2.,2.,1.,2.,-3.,2./), (/3,3/))

   call sgemv('N',n,n,1.0,a,n,x,1,0.0,b,1)   ! b = Ax
   write(*,*)'The RHS, b = Ax :'
   write(*,*)b
   write(*,*)'The solution, x :'
   write(*,*)x
   call sgetrf(n,n,a,n,ipiv,info)            ! factor A
   call sgetrs('N',n,1,a,n,ipiv,b,n,info)    ! solve for x

   write(*,*)'The solution from sgetrs :'
   write(*,*)b
   deallocate(a, ipiv, x, b)
end program sgetrf_example