Changeset 9581
 Timestamp:
 Jan 20, 2017 4:09:03 PM (3 years ago)
 Location:
 trunk/RuntimeLibrary/src/jmi
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.c
r9293 r9581 313 313 } 314 314 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 335 315 int jmi_kin_setup_column_partition(jmi_block_solver_t * block) { 336 316 jmi_kinsol_solver_t* solver = (jmi_kinsol_solver_t*)block>solver; … … 338 318 int leftColumns = N; 339 319 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); 342 323 while(leftColumns > 0 && group <N) { 343 324 group++; 344 N_VConst _Serial(0.0, solver>work_vector2);325 N_VConst(0.0, solver>work_vector2); 345 326 for(i = 0; i<N; i++) { 346 327 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 349 333 solver>jac_compression_groups[NleftColumns]=group; 350 334 solver>jac_compression_group_index[NleftColumns]=i; … … 355 339 } 356 340 } 341 N_VDestroy(col_i); 342 357 343 if(leftColumns != 0) { 358 344 jmi_log_node(block>log, logWarning, "ColumnPartitioning", "Column partitioning error, dependency data may be corrupt."); … … 638 624 N_VSetArrayPointer(DENSE_COL(J, solver>jac_compression_group_index[k]), jthCol); 639 625 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 } 642 631 /* Restore original array pointer in tmp2 */ 643 632 N_VSetArrayPointer(tmp2_data, tmp2); … … 1800 1789 1801 1790 if(block>force_rescaling) { 1802 1803 1791 if(block>options>residual_equation_scaling_mode != jmi_residual_scaling_none) 1792 kin_char_log(solver, 's'); 1804 1793 jmi_update_f_scale(block); 1805 1794 jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block); … … 2134 2123 jmi_log_node_t node = jmi_log_enter_fmt(block>log, logInfo, "AggressiveResidualScalingUpdate", "Updating f_scale aggressively"); 2135 2124 N_VScale(1.0, block>f_scale, solver>work_vector); 2136 2125 kin_char_log(solver, 's'); 2137 2126 jmi_update_f_scale(block); 2138 2127 jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block); … … 2251 2240 } else if (solver>handling_of_singular_jacobian_flag == JMI_MINIMUM_NORM) { 2252 2241 /* 2253 * DGELSS  compute the minimum norm solution to 2242 * DGELSS  compute the minimum norm solution to a real 2254 2243 * linear least squares problem 2255 2244 * … … 2459 2448 solver>work_vector = N_VNew_Serial(n); 2460 2449 solver>work_vector2 = N_VNew_Serial(n); 2450 solver>work_vector3 = N_VNewEmpty_Serial(n); 2461 2451 solver>lapack_work = (realtype*)calloc(4*(n+1),sizeof(realtype)); 2462 2452 solver>lapack_iwork = (int *)calloc(n+2, sizeof(int)); … … 2578 2568 N_VDestroy_Serial(solver>work_vector); 2579 2569 N_VDestroy_Serial(solver>work_vector2); 2570 N_VDestroy_Serial(solver>work_vector3); 2580 2571 free(solver>lapack_work); 2581 2572 free(solver>lapack_iwork); … … 2822 2813 /* update the scaling only once per time step */ 2823 2814 if(block>init  (block>options>rescale_each_step_flag && (curtime > block>scale_update_time))  block>force_rescaling) { 2824 2825 2826 2815 if(block>options>residual_equation_scaling_mode != jmi_residual_scaling_none) 2816 kin_char_log(solver, 's'); 2817 jmi_update_f_scale(block); 2827 2818 jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block); 2828 2819 } … … 2898 2889 } 2899 2890 /* Update the scaling */ 2900 2891 kin_char_log(solver, 's'); 2901 2892 jmi_update_f_scale(block); 2902 2893 jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block); 
trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.h
r8512 r9581 110 110 N_Vector work_vector; /**< \brief work vector for vector operations */ 111 111 N_Vector work_vector2; /**< \brief work vector for vector operations */ 112 N_Vector work_vector3; /**< \brief work vector for vector operations */ 112 113 realtype* lapack_work; /**< \brief work vector for lapack */ 113 114 int * lapack_iwork; /**< \brief work vector for lapack */
Note: See TracChangeset
for help on using the changeset viewer.