Fortran 95 source code to solve a linear programming by simplex method.

simplex.f95

!****************************************************
! 線形計画法（シンプレックス法）
!
! * 入力はテキストファイルをパイプ処理
!     1行目:     行数  列数  変数の数
!     2行目以降: 1行に列数分の係数 * 行数
!
!   date          name            version
!   2018.12.05    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
module const
  ! SP: 単精度(4), DP: 倍精度(8)
  integer,     parameter :: SP = kind(1.0)
  integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
end module const

module simplex
  use const
  implicit none
  private
  public :: solve

contains
  ! 線形計画法
  !
  ! :param(in)    integer(4)       n_row: 元数
  ! :param(in)    integer(4)       n_col: 元数
  ! :param(inout) real(8) a(n_row,n_col): 係数配列
  subroutine solve(n_row, n_col, a)
    implicit none
    integer(SP), intent(in)    :: n_row, n_col
    real(DP),    intent(inout) :: a(n_row, n_col)
    integer(SP) :: x, y              ! 最小値の行・列インデックス
    real(DP)    :: c_min             ! 最小値
    real(DP)    :: a_tmp(n_row - 1)  ! 作業用配列
    integer(SP) :: i, j
    real(DP)    :: d

    do
      ! 列選択
      c_min = minval(a(n_row, :))
      y = minloc(a(n_row, :), 1)
      if (c_min >= 0.0_DP) exit

      ! 行選択
      a_tmp = a(1:n_row-1,n_col) / a(1:n_row-1, y)
      x = minloc(a_tmp, 1)

      ! ピボット係数を p で除算
      a(x, :) = a(x, :) / a(x, y)

      ! ピボット列の掃き出し
      do i = 1, n_row
        if (i == x) cycle
        d = a(i, y)
        a(i, :) = a(i, :) - d * a(x, :)
      end do
    end do
  end subroutine solve
end module simplex

program linear_programming
  use const
  use simplex
  implicit none
  integer(SP) :: n_row, n_col, n_var  ! 行数、列数、変数の数
  real(DP), allocatable :: a(:, :)    ! 係数行列
  character(20) :: f                  ! 書式文字列
  integer(SP)   :: i, j, flag
  real(DP)      :: v

  ! 行数、列数、変数の数の取得
  read (*, *) n_row, n_col, n_var
  print '(A)', "Number:"
  print '("  rows      = ", I0)', n_row
  print '("  colmuns   = ", I0)', n_col
  print '("  variables = ", I0)', n_var

  ! 係数行列メモリ確保
  allocate(a(n_row, n_col))

  ! 係数行列読み込み
  do i = 1, n_row
    read (*, *) a(i, :)
  end do
  print '(A)', "Coefficients:"
  write (f, '("(", I0, "(F8.4, 2X)", ")")') n_col
  print f, a(1:n_row, :)

  ! 線形計画法で解く
  call solve(n_row, n_col, a)

  ! 結果出力
  print '(A)', "Answer:"
  do j = 1, n_var
    flag = -1
    do i = 1, n_row
      if (a(i, j) == 1.0_DP) flag = i
    end do
    if (flag == -1) then
      v = 0.0_DP
    else
      v = a(flag, n_col)
    end if
    print '("  x", I0, " = ", F8.4)', j, v
  end do
  print '("  z  = ", F8.4)', a(n_row, n_col)

  ! 係数行列メモリ解放
  deallocate(a)
end program linear_programming

You have a very restrictive license used here. All rights reserved means nobody can use it. Was that your intent?

Nowadays, the practice of writing "All rights reserved" may not make sense.
However, the author name display has the advantage of being able to clearly indicate the right holder.
So I dare to write it.

• I'm sorry for the poor English.

Can you supply an example? An input file. Better would be a hardwired example called as a subroutine.

Please refer to https://www.mk-mode.com/blog/2019/03/02/fortran95-linear-programming-by-simplex/

What would be the text file for this example?

Maximize :
15 X1 + 17 X2 + 20 X3

Conditions:
0 X1 + 1 X2 - 1 X3 <= 2
3 X1 + 3 X2 + 5 X3 <= 15
3 X1 + 2 X2 + 1 X3 <= 8

0 <= X1
0 <= X2
0 <= X3

The text file for this example:
Imagine an algorithm)

(Filename: data.txt)

4 7 3
 0  1 -1 1 0 0  2
 3  3  5 0 1 0 15
 3  2  1 0 0 1  8
15 17 20 0 0 0  0

And then, run it like this:

cat data.txt | ./simplex

Then, it should look like this:

Number:
  rows      = 4
  colmuns   = 7
  variables = 3
Coefficients:
  0.0000    1.0000   -1.0000    1.0000    0.0000    0.0000    2.0000
  3.0000    3.0000    5.0000    0.0000    1.0000    0.0000   15.0000
  3.0000    2.0000    1.0000    0.0000    0.0000    1.0000    8.0000
 15.0000   17.0000   20.0000    0.0000    0.0000    0.0000    0.0000
Answer:
  x1 =   0.0000
  x2 =   2.0000
  x3 =   8.0000
  z  =   0.0000

• Depending on the compilation environment, the output direction of "Coefficients" may be switched.

WolframAlpha.com and Mathematica (on my computer) give entirely different answers than yours.
You can run this same example in WolframAlpha by clicking on this in url:

https://www.wolframalpha.com/input?i=%5B%2F%2Fmath%3Amaximize%2F%2F%5D+%5B%2F%2Fmath%3A15x_1%2B17x_2%2B20x_3%2F%2F%5D+s.t.+%5B%2F%2Fmath%3Ax_2+-+x_3+%3C%3D+2%2C%2F%2F%5D+%5B%2F%2Fmath%3A3x_1+%2B+3x_2+%2B5x_3+%3C%3D+15%2C%2F%2F%5D+%5B%2F%2Fmath%3A3x_1+%2B+2x_2+%2B1x_3+%3C%3D+8%2C%2F%2F%5D+%5B%2F%2Fmath%3Ax_1%3E%3D0%2C+x_2%3E%3D0%2Cx_3%3E%3D0%2F%2F%5D++%5B%2F%2Fmath%3A%2F%2F%5D+%5B%2F%2Fmath%3A%2F%2F%5D+%5B%2F%2Fmath%3A%2F%2F%5D

In the line and in the code:

colmuns = 7
print '(" colmuns = ", I0)', n_col

The correct spelling is:

columns

should it be negatives:
-15 -17 -20 0 0 0 0

Your code did not agree with Wolfram Alpha results or 2 other LPs
Here's what fixed it.
! 1.. Avoid division by zero: Add a small tolerance (eps) to prevent
! divisions by zero or very small values.
! 2.. Prevent cycling (Bland's Rule): Implement Bland's Rule to avoid
! cycles by selecting the first valid pivot instead of the
! smallest negative value.
! 3.. Iteration limit: Add a maximum iteration limit to prevent
! infinite loops.

module simplex
use const
implicit none
private
public :: solve

contains

! Linear Programming
!
! :param(in) integer(4) n_row: Number of rows
! :param(in) integer(4) n_col: Number of columns
! :param(inout) real(8) a(n_row,n_col): Coefficient Array

subroutine solve(n_row, n_col, a)

implicit none
integer(SP), intent(in)    :: n_row, n_col
real(DP),    intent(inout) :: a(n_row, n_col)
integer(SP) :: x, y              ! row/column index of minimum value
real(DP)    :: a_tmp(n_row - 1)  ! working array for ratios
integer(SP) :: i, iter            ! iteration counter
real(DP)    :: d, eps             ! pivot factor and epsilon
integer(SP), parameter :: max_iter = 1000  ! maximum number of iterations

eps = 1.0E-8  ! Small tolerance to avoid division by zero or very small values
iter = 0      ! Initialize iteration counter

do while (iter < max_iter)
  iter = iter + 1

  ! Column selection (Bland's Rule: Choose the first negative entry in the last row)
  y = 0
  do i = 1, n_col-1
    if (a(n_row, i) < -eps) then
      y = i
      exit
    end if
  end do
  if (y == 0) exit  ! Optimal solution found if no negative entries in the last row

  ! Line selection (choose pivot row with the smallest positive ratio)
  do i = 1, n_row-1
    if (a(i, y) > eps) then
      a_tmp(i) = a(i, n_col) / a(i, y)
    else
      a_tmp(i) = HUGE(1.0_DP)  ! Set to large value if division is invalid
    end if
  end do
  x = minloc(a_tmp, 1)

  ! Check for degenerate case where no valid pivot row is found
  if (a(x, y) < eps) exit

  ! Divide the pivot row by the pivot element
  a(x, :) = a(x, :) / a(x, y)

  ! Pivot Column Sweep
  do i = 1, n_row
    if (i == x) cycle  ! Skip the pivot row
    d = a(i, y)
    a(i, :) = a(i, :) - d * a(x, :)  ! Adjust all other rows
  end do
end do

if (iter >= max_iter) then
  print *, "Reached maximum number of iterations (", max_iter, "). Exiting loop."
end if

end subroutine solve

end module

Thank you for pointing that out.
If I modify it as you suggested, it will become your code, so I will not modify it.
(Currently, I do not have an environment where I can use Fortran95, so I cannot try it out...)

In the future, if you see this code, please keep in mind that this point was made.