Changeset 5726


Ignore:
Timestamp:
Jan 2, 2014 10:44:16 AM (6 years ago)
Author:
jakesson
Message:

Fixed computation of condition number in regularization check, see #3275.

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

Legend:

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

    r5720 r5726  
    727727
    728728/* Estimate condition number utilizing dgecon from LAPACK*/
    729 static realtype jmi_calculate_condition_number(jmi_block_solver_t * block, realtype* A) {
     729static realtype jmi_calculate_jacobian_condition_number(jmi_block_solver_t * block) {
    730730    jmi_kinsol_solver_t* solver = block->solver;
    731731    char norm = 'I';
    732732    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;
    734735    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?
    737752    if(Jcond < 0) Jcond = -Jcond;
    738753    if(Jcond == 0.0) Jcond = 1e-30;
    739     return 1.0/Jcond;
     754    */
     755    return 1.0/J_recip_cond;
    740756}
    741757
     
    10411057   
    10421058    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);
    10441060
    10451061        jmi_log_node(block->log, logInfo, "Regularization",
  • trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.h

    r5628 r5726  
    101101extern void dgecon_(char *norm, int *n, double *a, int *lda, double *anorm, double *rcond,
    102102             double *work, int *iwork, int *info);
    103 
     103extern double dlange_(char *norm, int *m, int *n, double *a, int *lda,
     104             double *work);
    104105extern int dgeequ_(int *m, int *n, double *a, int *
    105106    lda, double *r__, double *c__, double *rowcnd, double
Note: See TracChangeset for help on using the changeset viewer.