@Hendrik: Thank you for your contribution and the confirmation.

Closing report as fixed.

With those changes applied, I see the following for the original test case on Windows:

octave:1> glpk_bug
ctype = SSSSSUSUUUUUSUUUUUSUUUUUSUUUUU
ans =  NA
ans = 0.4099
status = 10
octave:2>

Is that the expected result?

I just realized that this patch has not been migrated.

Anybody from the core developer team with the appropriate rights  willing to review and migrate it ?

Have a good vacation.  I'll get around to updating the patch and pushing it in the next week.

@Hendrik: I agree that pivot tolerance should be updated to 1e-9.  I also think your choice of messaging level 3 is fine.  However, I don't think the presolver should be defaulted to OFF.  For the GLPK library, the documentation says that the presolver is disabled by default.  However, for the user space program which interacts with libglpk, that is glpsol.c, the default seems to be ON.  Given that Octave is attempting to act like a user space interface to libglpk, I think it should follow glpsol.c and keep the presolver enabled.  Could you check that I'm right on that?

The code I see is

      xprintf("   --presol          use presolver (default; assumes --s"
         "cale and --adv)\n");

and

      /* perform initialization */
      csa->prob = glp_create_prob();
      glp_get_bfcp(csa->prob, &csa->bfcp);
      glp_init_smcp(&csa->smcp);
      csa->smcp.presolve = GLP_ON;

Well, at least glpk (https://www.gnu.org/software/glpk/) has issues with the dodgy constraint matrix A. Below is the corresponding glpk CPLEX LP Format (plan.lp). Octave uses glpk, so calling directly glpk utilities using the CPLEX LP format avoids any potential issues coming from Octave interacting with glpk.

Minimize
value :
+1 x(2) +1 x(3) +1 x(4) +1 x(5) +1 x(6) +1 x(8) +1 x(9) +1 x(10) +1 x(11) +1 x(12) +1 x(14) +1 x(15) +1 x(16) +1 x(17) +1 x(18) +1 x(20) +1 x(21) +1 x(22) +1 x(23) +1x(24)

Subject To
+1 x(1)  +1 x(2)  +1 x(3)  +1 x(4)  +1 x(5)   = +1
+1 x(6)  +1 x(7)  +1 x(8)  +1 x(9)  +1 x(10)   = +1
+1 x(11)  +1 x(12)  +1 x(13)  +1 x(14)  +1 x(15)   = +1
+1 x(16)  +1 x(17)  +1 x(18)  +1 x(19)  +1 x(20)   = +1
+1 x(21)  +1 x(22)  +1 x(23)  +1 x(24)  +1 x(25)   = +1
+9.9999999999999967e-01 x(1)  -1.3877787807814457e-17 x(6)  -7.7023973792699711e-18 x(16)  +8.6248950136182310e-18 x(21)   <= +1.9322110003657564e+00
+9.9999999999999967e-01 x(2)  -1.3877787807814457e-17 x(7)  -7.7023973792699711e-18 x(17)  +8.6248950136182310e-18 x(22)   = +6.7788999634242716e-02
+9.9999999999999967e-01 x(3)  -1.3877787807814457e-17 x(8)  -7.7023973792699711e-18 x(18)  +8.6248950136182310e-18 x(23)   <= +6.9388939039072284e-18
+9.9999999999999967e-01 x(4)  -1.3877787807814457e-17 x(9)  -7.7023973792699711e-18 x(19)  +8.6248950136182310e-18 x(24)   <= -1.1553596068904957e-17
+9.9999999999999967e-01 x(5)  -1.3877787807814457e-17 x(10)  -7.7023973792699711e-18 x(20)  +8.6248950136182310e-18 x(25)   <= +6.4651290213686789e-18
+8.6442200073151410e-01 x(1)  +1.3557799926848546e-01 x(6)  +1.3877787807814457e-17 x(11)  -7.7037197775489434e-34 x(16)  -1.2944426998117335e-17 x(21)   <= +1.6851405085210678e+00
+8.6442200073151410e-01 x(2)  +1.3557799926848546e-01 x(7)  +1.3877787807814457e-17 x(12)  -7.7037197775489434e-34 x(17)  -1.2944426998117335e-17 x(22)   <= +2.8596532282743758e-01
+8.6442200073151410e-01 x(3)  +1.3557799926848546e-01 x(8)  +1.3877787807814457e-17 x(13)  -7.7037197775489434e-34 x(18)  -1.2944426998117335e-17 x(23)   = +2.8894168651493857e-02
+8.6442200073151410e-01 x(4)  +1.3557799926848546e-01 x(9)  +1.3877787807814457e-17 x(14)  -7.7037197775489434e-34 x(19)  -1.2944426998117335e-17 x(24)   <= -1.9255993448174939e-18
+8.6442200073151410e-01 x(5)  +1.3557799926848546e-01 x(10)  +1.3877787807814457e-17 x(15)  -7.7037197775489434e-34 x(20)  -1.2944426998117335e-17 x(25)   <= -2.1642235619263752e-17
+5.5403002969518678e-01 x(1)  +3.3039329569883769e-01 x(6)  +1.1557667460597534e-01 x(11)  -1.7527275501969224e-17 x(21)   <= +1.1079929963124644e+00
+5.5403002969518678e-01 x(2)  +3.3039329569883769e-01 x(7)  +1.1557667460597534e-01 x(12)  -1.7527275501969224e-17 x(22)   <= +6.0045746952612600e-01
+5.5403002969518678e-01 x(3)  +3.3039329569883769e-01 x(8)  +1.1557667460597534e-01 x(13)  -1.7527275501969224e-17 x(23)   <= +2.2323808686430061e-01
+5.5403002969518678e-01 x(4)  +3.3039329569883769e-01 x(9)  +1.1557667460597534e-01 x(14)  -1.7527275501969224e-17 x(24)   = +6.8311447297108668e-02
+5.5403002969518678e-01 x(5)  +3.3039329569883769e-01 x(10)  +1.1557667460597534e-01 x(15)  -1.7527275501969224e-17 x(25)   <= -2.9527020002483170e-17
+2.4336980634722233e-01 x(1)  +2.8389210464299264e-01 x(6)  +1.9949229982135033e-01 x(11)  +2.7324578918843467e-01 x(16)   <= +5.0356221694463010e-01
+2.4336980634722233e-01 x(2)  +2.8389210464299264e-01 x(7)  +1.9949229982135033e-01 x(12)  +2.7324578918843467e-01 x(17)   <= +5.0843648088919835e-01
+2.4336980634722233e-01 x(3)  +2.8389210464299264e-01 x(8)  +1.9949229982135033e-01 x(13)  +2.7324578918843467e-01 x(18)   <= +3.2813261838351937e-01
+2.4336980634722233e-01 x(4)  +2.8389210464299264e-01 x(9)  +1.9949229982135033e-01 x(14)  +2.7324578918843467e-01 x(19)   <= +4.0986868378265201e-01
+2.4336980634722233e-01 x(5)  +2.8389210464299264e-01 x(10)  +1.9949229982135033e-01 x(15)  +2.7324578918843467e-01 x(20)   = +2.5000000000000000e-01
+1 x(21)   <= +0
+1 x(22)   <= +0
+1 x(23)   <= +0
+1 x(24)   <= +0
+1 x(25)   <= +2

Bounds

integer

End

If I try to solve it via glpsol :

 glpsol  --lp --nopresol plan.lp  -o a.out ; cat a.out

...
Karush-Kuhn-Tucker optimality conditions:

KKT.PE: max.abs.err = 6.05e-02 on row 23
        max.rel.err = 5.71e-02 on row 23
        PRIMAL SOLUTION IS WRONG

KKT.PB: max.abs.err = 1.20e-17 on row 20
        max.rel.err = 1.20e-17 on row 20
        High quality

KKT.DE: max.abs.err = 1.11e-16 on column 1
        max.rel.err = 3.70e-17 on column 1
        High quality

KKT.DB: max.abs.err = 1.04e+00 on column 8
        max.rel.err = 1.04e+00 on column 8
        DUAL SOLUTION IS INFEASIBLE

and using xcheck

glpsol  --lp --nopresol --xcheck plan.lp

...
PROBLEM HAS NO FEASIBLE SOLUTION

and using presol

 glpsol  --lp --presol  plan.lp
...
PROBLEM HAS NO PRIMAL FEASIBLE SOLUTION

So as stated a few times, funny things can happen in numerical analysis when values go below eps, so it is good practice to check input and results and potentially cleans input data if results from a complex numerical algorithm look dodgy.
 
As already mentioned I am not an expert in LP nor in glpk, so feel free to report a bug in glpk (https://www.gnu.org/software/glpk/#bug).

I had another deeper look at the octave code and the glpk documentation and its sample code:

- Scaling the LP problem is a widely used approach to improve the numerical behavior of the dual simplex algorithm. So I think as well that msglev level 2 is appropriate to show when the scaling is requested.


- As seen below the constraints matrix with a number of values smaller than eps() is the root issue.
 
The (octave default) presolving converts the LP problem to an equivalent LP problem, which may be easier for solving with the simplex method than the original one. This is attained mainly due to reducing the problem size and improving its numeric properties (for example, by removing some inactive constraints or by fixing some non-basic variables).
Or for contraints matrix A with such dodgy values (< eps) presolving converts obviously the original (unsolvable) constraint matrix to a solvable albeit incorrect "equivalent" one (as already mentioned funny things can happen in numerical analysis for very small numbers (<eps)...). Therefore the return code giving back by glp_simplex did not show any issues.

Or Glpk offers sorting ("glp_sort_matrix") which should never hurt in a normal case, but which actually helps presolving to detect such ill formed contraints better. So it is done by default in their sample code before doing the solving (e.g. see examples/glpsol.c).

I tried this out calling glp_sort_matrix after loading it:

diff -r 583637c89c11 libinterp/dldfcn/__glpk__.cc
--- a/libinterp/dldfcn/__glpk__.cc
+++ b/libinterp/dldfcn/__glpk__.cc
@@ -172,6 +172,7 @@
     }

   glp_load_matrix (lp, nz, rn, cn, a);
+  glp_sort_matrix(lp);

   if (save_pb)
     {

It works and now the error is detected when presolving is active:
Return code status = 10 from glp_simplex, so "No primal feasible solution."

So I suggest doing this change and do always sorting (see above patch), as I cannot hurt, as it follows the glpk use case example and as it should help.

Some further suggested changes:
===============================
a) gplk does NOT activate presolv by default (int presolve (default: GLP OFF), so I suggest we do the same : presol (default: 0)

b) gplk has changed the default tol_piv value to 1E-9.  Currently octave still used the (former) default value 1e-10, so I suggest we change tolpiv to (default: 1e-9) as well.

On excessive messaging, I think the problem is with glpk.  Here is the Octave code for interfacing with GLPK.

  // scale the problem data
  if (! par.presol || lpsolver != 1)
    glp_scale_prob (lp, scale);

  // build advanced initial basis (if required)
  if (lpsolver == 1 && ! par.presol)
    glp_adv_basis (lp, 0);

  // For MIP problems without a presolver, a first pass with glp_simplex
  // is required
  if ((! isMIP && lpsolver == 1)
      || (isMIP && ! par.presol))
    {
      glp_smcp smcp;
      glp_init_smcp (&smcp);
      smcp.msg_lev = par.msglev;
      smcp.meth = par.dual;
      smcp.pricing = par.price;
      smcp.r_test = par.rtest;
      smcp.tol_bnd = par.tolbnd;
      smcp.tol_dj = par.toldj;
      smcp.tol_piv = par.tolpiv;
      smcp.obj_ll = par.objll;
      smcp.obj_ul = par.objul;
      smcp.it_lim = par.itlim;
      smcp.tm_lim = par.tmlim;
      smcp.out_frq = par.outfrq;
      smcp.out_dly = par.outdly;
      smcp.presolve = par.presol;
      errnum = glp_simplex (lp, &smcp);
    }

  if (isMIP)
    {
      glp_iocp iocp;
      glp_init_iocp (&iocp);
      iocp.msg_lev = par.msglev;
      iocp.br_tech = par.branch;
      iocp.bt_tech = par.btrack;
      iocp.tol_int = par.tolint;
      iocp.tol_obj = par.tolobj;
      iocp.tm_lim = par.tmlim;
      iocp.out_frq = par.outfrq;
      iocp.out_dly = par.outdly;
      iocp.presolve = par.presol;
      errnum = glp_intopt (lp, &iocp);
    }

  if (! isMIP && lpsolver == 2)
    {
      glp_iptcp iptcp;
      glp_init_iptcp (&iptcp);
      iptcp.msg_lev = par.msglev;
      errnum = glp_interior (lp, &iptcp);
    }

When lpsolver is 2, the last if statement is used and the msg_lev parameter is being passed to GLPK.  If I set the message level to 3 for maximum verbosity I get

Scaling...
 A: min|aij| =  7.704e-34  max|aij| =  1.000e+00  ratio =  1.298e+33
EQ: min|aij| =  8.912e-34  max|aij| =  1.000e+00  ratio =  1.122e+33
Original LP has 30 row(s), 25 column(s), and 115 non-zero(s)
Working LP has 30 row(s), 46 column(s), and 136 non-zero(s)
Matrix A has 136 non-zeros
Matrix S = A*A' has 165 non-zeros (upper triangle)
Approximate minimum degree ordering (AMD)...
Computing Cholesky factorization S = L*L'...
Matrix L has 185 non-zeros
Guessing initial point...
Optimization begins...
  0: obj =   1.392436989e+01; rpi =  1.5e+00; rdi =  1.9e+00; gap =  7.8e-01
  1: obj =   5.208283688e+00; rpi =  4.2e-01; rdi =  1.9e-01; gap =  1.3e+00
  2: obj =   2.433432928e+00; rpi =  8.8e-02; rdi =  7.1e-02; gap =  7.2e-01
  3: obj =   1.953120363e+00; rpi =  4.4e-02; rdi =  2.8e-02; gap =  3.6e-01
  4: obj =   2.053417812e+00; rpi =  2.5e-02; rdi =  2.1e-02; gap =  2.8e-01
  5: obj =   1.829223722e+00; rpi =  3.8e-03; rdi =  3.3e-03; gap =  4.7e-02
  6: obj =   1.791162054e+00; rpi =  3.9e-04; rdi =  3.3e-04; gap =  4.8e-03
  7: obj =   1.787311256e+00; rpi =  3.9e-05; rdi =  3.3e-05; gap =  4.8e-04
  8: obj =   1.786926127e+00; rpi =  3.9e-06; rdi =  3.3e-06; gap =  4.8e-05
  9: obj =   1.786887614e+00; rpi =  3.9e-07; rdi =  3.3e-07; gap =  4.8e-06
 10: obj =   1.786883762e+00; rpi =  3.9e-08; rdi =  3.3e-08; gap =  4.8e-07
 11: obj =   1.786883377e+00; rpi =  3.9e-09; rdi =  3.3e-09; gap =  4.8e-08
 12: obj =   1.786883339e+00; rpi =  3.9e-10; rdi =  3.3e-10; gap =  4.8e-09
OPTIMAL SOLUTION FOUND

If I switch to minimum verbosity (0) then I get

Scaling...
 A: min|aij| =  7.704e-34  max|aij| =  1.000e+00  ratio =  1.298e+33
EQ: min|aij| =  8.912e-34  max|aij| =  1.000e+00  ratio =  1.122e+33

The scaling messages appear on the terminal after the call to glpk_scale_prob() which is done before any of the solvers.

  // scale the problem data
  if (! par.presol || lpsolver != 1)
    glp_scale_prob (lp, scale);

I looked at the GLPK API and the glpk_scale_prob() doesn't take a structure argument with options to control it.  Maybe there is a way to set the global verbosity for glpk, but I don't know how to do it.  I think we might just have to live with it, unless you want to file an issue report with the GLPK project itself.

Also, this doesn't seem great, but the documentation on the return code from glpk says

Returns
0

The LP problem instance has been successfully solved. (This code does not nec-
essarily mean that the solver has found optimal solution. It only means that the
solution process was successful.)

I don't know how to interpret this.  Does this mean that GLPK has successfully found a solution which satisfies all constraints, but may not be the absolute, global mininum?  Or does it mean that the process simply completed without one of the specific failure modes such as exceeding maximum iterations?  Not great.

In any case, maybe you understand the GLPK algorithms better than I (I'm only an Octave coder).  Should we be pre-scaling problems that call glp_interior?

The previous comment was very useful for me. A bit of data cleansing seems to solve my problems.

This is surprising to me. Even without data cleansing the problem has a solution, as far as I can see, or if it hadn't I would have expected an error message (10. no primal feasible solution). I still don't fully understand how small rounding errors can have the effect they have. (Apparently I am even less of an expert on linear programming.) But I am happy that I can now use glpk.

A small side comment: if param.lpsolver = 2, then setting param.msglev = 1 does not have the desired effect (error messages only). If one uses glpk in a loop (as I do) this means a lot of annoying messages flooding the screen. But this is a minor issue and I can anyway now work with the default value param.lpsolver = 1.

I'm guessing that this is an upstream problem with glpk.  It's always possible that Octave is somehow not using the API correctly.  One way to test this would be to code this example in C or C++ and directly compile the example to link with libglpk and see if it gets the same results.

Just to check some other possibilities, I changed ctype(24) to 'S' and it does arrive at a solution where

(A*xmin)(24) == b(24)

This might be close to good enough since using 'U' would use '<=' constraint.  I also tried using 'F' constraint (free parameter, ignored), but the overall result of the minimization was worse.

c*xmin

For completeness, this is the result I get.
Since 0.4658 > 0.4099, the constraint is not satisfied.

octave:24> glpk_bug
ctype = SSSSSUSUUUUUSUUUUUSUUUUUSUUUUU
ans = 0.4658
ans = 0.4099
status = 0

For a moment I thought that I had solved my problems
by adding the line param.lpsolver = 2
which switches to the interior point method
but I have now found out this also sometimes
produces incorrect results.

So far I have only checked this on my own laptop. I am in the process of trying to check if installations of Octave on other hardware produce the same incorrect output.

Please paste the result you get.