Changeset 5726
 Timestamp:
 Jan 2, 2014 10:44:16 AM (6 years ago)
 Location:
 trunk/RuntimeLibrary/src/jmi
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.c
r5720 r5726 727 727 728 728 /* Estimate condition number utilizing dgecon from LAPACK*/ 729 static realtype jmi_calculate_ condition_number(jmi_block_solver_t * block, realtype* A) {729 static realtype jmi_calculate_jacobian_condition_number(jmi_block_solver_t * block) { 730 730 jmi_kinsol_solver_t* solver = block>solver; 731 731 char norm = 'I'; 732 732 int N = block>n; 733 double Jnorm = 1.0, Jcond = 1.0; 733 double J_norm = 1.0; 734 double J_recip_cond = 1.0; 734 735 int info; 735 736 dgecon_(&norm, &N, A, &N, &Jnorm, &Jcond, solver>lapack_work, solver>lapack_iwork,&info); 736 737 /* Copy Jacobian to factorization matrix */ 738 DenseCopy(solver>J, solver>J_LU); 739 /* Perform LU factorization to be used with dgecon */ 740 dgetrf_(&N, &N, solver>J_LU>data, &N, solver>lapack_ipiv, &info); 741 if (info != 0 ) { 742 /* If matrix i singular, return something very large to be evaluated*/ 743 return 1e100; 744 } 745 746 /* Compute infinity norm of J to be used with dgecon */ 747 J_norm = dlange_(&norm, &N, &N, solver>J>data, &N, solver>lapack_work); 748 749 /* Compute reciprocal condition number */ 750 dgecon_(&norm, &N, solver>J_LU>data, &N, &J_norm, &J_recip_cond, solver>lapack_work, solver>lapack_iwork,&info); 751 /* To be evaluated  why is this needed? Error handling due to J being used instead of J_LU? 737 752 if(Jcond < 0) Jcond = Jcond; 738 753 if(Jcond == 0.0) Jcond = 1e30; 739 return 1.0/Jcond; 754 */ 755 return 1.0/J_recip_cond; 740 756 } 741 757 … … 1041 1057 1042 1058 if (solver>using_max_min_scaling_flag) { 1043 realtype cond = jmi_calculate_ condition_number(block, solver>J>data);1059 realtype cond = jmi_calculate_jacobian_condition_number(block); 1044 1060 1045 1061 jmi_log_node(block>log, logInfo, "Regularization", 
trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.h
r5628 r5726 101 101 extern void dgecon_(char *norm, int *n, double *a, int *lda, double *anorm, double *rcond, 102 102 double *work, int *iwork, int *info); 103 103 extern double dlange_(char *norm, int *m, int *n, double *a, int *lda, 104 double *work); 104 105 extern int dgeequ_(int *m, int *n, double *a, int * 105 106 lda, double *r__, double *c__, double *rowcnd, double
Note: See TracChangeset
for help on using the changeset viewer.