Changeset 5569


Ignore:
Timestamp:
Nov 28, 2013 11:46:14 AM (6 years ago)
Author:
Christian Andersson
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

Legend:

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

    r5566 r5569  
    673673}
    674674
     675static realtype jmi_calculate_condition_number(jmi_block_solver_t * block, realtype* A) {
     676    jmi_kinsol_solver_t* solver = block->solver;
     677    char norm = 'I';
     678    int N = block->n;
     679    double Jnorm = 1.0, Jcond = 1.0;
     680    int info;
     681   
     682    dgecon_(&norm, &N, A, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);
     683   
     684    return 1.0/Jcond;
     685}
     686
    675687static int jmi_kin_lsetup(struct KINMemRec * kin_mem) {
    676688    jmi_block_solver_t *block = kin_mem->kin_user_data;
     
    685697    /* Evaluate Jacobian */
    686698    ret = kin_dF(N, solver->kin_y, kin_mem->kin_fval, solver->J, block, kin_mem->kin_vtemp1, kin_mem->kin_vtemp2);
     699    solver->updated_jacobian_flag = 1; /* The Jacobian is current */
    687700   
    688701    if(ret != 0 ) return ret;
     
    715728        jmi_kinsol_reg_matrix(block);
    716729        dgetrf_(  &N, &N, solver->JTJ->data, &N, solver->lapack_ipiv, &info);
    717     }
    718     else {
     730    } else {
     731        /* if (solver->using_max_min_scaling_flag) {
     732            realtype cond = jmi_calculate_condition_number(block, solver->J->data);
     733            jmi_log_node(block->log, logWarning, "JacobianConditioningNumber",
     734                             "<JacobianConditionEstimate:%E> large values may lead to convergence problems.", cond);
     735        }
     736        */
    719737        solver->J_is_singular_flag = 0;       
    720738    }
     
    735753    realtype*  bd = N_VGetArrayPointer(b);
    736754    realtype*  xd = N_VGetArrayPointer(x);
     755    jmi_log_node_t node;
    737756   
    738757    int N = block->n;
    739758    char trans = 'N';
    740759    int ret, i;
     760   
     761    solver->updated_jacobian_flag = 0; /* The Jacobian is no longer current */
    741762   
    742763    if(solver->force_new_J_flag) {       
     
    756777       sfdotJp is the dot product of the scaled f vector and the scaled
    757778       vector J*p, where the scaling uses fscale. */
    758  
     779    if((block->callbacks->log_options.log_level >= 6)) {
     780        node = jmi_log_enter_fmt(block->log, logInfo, "KinsolLinearSolver", "Solving the linear system in <block:%d>", block->id);
     781    }
    759782 
    760783    kin_mem->kin_sJpnorm = N_VWL2Norm(b,solver->kin_f_scale);
     
    813836        dgetrs_(&trans, &N, &i, solver->JTJ->data, &N, solver->lapack_ipiv, xd, &N, &ret);
    814837        solver->force_new_J_flag = 1;
     838       
     839        if((block->callbacks->log_options.log_level >= 6)) {
     840            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", solver->JTJ->data, N, N);
     841            jmi_log_leave(block->log, node);
     842        }
    815843    }
    816844    else {
    817845        N_VScale(ONE, b, x);
    818846        i = 1;
     847       
     848        if((block->callbacks->log_options.log_level >= 6)) {
     849            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", solver->J_LU->data, N, N);
     850            jmi_log_reals(block->log, node, logInfo, "b", xd, N);
     851        }
     852       
    819853        dgetrs_(&trans, &N, &i, solver->J_LU->data, &N, solver->lapack_ipiv, xd, &N, &ret);
     854       
     855        if((block->callbacks->log_options.log_level >= 6)) {
     856            jmi_log_reals(block->log, node, logInfo, "x", xd, N);
     857            jmi_log_leave(block->log, node);
     858        }
    820859    }
    821860   
     
    900939   
    901940    if(use_scaling_flag) {
     941        solver->using_max_min_scaling_flag = 0; /* NOT using max/min scaling */
    902942        /* check that scaling factors has reasonable magnitude */
    903943        for(i = 0; i < N; i++) {
    904944            if(scale_ptr[i] < tol) {
    905945                scale_ptr[i] = 1/tol; /* Singular Jacobian? */
     946                solver->using_max_min_scaling_flag = 1; /* Using maximum scaling */
    906947                jmi_log_node(block->log, logWarning, "Warning", "Using maximum scaling factor in <block: %d>, "
    907948                             "<equation: %d> Consider rescaling in the model or tighter tolerance.", block->id, i);
     
    909950            else if(scale_ptr[i] > 1/tol) {
    910951                scale_ptr[i] = tol;
     952                solver->using_max_min_scaling_flag = 1; /* Using minimum scaling */
    911953                jmi_log_node(block->log, logWarning, "Warning", "Using minimal scaling factor in <block: %d>, "
    912954                             "<equation: %d> Consider rescaling in the model or tighter tolerance.", block->id, i);
     
    925967        }
    926968    }
     969   
     970    if (solver->using_max_min_scaling_flag) {
     971        realtype cond = jmi_calculate_condition_number(block, solver->J->data);
     972       
     973        jmi_log_node(block->log, logInfo, "ConditionNumber",
     974                         "Calculated condition number in <block: %d>. Regularizing if <cond: %E> greater than <regtol: %E>", block->id, cond, solver->kin_reg_tol);
     975        if (cond > solver->kin_reg_tol) {
     976            int info;
     977            jmi_kinsol_reg_matrix(block);
     978            dgetrf_(  &block->n, &block->n, solver->JTJ->data, &block->n, solver->lapack_ipiv, &info);
     979            solver->J_is_singular_flag = 1;
     980        }
     981    }
     982   
    927983    /* estimate condition number of the scaled jacobian
    928984        and scale function tolerance with it. */
     
    9861042    solver->kin_ftol = block->options->min_tol;
    9871043    solver->kin_stol = block->options->min_tol;
     1044    solver->using_max_min_scaling_flag = 0; /* Not using max/min scaling */
    9881045   
    9891046    solver->J = NewDenseMat(n ,n);
  • trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.h

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