Ticket #3879: linear_solution_check.patch
File linear_solution_check.patch, 7.6 KB (added by , 5 years ago) 


RuntimeLibrary/src/jmi/jmi_linear_solver.h
33 33 extern void dgels_(char* TRANS, int* M, int* N, int* NRHS,double* A,int* LDA, double* B,int* LDB,double* WORK,int* LWORK,int* INFO ); 34 34 extern int dgeequ_(int *m, int *n, double *a, int * lda, double *r__, double *c__, double *rowcnd, double 35 35 *colcnd, double *amax, int *info); 36 extern void dgesdd_(char* JOBZ, int* M, int* N, double* A, int* LDA, 37 double* S, double* U, int* LDU, double* VT, int* LDVT, 38 double* WORK, int* LWORK, int* IWORK, int* INFO); 36 39 37 40 extern int dlaqge_(int *m, int *n, double *a, int * lda, double *r__, double *c__, double *rowcnd, double 38 41 *colcnd, double *amax, char *equed); 42 extern void dgecon_(char *norm, int *n, double *a, int *lda, double *anorm, double *rcond, 43 double *work, int *iwork, int *info); 44 extern double dlange_(char *norm, int *m, int *n, double *a, int *lda, 45 double *work); 46 extern double dlamch_(char *cmach); 39 47 40 48 typedef struct jmi_linear_solver_t jmi_linear_solver_t; 41 49 … … 54 62 jmi_real_t* factorization; /**< \brief Matrix for storing the Jacobian factorization */ 55 63 jmi_real_t* jacobian; /**< \brief Matrix for storing the Jacobian */ 56 64 jmi_real_t* jacobian_temp; /**< \brief Matrix for storing the Jacobian */ 65 jmi_real_t* jacobian_temp_augmented; /**< \brief Matrix for storing the temporary augmented Jacobian */ 57 66 jmi_real_t* singular_values; /**< \brief Vector for the singular values of the Jacobian */ 58 67 double* rScale; /**< \brief Row scaling of the Jacobian matrix */ 59 68 double* cScale; /**< \brief Column scaling of the Jacobian matrix */ … … 63 72 int iwork; 64 73 double* zero_vector; 65 74 jmi_real_t* rwork; 75 int* iwork_array; 76 double* dgesdd_work; /**< \brief Work vector for desdd */ 77 int dgesdd_lwork; /**< \brief Work vector for desdd */ 78 int* dgesdd_iwork; /**< \brief Work vector for desdd */ 66 79 }; 67 80 68 81 #endif /* _JMI_LINEAR_SOLVER_H */ 
RuntimeLibrary/src/jmi/jmi_linear_solver.c
23 23 #include "jmi_block_solver_impl.h" 24 24 #include "jmi_log.h" 25 25 26 27 /* Estimate condition number utilizing dgecon from LAPACK*/ 28 static double jmi_calculate_jacobian_condition_number(jmi_block_solver_t * block) { 29 jmi_linear_solver_t* solver = block>solver; 30 char norm = 'I'; 31 int N = block>n; 32 double jnorm = 1.0; 33 double rcond = 1.0; 34 int info; 35 36 /* Compute infinity norm of J to be used with dgecon */ 37 jnorm = dlange_(&norm, &N, &N, solver>jacobian, &N, solver>rwork); 38 39 /* Compute reciprocal condition number */ 40 dgecon_(&norm, &N, solver>factorization, &N, &jnorm, &rcond, solver>rwork, solver>iwork_array,&info); 41 42 return 1.0/rcond; 43 } 44 45 26 46 int jmi_linear_solver_new(jmi_linear_solver_t** solver_ptr, jmi_block_solver_t* block) { 27 47 jmi_linear_solver_t* solver= (jmi_linear_solver_t*)calloc(1,sizeof(jmi_linear_solver_t)); 28 48 /* jmi_t* jmi = block>jmi; … … 38 58 solver>factorization = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t)); 39 59 solver>jacobian = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t)); 40 60 solver>jacobian_temp = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t)); 61 solver>jacobian_temp_augmented = (jmi_real_t*)calloc(n_x*(n_x+1),sizeof(jmi_real_t)); 41 62 solver>singular_values = (jmi_real_t*)calloc(n_x,sizeof(jmi_real_t)); 42 63 solver>rScale = (double*)calloc(n_x,sizeof(double)); 43 64 solver>cScale = (double*)calloc(n_x,sizeof(double)); … … 46 67 47 68 solver>singular_jacobian = 0; 48 69 solver>iwork = 5*n_x; 70 solver>iwork_array = (int*)calloc(solver>iwork,sizeof(int)); 49 71 solver>rwork = (double*)calloc(solver>iwork,sizeof(double)); 50 72 solver>zero_vector = (double*)calloc(n_x,sizeof(double)); 73 74 solver>dgesdd_lwork = 10*n_x; 75 solver>dgesdd_work = (double*)calloc(solver>dgesdd_lwork,sizeof(double)); 76 solver>dgesdd_iwork = (int*)calloc(8*n_x,sizeof(int)); 51 77 52 78 for (i=0; i<n_x; i++) { 53 79 solver>zero_vector[i] = 0.0; … … 164 190 } 165 191 166 192 if((block>callbacks>log_options.log_level >= 5)) { 193 /* double rcond = jmi_calculate_jacobian_condition_number(block); */ 167 194 jmi_log_reals(block>log, destnode, logInfo, "initial_guess", block>initial, block>n); 168 jmi_log_reals(block>log, destnode, logInfo, "b", block>res, block>n); 195 jmi_log_reals(block>log, destnode, logInfo, "b", block>res, block>n); 196 /* jmi_log_reals(block>log, destnode, logInfo, "condition_number", &rcond, 1); */ 169 197 jmi_log_leave(block>log, destnode); 170 198 } 171 199 … … 174 202 i = 1; /* One rhs to solve for */ 175 203 176 204 if (solver>singular_jacobian == 1){ 205 char jobz = 'N', ceps = 'P', csfmin = 'S'; 206 int j,rank_aug = 0, m_x = n_x + 1; 207 double eps = dlamch_(&ceps); 208 double sfmin = dlamch_(&csfmin); 209 double threshold; 210 211 memcpy(solver>jacobian_temp_augmented, solver>jacobian, n_x*n_x*sizeof(jmi_real_t)); 212 memcpy(&solver>jacobian_temp_augmented[n_x*n_x], block>res, n_x*sizeof(jmi_real_t)); 177 213 /* 178 214 * DGELSS  compute the minimum norm solution to a real 179 215 * linear least squares problem … … 191 227 return 1; 192 228 } 193 229 230 /* Verify solution */ 231 /* 232 * SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, 233 * LWORK, IWORK, INFO ) 234 */ 235 dgesdd_(&jobz, &n_x, &m_x, solver>jacobian_temp_augmented, &n_x, 236 solver>singular_values, NULL, &n_x, NULL, &n_x, 237 solver>dgesdd_work, &solver>dgesdd_lwork, solver>dgesdd_iwork, &info); 238 239 if(info != 0) { 240 jmi_log_node(block>log, logError, "Error", "DGESDD failed to compute the SVD of the augmented system [A,b] in <block: %d> with error code <error: %s>", block>label, info); 241 return 1; 242 } 243 244 threshold = (eps*solver>singular_values[0] > sfmin)? eps*solver>singular_values[0]: sfmin; 245 for (j = 0; j < n_x; j++) 246 { 247 if (solver>singular_values[j] > threshold) { 248 rank_aug += 1; 249 } else { 250 break; 251 } 252 } 253 254 if (rank_aug != rank) { 255 jmi_log_node(block>log, logError, "Error", "Rank test failed, no solution exists in <block: %d>. Jacobian <rank: %d>, Augmented Jacobian <rank: %d>", block>label, rank, rank_aug); 256 return 1; 257 } 258 194 259 if(block>callbacks>log_options.log_level >= 5){ 195 260 jmi_log_node(block>log, logInfo, "Info", "Successfully calculated the minimum norm solution to the linear system in <block: %s>", block>label); 196 261 } … … 236 301 free(solver>jacobian); 237 302 free(solver>singular_values); 238 303 free(solver>jacobian_temp); 304 free(solver>jacobian_temp_augmented); 239 305 free(solver>rScale); 240 306 free(solver>cScale); 241 307 free(solver>rwork); 242 308 free(solver>zero_vector); 309 free(solver>iwork_array); 310 free(solver>dgesdd_work); 311 free(solver>dgesdd_iwork); 243 312 free(solver); 244 313 block>solver = 0; 245 314 }