Changeset 9581


Ignore:
Timestamp:
Jan 20, 2017 4:09:03 PM (3 years ago)
Author:
Christian Andersson
Message:

Moved vector operations to use the ones implemented on the NVector module. Related to ticket:4815

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

Legend:

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

    r9293 r9581  
    313313}
    314314
    315 static double kin_dot_product(double* y, double* x, int n) {
    316     double prod = 0;
    317     int i;
    318     for(i=0; i<n; i++)
    319         prod+= y[i]*x[i];
    320     return prod;
    321 }
    322 
    323 static void kin_sum_vectors(double* y, double* x, int n, double* z) {
    324     int i;
    325     for(i=0; i<n; i++)
    326         z[i] = y[i]+ x[i];
    327 }
    328 
    329 static void kin_prod_vectors(double* y, double* x, int n, double* z) {
    330     int i;
    331     for(i=0; i<n; i++)
    332         z[i] = y[i]*x[i];
    333 }
    334 
    335315int jmi_kin_setup_column_partition(jmi_block_solver_t * block) {
    336316    jmi_kinsol_solver_t* solver = (jmi_kinsol_solver_t*)block->solver;
     
    338318    int leftColumns = N;
    339319    int i, group = 0;
    340     realtype* rows_r = N_VGetArrayPointer(solver->work_vector2);
    341     N_VConst_Serial(0.0, solver->work_vector);
     320    N_Vector col_i = N_VCloneEmpty(solver->work_vector);
     321
     322    N_VConst(0.0, solver->work_vector);
    342323    while(leftColumns > 0 && group <N) {
    343324        group++;
    344         N_VConst_Serial(0.0, solver->work_vector2);
     325        N_VConst(0.0, solver->work_vector2);
    345326        for(i = 0; i<N; i++) {
    346327            if(Ith(solver->work_vector, i) == 0.0) {
    347                 if(kin_dot_product(rows_r, solver->J_Dependency->cols[i], N)==0.0) { /* Column can be added */
    348                     kin_sum_vectors(rows_r, solver->J_Dependency->cols[i], N, rows_r);
     328                N_VSetArrayPointer(DENSE_COL(solver->J_Dependency, i), col_i);
     329               
     330                if(N_VDotProd(solver->work_vector2, col_i) == 0.0) { /* Column can be added */                 
     331                    N_VLinearSum(1.0, solver->work_vector2, 1.0, col_i, solver->work_vector2);
     332
    349333                    solver->jac_compression_groups[N-leftColumns]=group;
    350334                    solver->jac_compression_group_index[N-leftColumns]=i;
     
    355339        }
    356340    }
     341    N_VDestroy(col_i);
     342   
    357343    if(leftColumns != 0) {
    358344        jmi_log_node(block->log, logWarning, "ColumnPartitioning", "Column partitioning error, dependency data may be corrupt.");
     
    638624                    N_VSetArrayPointer(DENSE_COL(J, solver->jac_compression_group_index[k]), jthCol);
    639625                    N_VLinearSum(inc_inv, ftemp, -inc_inv, fu, jthCol);
    640                     kin_prod_vectors(solver->J_Dependency->cols[solver->jac_compression_group_index[k]], N_VGetArrayPointer(jthCol), N, N_VGetArrayPointer(jthCol));
    641                 }
     626                   
     627                    N_VSetArrayPointer(DENSE_COL(solver->J_Dependency, solver->jac_compression_group_index[k]), solver->work_vector3);
     628                    N_VProd(solver->work_vector3, jthCol, jthCol);
     629                   
     630                }
    642631                /* Restore original array pointer in tmp2 */
    643632                N_VSetArrayPointer(tmp2_data, tmp2);
     
    18001789
    18011790    if(block->force_rescaling) {
    1802         if(block->options->residual_equation_scaling_mode != jmi_residual_scaling_none)
    1803             kin_char_log(solver, 's');
     1791        if(block->options->residual_equation_scaling_mode != jmi_residual_scaling_none)
     1792            kin_char_log(solver, 's');
    18041793        jmi_update_f_scale(block);
    18051794        jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block);
     
    21342123        jmi_log_node_t node = jmi_log_enter_fmt(block->log, logInfo, "AggressiveResidualScalingUpdate", "Updating f_scale aggressively");
    21352124        N_VScale(1.0, block->f_scale, solver->work_vector);
    2136         kin_char_log(solver, 's');
     2125        kin_char_log(solver, 's');
    21372126        jmi_update_f_scale(block);
    21382127        jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block);
     
    22512240        } else if (solver->handling_of_singular_jacobian_flag == JMI_MINIMUM_NORM) {
    22522241            /*
    2253              *   DGELSS - compute the minimum norm solution to  a real
     2242             *   DGELSS - compute the minimum norm solution to  a real
    22542243             *   linear least squares problem
    22552244             *
     
    24592448    solver->work_vector = N_VNew_Serial(n);
    24602449    solver->work_vector2 = N_VNew_Serial(n);
     2450    solver->work_vector3 = N_VNewEmpty_Serial(n);
    24612451    solver->lapack_work = (realtype*)calloc(4*(n+1),sizeof(realtype));
    24622452    solver->lapack_iwork = (int *)calloc(n+2, sizeof(int));
     
    25782568    N_VDestroy_Serial(solver->work_vector);
    25792569    N_VDestroy_Serial(solver->work_vector2);
     2570    N_VDestroy_Serial(solver->work_vector3);
    25802571    free(solver->lapack_work);
    25812572    free(solver->lapack_iwork);
     
    28222813    /* update the scaling only once per time step */
    28232814    if(block->init || (block->options->rescale_each_step_flag && (curtime > block->scale_update_time)) || block->force_rescaling) {
    2824         if(block->options->residual_equation_scaling_mode != jmi_residual_scaling_none)
    2825             kin_char_log(solver, 's');
    2826         jmi_update_f_scale(block);
     2815        if(block->options->residual_equation_scaling_mode != jmi_residual_scaling_none)
     2816            kin_char_log(solver, 's');
     2817        jmi_update_f_scale(block);
    28272818        jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block);
    28282819    }
     
    28982889        }
    28992890        /* Update the scaling  */
    2900         kin_char_log(solver, 's');
     2891        kin_char_log(solver, 's');
    29012892        jmi_update_f_scale(block);
    29022893        jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block);
  • trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.h

    r8512 r9581  
    110110    N_Vector work_vector;           /**< \brief work vector for vector operations */
    111111    N_Vector work_vector2;           /**< \brief work vector for vector operations */
     112    N_Vector work_vector3;           /**< \brief work vector for vector operations */
    112113    realtype* lapack_work;         /**< \brief work vector for lapack */
    113114    int * lapack_iwork;            /**< \brief work vector for lapack */
Note: See TracChangeset for help on using the changeset viewer.