Allow QPTDualize with only lin. inequality

[jkruzik]

+ fix unused variable warnings in nest

[mpecha]

Looks great. I've added just a few comments.

[haplav]

@jkruzik Looking at the output of the new test, even with very high rtol

mpirun -n 3 ./ex3 -n 100 -qps_view_convergence -qp_chain_view_kkt -qps_rtol 1e-14

I get quite a high residual for QP #0

r = ||A*x - b + B'*lambda|| = 9.44e-02    rO/||b|| = 3.12e+00

I'm quite afraid whether everything is mathematically correct in this case.
There might be a problem in the transform or reconstruction.

I'm investigating further.

[jkruzik]

With n=100, atol=1e-14, the last component of lambda is -0.0944291, when multiplied by Kplus this component is 3.27418e-14 which yields the correct solution (zero Dirichlet in the last component. However, Ax -b + B'lambda is almost zero everywhere except for the last component.

The examples solves the same problem as ex1. The difference between solutions is 1e-12 in the Euclidean norm.

[haplav]

With n=100, atol=1e-14, the last component of lambda is -0.0944291, when multiplied by Kplus this component is 3.27418e-14 which yields the correct solution (zero Dirichlet in the last component. However, Ax -b + B'lambda is almost zero everywhere except for the last component.

The examples solves the same problem as ex1. The difference between solutions is 1e-12 in the Euclidean norm.

@jkruzik Should I interpret this that you deem this behavior correct?

[jkruzik]

With n=100, atol=1e-14, the last component of lambda is -0.0944291, when multiplied by Kplus this component is 3.27418e-14 which yields the correct solution (zero Dirichlet in the last component. However, Ax -b + B'lambda is almost zero everywhere except for the last component.
The examples solves the same problem as ex1. The difference between solutions is 1e-12 in the Euclidean norm.

@jkruzik Should I interpret this that you deem this behavior correct?

I would say that the primal solution is close to correct.

I thought that the problem is with either the algorithm or Hessian/Hessian inversion. Turns out the latter is the case as there is an issue with MUMPS:

|| x|| = 1.51709423e+01    max( x) = 1.17e-11 =  x(99)    min( x) = -2.17e+00 =  x(41)    4578a3b0
r = ||A*x - b + B'*lambda|| = 9.44e-02    rO/||b|| = 3.12e+00
r = ||max(BI*x-cI,0)||   = 2.58e-10    r/||b|| = 8.52e-09
r = ||min(lambda_I,0)||  = 0.00e+00    r/||b|| = 0.00e+00
r = |lambda_I'*(BI*x-cI)|= 1.02e-12    r/||b|| = 3.35e-11

SuperLU_DIST (-dual_mat_inv_pc_factor_mat_solver_type superlu_dist):

|| x|| = 1.51709423e+01    max( x) = 0.00e+00 =  x(0)    min( x) = -2.17e+00 =  x(41)    75f3f020
r = ||A*x - b + B'*lambda|| = 1.61e-15    rO/||b|| = 5.32e-14
r = ||max(BI*x-cI,0)||   = 2.50e-10    r/||b|| = 8.24e-09
r = ||min(lambda_I,0)||  = 0.00e+00    r/||b|| = 0.00e+00
r = |lambda_I'*(BI*x-cI)|= 5.15e-13    r/||b|| = 1.70e-11

CG (-dual_mat_inv_ksp_type cg -dual_mat_inv_pc_type none):

|| x|| = 1.51709423e+01    max( x) = 0.00e+00 =  x(0)    min( x) = -2.17e+00 =  x(41)    1c8514c0
r = ||A*x - b + B'*lambda|| = 1.36e-14    rO/||b|| = 4.47e-13
r = ||max(BI*x-cI,0)||   = 2.50e-10    r/||b|| = 8.26e-09
r = ||min(lambda_I,0)||  = 0.00e+00    r/||b|| = 0.00e+00
r = |lambda_I'*(BI*x-cI)|= 5.32e-13    r/||b|| = 1.76e-11

Stopping criteria is on atol=1e-8. Note that the norm of the solution is the same in all cases. Maybe MUMPS has some problem with the enforcement of the Dirichlet BCs in Hessian (it seems there is the largest difference between the solutions)?

[haplav]

OK. Since something similar happens for me in ex4 in #20, I think it's probably not related to this pull request. So I think we can merge this and create a separate issue - it definitely requires further attention.

[jkruzik]

Works OK with MUMPS LU (-dual_mat_inv_pc_type lu):

|| x|| = 1.51709423e+01    max( x) = 0.00e+00 =  x(0)    min( x) = -2.17e+00 =  x(41)
r = ||A*x - b + B'*lambda|| = 2.19e-15    rO/||b|| = 7.24e-14
r = ||max(BI*x-cI,0)||   = 2.50e-10    r/||b|| = 8.24e-09
r = ||min(lambda_I,0)||  = 0.00e+00    r/||b|| = 0.00e+00
r = |lambda_I'*(BI*x-cI)|= 5.17e-13    r/||b|| = 1.71e-11

I will switch to LU, add require: mumps and merge this.

[haplav]

Works OK with MUMPS LU (-dual_mat_inv_pc_type lu):

|| x|| = 1.51709423e+01    max( x) = 0.00e+00 =  x(0)    min( x) = -2.17e+00 =  x(41)
r = ||A*x - b + B'*lambda|| = 2.19e-15    rO/||b|| = 7.24e-14
r = ||max(BI*x-cI,0)||   = 2.50e-10    r/||b|| = 8.24e-09
r = ||min(lambda_I,0)||  = 0.00e+00    r/||b|| = 0.00e+00
r = |lambda_I'*(BI*x-cI)|= 5.17e-13    r/||b|| = 1.71e-11

I will switch to LU, add require: mumps and merge this.

Nice! So a problem with Cholesky - could be a symmetry broken or something?

[jkruzik]

OK, I am stupid :(

The assembly of the tridiag matrix does not enforce Dirichlet correctly - which breaks symmetry.

[haplav]

Nice, so it seems #22 gets solved before I really got to that ;-)

[haplav]

I'm just a bit surprised MUMPS proceeds without any error/warning...

[jkruzik]

Fixed. I am going to fix ex1 and ex2 as well in a separate PR.

[haplav]

Sounds good.