[numerics] gfc3d_ipm : test if the Jacobian matrix contains NaN
- only the 3x3_QpH case
- to avoid MA57 to run into troubles with incinsistent matrix
- shoul be also extended to other cases for JR
@@ -1271,7 +1271,7 @@ void gfc3d_IPM(GlobalFrictionContactProblem* restrict problem, double* restrict
/* r_rhs = [ ] */
/* [ -Qp*(H*v + w) ] nd */
-
+ int jacobian_is_nan = 0;
switch ( options->iparam[SICONOS_FRICTION_3D_IPM_IPARAM_LS_FORM] )
{
case SICONOS_FRICTION_3D_IPM_IPARAM_LS_3X3_NOSCAL:
@@ -1371,6 +1371,7 @@ void gfc3d_IPM(GlobalFrictionContactProblem* restrict problem, double* restrict
NM_insert(J, minusQpHt, 0, m + nd);
NM_insert(J, eye_nd, m, m + nd);
+
NM_free(minusQpH);
NM_free(minusQpHt);
@@ -1390,6 +1391,12 @@ void gfc3d_IPM(GlobalFrictionContactProblem* restrict problem, double* restrict
free(r_check);
free(Qp_primalConstraint);
+ jacobian_is_nan = NM_isnan(J);
+ if (jacobian_is_nan)
+ {
+ numerics_printf_verbose(0, "The Jacobian matrix J contains NaN");
+ break;
+ }
NSM_linearSolverParams(J)->solver = NSM_HSL;
if ( options->iparam[SICONOS_FRICTION_3D_IPM_IPARAM_REFINEMENT] == SICONOS_FRICTION_3D_IPM_IPARAM_REFINEMENT_YES )
@@ -1401,6 +1408,8 @@ void gfc3d_IPM(GlobalFrictionContactProblem* restrict problem, double* restrict
}
else
{
+
+
NM_LDLT_solve(J, rhs, 1);
}
@@ -1564,6 +1573,12 @@ void gfc3d_IPM(GlobalFrictionContactProblem* restrict problem, double* restrict
printf("ERROR\n");
}
}
+ if (jacobian_is_nan)
+ {
+ hasNotConverged = 2;
+ J = NM_free(J);
+ break;
+ }
/* computing the affine step-length */
alpha_primal = getStepLength(velocity, d_velocity, nd, n, 1.0); //gmm);
@@ -2206,7 +2221,7 @@ void gfc3d_ipm_set_default(SolverOptions* options)
options->iparam[SICONOS_FRICTION_3D_IPM_IPARAM_UPDATE_S] = 0;
- options->iparam[SICONOS_FRICTION_3D_IPM_IPARAM_LS_FORM] = SICONOS_FRICTION_3D_IPM_IPARAM_LS_2X2_QPH;
+ options->iparam[SICONOS_FRICTION_3D_IPM_IPARAM_LS_FORM] = SICONOS_FRICTION_3D_IPM_IPARAM_LS_3X3_QPH;
//options->iparam[SICONOS_FRICTION_3D_IPARAM_RESCALING] = SICONOS_FRICTION_3D_RESCALING_BALANCING_MHHT;
options->iparam[SICONOS_FRICTION_3D_IPARAM_RESCALING] = SICONOS_FRICTION_3D_RESCALING_NO;
