Hello,
I have implemented a mex c++ interface to a nonlinear programming solver and until now it worked great. I am trying to cover special cases and when trying to deal with vanishing Hessians it seems to have a werid issue and I cannot detemine where it comes from. I am getting confusing information about where the problem occurs but it seems somehow related to sparse matrices. So I attach my class members that interface to the Jacobian and Hessian.
When I attach the debugger (XCode) for some reason it crashes matlab, without the debugger an exception is thrown and a regular error tells me Can't index into an empty array – obviously not having the debugger working is not ideal. According to my printout the last call I get on the output comes from the Jacobian evaluation, but after the succesfull passing of the numerical values. For the one with the debugger attached, the crashdump points to the Hessian evaluation.
Is there a way to determine where the exception was thrown in the c++ stack? Currently I only have this info:
MException with properties: identifier: 'MATLAB:mex:CppMexException' message: 'Can't index into an empty array' cause: {} stack: [2×1 struct] Correction: []
and the stack is on the Matlab side.
Thank you,
Manuel
bool eval_jac_g(Index n, const Number* x, bool new_x, Index m, Index nele_jac, Index* iRow, Index *jCol, Number* values) {#if defined(DBG) std::cout << "Called eval_jac_g" << std::endl;#endif if (nullptr == values){ std::vector<Array> args(0); std::vector<Array> jacStrOut = feval(funcs[0]["jacobianstructure"], 1, args); SparseArray<double> JacobianStr(std::move(jacStrOut[0]));#if defined(DBG) std::cout << " --- Jacobian Structure ---" << std::endl;#endif auto i = 0; for (TypedIterator<double> it = JacobianStr.begin(); it != JacobianStr.end(); it++){ iRow[i] = JacobianStr.getIndex(it).first; jCol[i] = JacobianStr.getIndex(it).second;#if defined(DBG) std::cout << "iRow[" << i << "] = " << iRow[i] << " jCol[" << i << "] = " << jCol[i] << std::endl;#endif i++; } } else { if (new_x) updateX(x);#if defined(DBG) std::cout << " --- Jacobian values ---" << std::endl;#endif SparseArray<double> Jacobian = fevalWithX(funcs[0]["jacobian"]); auto i = 0; for (auto& elem : Jacobian){#if defined(DBG) std::cout << "values[" << i << "] = " << elem << std::endl;#endif values[i] = elem; i++; } } return true; };
bool eval_h(Index n, const Number* x, bool new_x, Number obj_factor, Index m, const Number* lambda, bool new_lambda, Index nele_hess, Index* iRow, Index* jCol, Number* values) {#if defined(DBG) std::cout << "Called eval_h" << std::endl;#endif assert(!returnHessian); if (nullptr == values){ Array hessianStructure = funcs[0]["hessianstructure"]; std::vector<Array> args(0); std::vector<Array> hesStrOut = feval(hessianStructure, 1, args); SparseArray<double> HessianStr(std::move(hesStrOut[0])); auto i = 0; for (TypedIterator<double> it = HessianStr.begin(); it != HessianStr.end(); it++){ iRow[i] = HessianStr.getIndex(it).first; jCol[i] = HessianStr.getIndex(it).second; i++; } } else { if (new_x) updateX(x); std::vector<Array> hessianArgs = { args[0], factory.createScalar(obj_factor), factory.createArray<double>({static_cast<size_t>(m),1}, lambda, lambda + m) }; std::vector<Array> retVals = feval(funcs[0]["hessian"], 1, hessianArgs); SparseArray<double> Hessian = std::move(retVals[0]); auto i = 0; for (auto& elem : Hessian){ values[i] = elem; i++; } } return true; };
Best Answer