### 2.36 Solve $$Ax=b$$

Solve for x in the following system of equations

[ 1 2 3 ][x1]   [ 1 ]
[ 7 5 6 ][x2]=  [ 2 ]
[ 7 8 9 ][x3]   [ 3 ]

 Mathematica mat={{1,2,3}, {7,5,6}, {7,8,9}}; b = {1,2,3} LinearSolve[mat,b] {0,0,1/3}

 Matlab format long A=[1 2 3; 7 5 6; 7 8 9]; b=[1;2; 3]; A\b ans = 0 0 0.333333333333333

 Maple restart; A:=Matrix([[1,2,3],[7,5,6],[7,8,9]]); b:=Vector([1,2,3]); LinearAlgebra:-LinearSolve(A,b) [0, 0, 1/3]

Fortran

program t2
implicit none
integer, parameter :: N=3
real(8),   DIMENSION(N, N) :: &
A=DBLE(reshape([ 1.0, 2.0, 3.0,&
7.0, 5.0, 6.0,&
7.0, 8.0, 9.0], shape(A), order=[2,1]))

real(8), parameter, DIMENSION(N, 1) :: &
b=DBLE(reshape([1.0, 2.0, 3.0], shape(b)))

real(8), DIMENSION(N, 1) :: IPIV
integer :: INFO

CALL DGETRF( N, N, A, N, IPIV, INFO )
if (INFO .eq. 0 ) then
CALL DGETRS( 'N', N,1, A, N, IPIV, b , N, INFO)
if ( INFO .eq. 0 ) then
print *, 'solution is',b
else
print *, 'failed DGETRS, INFO=',INFO
end if
else
print *, 'failed DGETRF, INFO=',INFO
end if
end program

compile and run

$gfortran -std=f2003 -Wextra -Wall -pedantic -funroll-loops -ftree-vectorize -march=native -Wsurprising -Wconversion t2.f90 /usr/lib/liblapack.a /usr/lib/libblas.a$ ./a.exe
solution is
-5.28677630773884192E-018 -3.70074341541718826E-017  0.33333333333333337