wip

Author: amontoison

src/moi_nlp_model.jl

@@ -311,7 +311,7 @@ function NLPModels.jprod_nln!(
   increment!(nlp, :neval_jprod_nln)
   if nlp.quadcon.nquad > 0
     for i = 1:(nlp.quadcon.nquad)
-      # Jv[i] = (Aᵢ * x + bᵢ)ᵀ * v
+      # Jv[i] += (Aᵢ * x + bᵢ)ᵀ * v
       qcon = nlp.quadcon.constraints[i]
       Jv[i] = coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
     end
@@ -333,19 +333,17 @@ function NLPModels.jprod!(
   if nlp.meta.nlin > 0
     jprod_lin!(nlp, nlp.lincon.jacobian.rows, nlp.lincon.jacobian.cols, nlp.lincon.jacobian.vals, v, Jv)
   end
-  if nlp.meta.nnln > 0
-    if nlp.quadcon.nquad > 0
-      for i = 1:(nlp.quadcon.nquad)
-        # Jv[i] = (Aᵢ * x + bᵢ)ᵀ * v
-        qcon = nlp.quadcon.constraints[i]
-        Jv[nlp.meta.nlin + i] = coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
-      end
-    end
-    if nlp.meta.nnln > nlp.quadcon.nquad
-      ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon)
-      MOI.eval_constraint_jacobian_product(nlp.eval, view(Jv, ind_nnln), x, v)
+  if nlp.quadcon.nquad > 0
+    for i = 1:(nlp.quadcon.nquad)
+      # Jv[i] = (Aᵢ * x + bᵢ)ᵀ * v
+      qcon = nlp.quadcon.constraints[i]
+      Jv[nlp.meta.nlin + i] = coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
     end
   end
+  if nlp.meta.nnln > nlp.quadcon.nquad
+    ind_nnln = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon)
+    MOI.eval_constraint_jacobian_product(nlp.eval, view(Jv, ind_nnln), x, v)
+  end
   return Jv
 end
 
@@ -378,7 +376,7 @@ function NLPModels.jtprod_nln!(
     ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln)
     MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, view(v, ind_nnln))
   end
-  (nlp.meta.nnln == nlp.quadcon.nquad) && (Jtv .= 0.0)
+  (nlp.quadcon.nquad > 0) && (Jtv .= 0.0)
   if nlp.quadcon.nquad > 0
     for i = 1:(nlp.quadcon.nquad)
       # Jtv += v[i] * (Aᵢ * x + bᵢ)
@@ -401,10 +399,11 @@ function NLPModels.jtprod!(
     ind_nnln = (nlp.quadcon.nquad + 1):(nlp.meta.nnln)
     MOI.eval_constraint_jacobian_transpose_product(nlp.eval, Jtv, x, view(v, ind_nnln))
   end
+  (nlp.quadcon.nquad > 0 || nlp.meta.nlin > 0) && (Jtv .= 0.0)
   if nlp.meta.nlin > 0
-    jtprod_lin!(nlp, nlp.lincon.jacobian.rows, nlp.lincon.jacobian.cols, nlp.lincon.jacobian.vals, v, Jtv)
+    transpose = true
+    coo_unsym_add_mul!(transpose, nlp.lincon.jacobian.rows, nlp.lincon.jacobian.cols, nlp.lincon.jacobian.vals, v, Jtv, 1.0)
   end
-  (nlp.meta.nnln == nlp.quadcon.nquad) && (Jtv .= 0.0)
   if nlp.quadcon.nquad > 0
     for i = 1:(nlp.quadcon.nquad)
       # Jtv += v[i] * (Aᵢ * x + bᵢ)

src/utils.jl

@@ -116,10 +116,37 @@ function replace!(ex, x)
   end
 end
 
+"""
+    coo_unsym_add_mul!(transpose, rows, cols, vals, x, y, α)
+
+Performs the update `y ← y + α * op(A) * x`, where `A` is an unsymmetric matrix in COO format given by `(rows, cols, vals)`.
+If `transpose == true`, then `op(A) = Aᵀ`; otherwise, `op(A) = A`.
+"""
+function coo_unsym_add_mul!(
+  transpose::Bool,
+  rows::AbstractVector{<:Integer},
+  cols::AbstractVector{<:Integer},
+  vals::AbstractVector,
+  x::AbstractVector,
+  y::AbstractVector,
+  α::Float64,
+)
+  nnz = length(vals)
+  @inbounds for k = 1:nnz
+    i, j, c = rows[k], cols[k], vals[k]
+    if transpose
+      y[j] += α * c * x[i]
+    else
+      y[i] += α * c * x[j]
+    end
+  end
+  return y
+end
+
 """
     coo_sym_add_mul!(rows, cols, vals, x, y, α)
 
-Update of the form `y ← y + αAx` where `A` is a symmetric matrix given by `(rows, cols, vals)`.
+Perform the update `y ← y + α * A * x` where `A` is a symmetric matrix in COO format given by `(rows, cols, vals)`.
 Only one triangle of `A` should be passed.
 """
 function coo_sym_add_mul!(