# Changeset 5569

Ignore:
Timestamp:
Nov 28, 2013 11:46:14 AM (6 years ago)
Message:

Added new heuristics for regularization. If the scaling factors are set to max then calculate the condition number and if the condition number is larger than regularization_tolerance (parameter), then regularize. Related to ticket:3275

Location:
trunk/RuntimeLibrary/src/jmi
Files:
2 edited

Unmodified
Removed
• ## trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.c

 r5566 } static realtype jmi_calculate_condition_number(jmi_block_solver_t * block, realtype* A) { jmi_kinsol_solver_t* solver = block->solver; char norm = 'I'; int N = block->n; double Jnorm = 1.0, Jcond = 1.0; int info; dgecon_(&norm, &N, A, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info); return 1.0/Jcond; } static int jmi_kin_lsetup(struct KINMemRec * kin_mem) { jmi_block_solver_t *block = kin_mem->kin_user_data; /* Evaluate Jacobian */ ret = kin_dF(N, solver->kin_y, kin_mem->kin_fval, solver->J, block, kin_mem->kin_vtemp1, kin_mem->kin_vtemp2); solver->updated_jacobian_flag = 1; /* The Jacobian is current */ if(ret != 0 ) return ret; jmi_kinsol_reg_matrix(block); dgetrf_(  &N, &N, solver->JTJ->data, &N, solver->lapack_ipiv, &info); } else { } else { /* if (solver->using_max_min_scaling_flag) { realtype cond = jmi_calculate_condition_number(block, solver->J->data); jmi_log_node(block->log, logWarning, "JacobianConditioningNumber", " large values may lead to convergence problems.", cond); } */ solver->J_is_singular_flag = 0; } realtype*  bd = N_VGetArrayPointer(b); realtype*  xd = N_VGetArrayPointer(x); jmi_log_node_t node; int N = block->n; char trans = 'N'; int ret, i; solver->updated_jacobian_flag = 0; /* The Jacobian is no longer current */ if(solver->force_new_J_flag) { sfdotJp is the dot product of the scaled f vector and the scaled vector J*p, where the scaling uses fscale. */ if((block->callbacks->log_options.log_level >= 6)) { node = jmi_log_enter_fmt(block->log, logInfo, "KinsolLinearSolver", "Solving the linear system in ", block->id); } kin_mem->kin_sJpnorm = N_VWL2Norm(b,solver->kin_f_scale); dgetrs_(&trans, &N, &i, solver->JTJ->data, &N, solver->lapack_ipiv, xd, &N, &ret); solver->force_new_J_flag = 1; if((block->callbacks->log_options.log_level >= 6)) { jmi_log_real_matrix(block->log, node, logInfo, "jacobian", solver->JTJ->data, N, N); jmi_log_leave(block->log, node); } } else { N_VScale(ONE, b, x); i = 1; if((block->callbacks->log_options.log_level >= 6)) { jmi_log_real_matrix(block->log, node, logInfo, "jacobian", solver->J_LU->data, N, N); jmi_log_reals(block->log, node, logInfo, "b", xd, N); } dgetrs_(&trans, &N, &i, solver->J_LU->data, &N, solver->lapack_ipiv, xd, &N, &ret); if((block->callbacks->log_options.log_level >= 6)) { jmi_log_reals(block->log, node, logInfo, "x", xd, N); jmi_log_leave(block->log, node); } } if(use_scaling_flag) { solver->using_max_min_scaling_flag = 0; /* NOT using max/min scaling */ /* check that scaling factors has reasonable magnitude */ for(i = 0; i < N; i++) { if(scale_ptr[i] < tol) { scale_ptr[i] = 1/tol; /* Singular Jacobian? */ solver->using_max_min_scaling_flag = 1; /* Using maximum scaling */ jmi_log_node(block->log, logWarning, "Warning", "Using maximum scaling factor in , " " Consider rescaling in the model or tighter tolerance.", block->id, i); else if(scale_ptr[i] > 1/tol) { scale_ptr[i] = tol; solver->using_max_min_scaling_flag = 1; /* Using minimum scaling */ jmi_log_node(block->log, logWarning, "Warning", "Using minimal scaling factor in , " " Consider rescaling in the model or tighter tolerance.", block->id, i); } } if (solver->using_max_min_scaling_flag) { realtype cond = jmi_calculate_condition_number(block, solver->J->data); jmi_log_node(block->log, logInfo, "ConditionNumber", "Calculated condition number in . Regularizing if greater than ", block->id, cond, solver->kin_reg_tol); if (cond > solver->kin_reg_tol) { int info; jmi_kinsol_reg_matrix(block); dgetrf_(  &block->n, &block->n, solver->JTJ->data, &block->n, solver->lapack_ipiv, &info); solver->J_is_singular_flag = 1; } } /* estimate condition number of the scaled jacobian and scale function tolerance with it. */ solver->kin_ftol = block->options->min_tol; solver->kin_stol = block->options->min_tol; solver->using_max_min_scaling_flag = 0; /* Not using max/min scaling */ solver->J = NewDenseMat(n ,n);
• ## trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.h

 r5566 int use_steepest_descent_flag;  /**< \brief A flag indicating that steepest descent and not Newton direction should be used */ int force_new_J_flag;           /**< \brief A flag indicating that J needs to be recalculated */ int using_max_min_scaling_flag; /**< \brief A flag indicating if either the maximum scaling is used of the minimum */ int updated_jacobian_flag;      /**< \brief A flag indicating if an updated Jacobian is used to solve the system */ DlsMat J_LU;                    /**< \brief Jacobian matrix/it's LU decomposition */ DlsMat J_scale;                 /**< \brief Jacobian matrix scaled with xnorm for used for fnorm calculation */
Note: See TracChangeset for help on using the changeset viewer.