Changeset 9277


Ignore:
Timestamp:
Sep 27, 2016 5:21:14 PM (3 years ago)
Author:
efredriksson
Message:

#4815 Refactored out the "IllConditionedJacobian" logging such that it's in its own function.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/20160824.x/RuntimeLibrary/src/jmi/jmi_kinsol_solver.c

    r9206 r9277  
    744744}
    745745
     746/* Local helper that will log the condition number of the Jacobian */
     747static void jmi_kinsol_log_jacobian_cond_nbr(jmi_block_solver_t *block) {
     748    jmi_kinsol_solver_t* solver = block->solver;
     749    int info, N = block->n;
     750    realtype tol = solver->kin_stol;
     751
     752    dgetrf_(&N, &N, block->J_scale->data, &N, solver->lapack_iwork, &info);
     753    if(info > 0) {
     754        jmi_log_node(block->log, logWarning, "SingularJacobian",
     755            "Singular Jacobian detected when checking condition number in <block:%s>. Solver may fail to converge.", block->label);
     756    }
     757    else {
     758        char norm = 'I';
     759        double Jnorm = 1.0, Jcond = 1.0;
     760        dgecon_(&norm, &N, block->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
     761
     762        if(tol * Jcond < UNIT_ROUNDOFF) {
     763            jmi_log_node_t node = jmi_log_enter_fmt(block->log, logWarning, "IllConditionedJacobian",
     764                "<JacobianInverseConditionEstimate:%E> Solver may fail to converge in <block: %s>.", Jcond, block->label);
     765            if (block->callbacks->log_options.log_level >= 4) {
     766                jmi_log_reals(block->log, node, logWarning, "ivs", N_VGetArrayPointer(solver->kin_y), block->n);
     767            }
     768            jmi_log_leave(block->log, node);
     769
     770        }
     771        else {
     772            jmi_log_node(block->log, logInfo, "JacobianCondition",
     773                "<JacobianInverseConditionEstimate:%E>", Jcond);
     774        }
     775    }
     776}
     777
    746778/* Logging callback for KINSOL used to report on errors during solution */
    747779void kin_err(int err_code, const char *module, const char *function, char *msg, void *eh_data){
     
    779811
    780812        if(block->options->check_jac_cond_flag) {
    781             char norm = 'I';
    782             double Jnorm = 1.0, Jcond = 1.0;
    783             int i, info, N = block->n;
     813            int i, N = block->n;
    784814            realtype tol = solver->kin_stol;
    785815            realtype *scale_ptr = N_VGetArrayPointer(block->f_scale);
     
    800830                DenseCopy(block->J, block->J_scale);
    801831            }
    802             dgetrf_(  &N, &N, block->J_scale->data, &N, solver->lapack_iwork, &info);
    803             if(info > 0) {
    804                 jmi_log_node(block->log, logWarning, "SingularJacobian",
    805                     "Singular Jacobian detected when checking condition number in <block:%s>. Solver may fail to converge.", block->label);
    806             }
    807             else {
    808                 dgecon_(&norm, &N, block->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
    809 
    810                 if(tol * Jcond < UNIT_ROUNDOFF) {
    811                     jmi_log_node_t node = jmi_log_enter_fmt(block->log, logWarning, "IllConditionedJacobian",
    812                         "<JacobianInverseConditionEstimate:%E> Solver may fail to converge in <block: %s>.", Jcond, block->label);
    813                     if (block->callbacks->log_options.log_level >= 4) {
    814                         jmi_log_reals(block->log, node, logWarning, "ivs", N_VGetArrayPointer(solver->kin_y), block->n);
    815                     }
    816                     jmi_log_leave(block->log, node);
    817 
    818                 }
    819                 else {
    820                     jmi_log_node(block->log, logInfo, "JacobianCondition",
    821                         "<JacobianInverseConditionEstimate:%E>", Jcond);
    822                 }
    823             }
     832           
     833            jmi_kinsol_log_jacobian_cond_nbr(block);
    824834        }
    825835    }
     
    16581668    jmi_block_solver_options_t* bsop = block->options;
    16591669    int use_scaling_flag = bsop->residual_equation_scaling_mode;
    1660     realtype tol = solver->kin_stol;
    16611670    realtype* scale_ptr = N_VGetArrayPointer(block->f_scale);
    16621671
     
    16811690    if((N > 1) && bsop->check_jac_cond_flag){
    16821691        realtype* scaled_col_ptr;
    1683         char norm = 'I';
    1684         double Jnorm = 1.0, Jcond = 1.0;
    1685         int info;
     1692
    16861693        if(use_scaling_flag) {
    16871694            for(i = 0; i < N; i++){
     
    16961703        }
    16971704
    1698         dgetrf_(  &N, &N, block->J_scale->data, &N, solver->lapack_iwork, &info);
    1699         if(info > 0) {
    1700             jmi_log_node(block->log, logWarning, "SingularJacobian",
    1701                 "Singular Jacobian detected when checking condition number in <block:%s>. Solver may fail to converge.", block->label);
    1702         }
    1703         else {
    1704             dgecon_(&norm, &N, block->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
    1705 
    1706             if(tol * Jcond < UNIT_ROUNDOFF) {
    1707                 jmi_log_node_t node = jmi_log_enter_fmt(block->log, logWarning, "IllConditionedJacobian",
    1708                     "<JacobianInverseConditionEstimate:%E> Solver may fail to converge in <block: %s>.", Jcond, block->label);
    1709                 if (block->callbacks->log_options.log_level >= 4) {
    1710                     jmi_log_reals(block->log, node, logWarning, "ivs", N_VGetArrayPointer(solver->kin_y), block->n);
    1711                 }
    1712                 jmi_log_leave(block->log, node);
    1713 
    1714             }
    1715             else {
    1716                 jmi_log_node(block->log, logInfo, "JacobianCondition",
    1717                     "<JacobianInverseConditionEstimate:%E>", Jcond);
    1718             }
    1719         }
     1705        jmi_kinsol_log_jacobian_cond_nbr(block);
    17201706    }
    17211707}
     
    27452731
    27462732        if(block->options->check_jac_cond_flag) {
    2747             int i, info, N = block->n;
     2733            int i, N = block->n;
    27482734            realtype* scale_ptr = N_VGetArrayPointer(block->f_scale);
    27492735            realtype tol = solver->kin_stol;
     
    27652751            }
    27662752
    2767             dgetrf_(  &N, &N, block->J_scale->data, &N, solver->lapack_iwork, &info);
    2768             if(info > 0) {
    2769                 jmi_log_node(block->log, logWarning, "SingularJacobian",
    2770                     "Singular Jacobian detected when checking condition number in <block:%s>. Solver may fail to converge.", block->label);
    2771             }
    2772             else {
    2773                 char norm = 'I';
    2774                 double Jnorm = 1.0, Jcond = 1.0;
    2775                 dgecon_(&norm, &N, block->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
    2776 
    2777                 if(tol * Jcond < UNIT_ROUNDOFF) {
    2778                     jmi_log_node_t node = jmi_log_enter_fmt(block->log, logWarning, "IllConditionedJacobian",
    2779                         "<JacobianInverseConditionEstimate:%E> Solver may fail to converge in <block: %s>.", Jcond, block->label);
    2780                     if (block->callbacks->log_options.log_level >= 4) {
    2781                         jmi_log_reals(block->log, node, logWarning, "ivs", N_VGetArrayPointer(solver->kin_y), block->n);
    2782                     }
    2783                     jmi_log_leave(block->log, node);
    2784 
    2785                 }
    2786                 else {
    2787                     jmi_log_node(block->log, logInfo, "JacobianCondition",
    2788                         "<JacobianInverseConditionEstimate:%E>", Jcond);
    2789                 }
    2790             }
     2753            jmi_kinsol_log_jacobian_cond_nbr(block);
    27912754        }
    27922755    }
Note: See TracChangeset for help on using the changeset viewer.