Changeset 9636


Ignore:
Timestamp:
Feb 13, 2017 7:08:05 PM (3 years ago)
Author:
Christian Andersson
Message:

Minor cleanup to Kinsol. Related to ticket:4815

File:
1 edited

Legend:

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

    r9597 r9636  
    16451645
    16461646/* Perform LU factorization with different linear algebra packages */
    1647 static int jmi_LU_factorization(jmi_block_solver_t * block, DlsMat matrix, int lin_alg_package) {
     1647static int jmi_LU_factorization(jmi_block_solver_t * block, DlsMat matrix) {
    16481648    jmi_kinsol_solver_t* solver = block->solver;
    16491649    int info = 0, N = block->n;
     1650    int lin_alg_package = block->options->experimental_mode & jmi_block_solver_experimental_LU_through_sundials ? 1:0;
     1651   
    16501652    if(lin_alg_package == 0) {
    16511653        dgetrf_(  &N, &N, matrix->data, &N, solver->lapack_ipiv, &info);
     
    16581660
    16591661/* Solve with an LU factorized matrix with different linear algebra packages */
    1660 static int jmi_LU_solve(jmi_block_solver_t * block, DlsMat matrix, realtype* xd, int lin_alg_package) {
     1662static int jmi_LU_solve(jmi_block_solver_t * block, DlsMat matrix, realtype* xd) {
    16611663    jmi_kinsol_solver_t* solver = block->solver;
    16621664    int ret = 0, N = block->n;
     1665    int lin_alg_package = block->options->experimental_mode & jmi_block_solver_experimental_LU_through_sundials ? 1:0;
     1666   
    16631667    if(lin_alg_package == 1) {
    16641668        DenseGETRS(matrix, solver->sundials_permutationwork, xd);
     
    16881692        if (cond > solver->kin_reg_tol) {
    16891693            if(N > 1 && solver->handling_of_singular_jacobian_flag == JMI_REGULARIZATION) {
    1690                 int info;
    16911694                jmi_kinsol_reg_matrix(block);
    1692                 dgetrf_(  &block->n, &block->n, solver->JTJ->data, &block->n, solver->lapack_ipiv, &info);
     1695                jmi_LU_factorization(block, solver->JTJ);
    16931696            }
    16941697            solver->J_is_singular_flag = 1;
     
    17331736    DenseCopy(block->J, solver->J_LU);
    17341737    /* Perform LU factorization to be used with dgecon */
    1735     dgetrf_(&N, &N, solver->J_LU->data, &N, solver->lapack_ipiv, &info);
     1738    info = jmi_LU_factorization(block, solver->J_LU);
    17361739    if (info != 0 ) {
    17371740        /* If matrix i singular, return something very large to be evaluated*/
     
    17401743    /* Compute reciprocal condition number */
    17411744    dgecon_(&norm, &N, solver->J_LU->data, &N, &J_norm, &J_recip_cond, solver->lapack_work, solver->lapack_iwork,&info);
    1742     /* To be evaluated - why is this needed? Error handling due to J being used instead of J_LU?
    1743     if(Jcond < 0) Jcond = -Jcond;
    1744     if(Jcond == 0.0) Jcond = 1e-30;
    1745     */
     1745
    17461746    return 1.0/J_recip_cond;
    17471747}
     
    18471847        jmi_log_node_t node = jmi_log_enter_fmt(block->log, logInfo, "BroydenJacobianUpdate", "<block:%s>", block->label);
    18481848        if (block->callbacks->log_options.log_level >= 6) {
    1849 #if 0
    1850             jmi_log_reals(block->log, node, logInfo, "last_residual", N_VGetArrayPointer(solver->last_residual), N);
    1851             jmi_log_reals(block->log, node, logInfo, "current_residual", N_VGetArrayPointer(b), N);
    1852             jmi_log_reals(block->log, node, logInfo, "last_step", N_VGetArrayPointer(kin_mem->kin_pp), N);
    1853             jmi_log_fmt(block->log, node,logInfo, "Update denominator <step_norm: %g>", denom);
    1854 #endif
    18551849            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", block->J->data, N, N);
    18561850        }
     
    19051899        jmi_log_node_t node = jmi_log_enter_fmt(block->log, logInfo, "SparseBroydenJacobianUpdate", "<block:%s>", block->label);
    19061900        if (block->callbacks->log_options.log_level >= 6) {
    1907 #if 0
    1908             jmi_log_reals(block->log, node, logInfo, "last_residual", N_VGetArrayPointer(solver->last_residual), N);
    1909             jmi_log_reals(block->log, node, logInfo, "current_residual", N_VGetArrayPointer(b), N);
    1910             jmi_log_reals(block->log, node, logInfo, "last_step", N_VGetArrayPointer(kin_mem->kin_pp), N);
    1911             jmi_log_reals(block->log, node,logInfo, "denoms", N_VGetArrayPointer(solver->work_vector2), N);
    1912             jmi_log_reals(block->log, node,logInfo, "res_JacStep", N_VGetArrayPointer(solver->work_vector), N);
    1913 #endif
    19141901            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", block->J->data, N, N);
    19151902        }
     
    19591946        jmi_log_node_t node = jmi_log_enter_fmt(block->log, logInfo, "ModifiedBFGSJacobianUpdate", "<block:%s>", block->label);
    19601947        if (block->callbacks->log_options.log_level >= 6) {
    1961 #if 0
    1962             jmi_log_reals(block->log, node, logInfo, "last_residual", N_VGetArrayPointer(solver->last_residual), N);
    1963             jmi_log_reals(block->log, node, logInfo, "Jac_step", N_VGetArrayPointer(solver->work_vector), N);
    1964             jmi_log_reals(block->log, node, logInfo, "res_diff_f_scale2_Jac", N_VGetArrayPointer(solver->work_vector2), N);
    1965             jmi_log_reals(block->log, node, logInfo, "current_residual", N_VGetArrayPointer(b), N);
    1966             jmi_log_reals(block->log, node, logInfo, "last_step", N_VGetArrayPointer(kin_mem->kin_pp), N);
    1967             jmi_log_reals(block->log, node, logInfo, "y_scale", N_VGetArrayPointer(solver->kin_y_scale), N);
    1968             jmi_log_reals(block->log, node, logInfo, "f_scale", N_VGetArrayPointer(solver->kin_f_scale), N);
    1969             jmi_log_fmt(block->log, node,logInfo, "Update denominator <step_norm: %g>", denom1);
    1970             jmi_log_fmt(block->log, node,logInfo, "Update denominator <denom2: %g>", denom2);
    1971 #endif
    19721948            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", block->J->data, N, N);
    19731949        }
     
    20121988    }
    20131989   
    2014     info = jmi_LU_factorization(block, solver->J_LU,block->options->experimental_mode & jmi_block_solver_experimental_LU_through_sundials ? 1:0);
     1990    info = jmi_LU_factorization(block, solver->J_LU);
    20151991    if(info != 0) {
    20161992        if(jmi_kinsol_zero_column_jacobian_handling(block)) {
    20171993            DenseCopy(block->J, solver->J_LU); /* make a copy of the Jacobian that will be used for LU factorization */
    2018             info = jmi_LU_factorization(block, solver->J_LU,block->options->experimental_mode & jmi_block_solver_experimental_LU_through_sundials ? 1:0);
     1994            info = jmi_LU_factorization(block, solver->J_LU);
    20191995        }
    20201996    }
     
    20482024                if(N > 1) {
    20492025                    jmi_kinsol_reg_matrix(block);
    2050                     info = jmi_LU_factorization(block, solver->JTJ, block->options->experimental_mode & jmi_block_solver_experimental_LU_through_sundials? 1:0);
     2026                    info = jmi_LU_factorization(block, solver->JTJ);
    20512027                }
    20522028            } else if (solver->handling_of_singular_jacobian_flag == JMI_MINIMUM_NORM) {
     
    20932069    solver->current_nni = nniters;
    20942070
    2095     if(block->options->jacobian_update_mode == jmi_broyden_jacobian_update_mode) {
    2096         if(!solver->updated_jacobian_flag) {
    2097             /* Never do update in the first iteration */
    2098             if (nniters > 0) {       
    2099                 ret = jmi_kin_make_Broyden_update(block, b);
    2100                 if(ret != 0) return ret;
    2101             }
    2102         }
    2103     }
    2104 
    2105     if(block->options->experimental_mode & jmi_block_solver_experimental_use_modifiedBFGS) {
    2106         if(!solver->updated_jacobian_flag) {
    2107             /* Never do update in the first iteration */
    2108             if (nniters > 0) {       
    2109                 ret = jmi_kin_make_modifiedBFGS_update(block, b);
    2110                 if(ret != 0) return ret;
    2111             }
    2112         }
    2113     }
    2114 
    2115     if(block->options->experimental_mode & jmi_block_solver_experimental_Sparse_Broyden) {
    2116         if(!solver->updated_jacobian_flag) {
    2117             /* Never do update in the first iteration */
    2118             if (nniters > 0) {       
    2119                 ret = jmi_kin_make_sparse_Broyden_update(block, b);
    2120                 if(ret != 0) return ret;
    2121             }
    2122         }
     2071   
     2072    if (!solver->updated_jacobian_flag && nniters > 0) {
     2073        if(block->options->jacobian_update_mode == jmi_broyden_jacobian_update_mode)
     2074            ret = jmi_kin_make_Broyden_update(block, b);
     2075        if (block->options->experimental_mode & jmi_block_solver_experimental_use_modifiedBFGS)
     2076            ret = jmi_kin_make_modifiedBFGS_update(block, b);
     2077        if (block->options->experimental_mode & jmi_block_solver_experimental_Sparse_Broyden)
     2078            ret = jmi_kin_make_sparse_Broyden_update(block, b);
     2079
     2080        if(ret != 0) return ret;
    21232081    }
    21242082
     
    22392197            /* Back-solve and get solution in x */
    22402198            t = jmi_block_solver_start_clock(block);
    2241             ret = jmi_LU_solve(block, solver->JTJ, xd, block->options->experimental_mode & jmi_block_solver_experimental_LU_through_sundials ? 1:0);
     2199            ret = jmi_LU_solve(block, solver->JTJ, xd);
    22422200            kin_char_log(solver, 'r');
    22432201            solver->force_new_J_flag = 1;
     
    23052263        }
    23062264       
    2307         ret = jmi_LU_solve(block, solver->J_LU, xd, block->options->experimental_mode & jmi_block_solver_experimental_LU_through_sundials ? 1:0);
     2265        ret = jmi_LU_solve(block, solver->J_LU, xd);
    23082266
    23092267        if((block->callbacks->log_options.log_level >= 6)) {
     
    28262784    }
    28272785   
    2828     /* Brent is called for 1D to get higher accuracy. It is called independently on the KINSOL success */
    2829     if((block->n == 1) && block->options->use_Brent_in_1d_flag) {
    2830         jmi_log_node(log, logInfo, "Brent", "Trying Brent's method in <block: %s>", block->label);
    2831         if(( solver->f_pos_min_1d != BIG_REAL) &&
    2832                 ( solver->f_neg_max_1d != -BIG_REAL)) {
    2833            
    2834             realtype u, f;
    2835             if(!jmi_brent_search(brentf, solver->y_neg_max_1d,  solver->y_pos_min_1d,
    2836                                  solver->f_neg_max_1d,  solver->f_pos_min_1d, 0, &u, &f,block)) {
    2837                 block->x[0] = u;
    2838                 flag = KIN_SUCCESS;
    2839             }               
    2840         }
    2841         if(flag != KIN_SUCCESS) {
    2842             jmi_log_node(log, logError, "Error", "Could neither iterate to required accuracy "
    2843                          "nor bracket the root of 1D equation in <block: %s>", block->label);
    2844         }
    2845     }
    2846    
     2786    /* If Kinsol failed, try to solve with steepest descent instead */
    28472787    if((block->options->experimental_mode & jmi_block_solver_experimental_steepest_descent) &&
    28482788        (flag != KIN_SUCCESS)) {
    2849         /* try to solve with steepest descent instead */
    2850 
     2789       
    28512790        jmi_log_node(log, logInfo, "Progress", "<source:%s><block:%s><message:%s>",
    28522791                     "jmi_kinsol_solver", block->label, "Attempting steepest descent iterations");       
     
    29142853                }
    29152854            }
    2916 #ifdef JMI_KINSOL_PRINT_ON_FAIL
    2917             {
    2918                 realtype* x = block->x;
    2919                 int i;
    2920                 struct KINMemRec * kin_mem = solver->kin_mem;
    2921                
    2922                 printf("Could not converge in block %s,KINSOL error:%s,fnorm=%g,stol=%g,ftol=%g:\n"
    2923                        "x = N.array([\n",block->label, jmi_kinsol_flag_to_name(flag), fnorm, solver->kin_stol,solver->kin_ftol);
    2924                 for (i=0;i<block->n;i++) {
    2925                     printf("%.16e\n",x[i]);
    2926                 }
    2927                 printf("]);\n");
    2928                 printf("x_scale = N.array([\n");
    2929                 for (i=0;i<block->n;i++) {
    2930                     printf("%.16e\n",Ith(solver->kin_y_scale,i));
    2931                 }
    2932                 printf("]);\n");
    2933                 printf("f = N.array([\n");
    2934                 for (i=0;i<block->n;i++) {
    2935                     printf("%.16e\n",Ith(kin_mem->kin_fval,i));
    2936                 }
    2937                 printf("]);\n");
    2938                 N_VScale( 1+2*UNIT_ROUNDOFF,solver->kin_y,solver->kin_y);
    2939                 kin_f(solver->kin_y,kin_mem->kin_fval,block);
    2940                 printf("f_inc = N.array([\n");
    2941                 for (i=0;i<block->n;i++) {
    2942                     printf("%.16e\n",Ith(kin_mem->kin_fval,i));
    2943                 }
    2944                 printf("]);\n");
    2945                 N_VScale( 1-4*UNIT_ROUNDOFF,solver->kin_y,solver->kin_y);
    2946                 kin_f(solver->kin_y,kin_mem->kin_fval,block);
    2947                 printf("f_dec = N.array([\n");
    2948                 for (i=0;i<block->n;i++) {
    2949                     printf("%.16e\n",Ith(kin_mem->kin_fval,i));
    2950                 }
    2951                 printf("]);\n");
    2952 
    2953                 printf("f_scale = N.array([\n");
    2954                 for (i=0;i<block->n;i++) {
    2955                     printf("%.16e\n",Ith(solver->kin_f_scale,i));
    2956                 }
    2957                 printf("]);\n");
    2958             }
    2959 #endif
    29602855        }
    29612856    }
Note: See TracChangeset for help on using the changeset viewer.