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