Changeset 8512


Ignore:
Timestamp:
Mar 1, 2016 9:15:21 AM (4 years ago)
Author:
aramle
Message:

#3879 Added a check for a valid minimum norm solution in case of a singular Jacobian.
#4475 In case we have an entire row in the Jacobian equal to 0, use heuristic residual scaling instead of the maximum.

Location:
trunk
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/Python/src/tests_jmodelica/files/Modelica/Singular.mo

    r8065 r8512  
    205205        (sin(x*y*z)+z) = 0;
    206206    end NonLinear5;
     207   
     208    model NoMinimumNormSolution
     209        Real x,y,z,v;
     210        parameter Real a11 = 1;
     211        parameter Real a12 = 2;
     212        parameter Real a13 = 3;
     213        parameter Real a21 = 1;
     214        parameter Real a22 = 1;
     215        parameter Real a23 = 1;
     216        Real a31(start=0);
     217        Real a32(start=0);
     218        Real a33(start=0);
     219        parameter Real b[3] = {1,2,3};
     220    equation
     221        a33=a31+a32;
     222        a33=a31^2+a32^2;
     223        a33=a31-a32;
     224        a11*x+a12*y+a13*z = b[1];
     225        a21*x+a22*y+a23*z = b[2];
     226        a31*x+a32*y+a33*z = b[3];
     227        der(v) = time;
     228    end NoMinimumNormSolution;
    207229
    208230end Singular;
  • trunk/Python/src/tests_jmodelica/simulation/test_assimulo_interface_fmi.py

    r8507 r8512  
    621621        compile_fmu("Singular.NonLinear4", file_name)
    622622        compile_fmu("Singular.NonLinear5", file_name)
     623        compile_fmu("Singular.NoMinimumNormSolution", file_name)
    623624   
    624625    @testattr(stddist = True)
     
    749750        nose.tools.assert_almost_equal(res["z"][0] ,0.000000000)
    750751        nose.tools.assert_almost_equal(res["z"][-1] ,-1.000000000)
     752   
     753    @testattr(stddist = True)
     754    def test_no_valid_minimum_norm_sol(self):
     755        model = load_fmu("NoMinimumNormSolution.fmu", log_level=3)
     756        model.set("_log_level", 3)
     757        model.set_log_level(log_level)
     758        nose.tools.assert_raises(FMUException, model.initialize)
    751759
    752760class Test_FMI_ODE_CS_2:
  • trunk/RuntimeLibrary/src/jmi/jmi_block_solver.c

    r8507 r8512  
    2828#include <string.h>
    2929#include <math.h>
     30
     31#include <sundials/sundials_math.h>
     32#include <sundials/sundials_direct.h>
     33#include <nvector/nvector_serial.h>
     34#include <kinsol/kinsol_direct.h>
     35#include <kinsol/kinsol_impl.h>
     36#include <sundials/sundials_dense.h>
    3037
    3138#include "jmi_log.h"
     
    5663                           jmi_block_solver_options_t* options,
    5764                           void* problem_data){
     65    int i;
    5866    jmi_block_solver_t* block_solver = (jmi_block_solver_t*)calloc(1, sizeof(jmi_block_solver_t));
    5967    if(!block_solver) {
     
    7482
    7583    block_solver->dx=(jmi_real_t*)calloc(n,sizeof(jmi_real_t));;                /**< \brief Work vector for the seed vector */
     84
     85    block_solver->f_scale = N_VNew_Serial(n);
     86    block_solver->scale_update_time = -1.0;
     87    if(n>0) {
     88        block_solver->J = NewDenseMat(n ,n);
     89        block_solver->J_scale = NewDenseMat(n ,n);
     90    }
    7691
    7792    block_solver->res = (jmi_real_t*)calloc(n,sizeof(jmi_real_t));
     
    8499    block_solver->max = (jmi_real_t*)calloc(n,sizeof(jmi_real_t));
    85100    block_solver->nominal = (jmi_real_t*)calloc(n,sizeof(jmi_real_t));
     101    block_solver->residual_nominal = (jmi_real_t*)calloc(n+1,sizeof(jmi_real_t));
     102    block_solver->residual_heuristic_nominal = (jmi_real_t*)calloc(n+1,sizeof(jmi_real_t));
    86103    block_solver->initial = (jmi_real_t*)calloc(n,sizeof(jmi_real_t));
    87104    block_solver->start_set = (jmi_real_t*)calloc(n,sizeof(jmi_real_t));
     
    97114
    98115    block_solver->cur_time = 0;
     116    block_solver->using_max_min_scaling_flag = 0; /* Not using max/min scaling */
     117
     118        /* Initial equation scaling is 1.0 */
     119    N_VConst_Serial(1.0,block_solver->f_scale);
     120
     121    /* Make sure that residual_nominals and heuristic nominals always have values */
     122    for(i=0; i<n; i++) {
     123        block_solver->residual_nominal[i] = 1.0;
     124        block_solver->residual_heuristic_nominal[i] = 1.0;
     125    }
     126
    99127    switch(options->solver) {
    100128        case JMI_KINSOL_SOLVER: {
     
    196224    free(block_solver->dx);
    197225
     226    N_VDestroy_Serial(block_solver->f_scale);
     227    if(block_solver->n > 0) {
     228        DestroyMat(block_solver->J);
     229        DestroyMat(block_solver->J_scale);
     230    }
     231
    198232    free(block_solver->res);
    199233    free(block_solver->dres);
     
    205239    free(block_solver->max);
    206240    free(block_solver->nominal);
     241    free(block_solver->residual_nominal);
     242    free(block_solver->residual_heuristic_nominal);
    207243    free(block_solver->initial);
    208244    free(block_solver->value_references);
     
    297333        block_solver->F(block_solver->problem_data,real_vrs,block_solver->res,JMI_BLOCK_VALUE_REFERENCE);
    298334        block_solver->F(block_solver->problem_data,block_solver->start_set,block_solver->res,JMI_BLOCK_START_SET);
     335        jmi_setup_f_residual_scaling(block_solver);
    299336
    300337        for (i=0;i<block_solver->n;i++) {
     
    710747}
    711748
     749/*
     750    Compute appropriate equation scaling and function tolerance based on Jacobian J,
     751    nominal values (block->nominal) and current point (block->x).
     752    Store result in block->f_scale.
     753*/
     754void jmi_update_f_scale(jmi_block_solver_t *block) {
     755    int i, N = block->n;
     756    realtype curtime = block->cur_time;
     757    realtype* scale_ptr = 0;
     758    jmi_block_solver_options_t* bsop = block->options;
     759    int use_scaling_flag = bsop->residual_equation_scaling_mode;
     760    realtype tol = bsop->res_tol;
     761
     762    block->scale_update_time = curtime;
     763    block->force_rescaling = 0;
     764   
     765    if(bsop->residual_equation_scaling_mode != jmi_residual_scaling_none) {
     766        /* Zero out the scales initially if we're modify this. */
     767        N_VConst_Serial(0,block->f_scale);
     768    }
     769
     770    scale_ptr = N_VGetArrayPointer(block->f_scale);
     771
     772    /* Form scaled Jacobian as needed for automatic scaling and condition number checking*/
     773    if (bsop->residual_equation_scaling_mode == jmi_residual_scaling_auto ||
     774        bsop->residual_equation_scaling_mode == jmi_residual_scaling_aggressive_auto ||
     775        bsop->residual_equation_scaling_mode == jmi_residual_scaling_full_jacobian_auto ||
     776        bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid ||
     777        bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual)
     778    {
     779        for (i = 0; i < N; i++) {
     780            int j;
     781            /* column scaling is formed by max(abs(nominal), abs(actual_value)) */
     782            realtype xscale = RAbs(block->nominal[i]);
     783            realtype x = RAbs(block->x[i]);
     784            realtype* col_ptr;
     785            realtype* scaled_col_ptr;
     786            if(x < xscale) x = xscale;
     787            if(x < tol) x = tol;
     788            col_ptr = DENSE_COL(block->J, i);
     789            scaled_col_ptr = DENSE_COL(block->J_scale, i);
     790
     791            /* row scaling is product of Jac entry and column scaling */
     792            for (j = 0; j < N; j++) {
     793                realtype dres = col_ptr[j];
     794                realtype fscale;
     795                fscale = dres * x;
     796                scaled_col_ptr[j] = fscale;
     797                scale_ptr[j] = MAX(scale_ptr[j], RAbs(fscale));
     798            }
     799        }
     800    }
     801
     802    /* put in manual non-zero scales in hybrid & manual cases */
     803    if(    (bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual)
     804        || (bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid))    {
     805            for(i = 0; i < N; i++) {
     806                if(block->residual_nominal[i] != 0) {
     807                    scale_ptr[i] = block->residual_nominal[i];
     808                }
     809            }
     810    }
     811
     812    if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_auto
     813        || bsop->residual_equation_scaling_mode == jmi_residual_scaling_aggressive_auto
     814        || bsop->residual_equation_scaling_mode == jmi_residual_scaling_full_jacobian_auto
     815        || bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid) {
     816        block->using_max_min_scaling_flag = 0; /* NOT using max/min scaling */
     817        /* check that scaling factors have reasonable magnitude */
     818        for(i = 0; i < N; i++) {
     819            if(scale_ptr[i] == 0.0) {
     820                jmi_log_node(block->log, logInfo, "JacobianZeroRow",
     821                    "The derivatives for <residual: %d> are all equal 0. Using heuristic residual nominals for scaling in <block: %s>", i, block->label);
     822                scale_ptr[i] = block->residual_heuristic_nominal[i];
     823            } else if(scale_ptr[i] < 1/bsop->max_residual_scaling_factor) {
     824                int j, maxInJacIndex = 0;
     825                realtype **jac = block->J_scale->cols;
     826                realtype maxInJac = RAbs(jac[0][i]);
     827                block->using_max_min_scaling_flag = 1; /* Using maximum scaling, singular Jacobian? */
     828                for(j = 0; j < N; j++) {
     829                    if(RAbs(maxInJac) < RAbs(jac[j][i])) {
     830                        maxInJac = jac[j][i];
     831                        maxInJacIndex = j;
     832                    }
     833                }
     834               
     835                scale_ptr[i] = bsop->max_residual_scaling_factor;
     836                jmi_log_node(block->log, logWarning, "MaxScalingUsed", "Poor scaling in <block: %s>. "
     837                    "Consider changes in the model. Partial derivative of <equation: %I> with respect to <Iter: #r%d#> is <dRes_dIter: %g> (<dRes_dIter_scaled: %g>).",
     838                    block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(block->J, i, maxInJacIndex) ,maxInJac);
     839            }
     840            else if(scale_ptr[i] > 1/bsop->min_residual_scaling_factor) {
     841                int j, maxInJacIndex = 0;
     842                realtype **jac = block->J_scale->cols;
     843                realtype maxInJac = RAbs(jac[0][i]);
     844                /* Likely not a problem: solver->using_max_min_scaling_flag = 1; -- Using minimum scaling */
     845                for(j = 0; j < N; j++) {
     846                    if(RAbs(maxInJac) < RAbs(jac[j][i])) {
     847                        maxInJac = jac[j][i];
     848                        maxInJacIndex = j;
     849                    }
     850                }
     851               
     852                if (scale_ptr[i] > block->residual_heuristic_nominal[i]/bsop->min_residual_scaling_factor) {
     853                    scale_ptr[i] = bsop->min_residual_scaling_factor/block->residual_heuristic_nominal[i];
     854                    jmi_log_node(block->log, logWarning, "MinScalingUsed", "Poor scaling in <block: %s>. "
     855                        "Consider changes in the model. Partial derivative of <equation: %I> with respect to <Iter: #r%d#> is <dRes_dIter: %g> (<dRes_dIter_scaled: %g>).",
     856                        block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(block->J, i, maxInJacIndex) ,maxInJac);
     857                } else {
     858                    scale_ptr[i] = 1/scale_ptr[i];
     859                    jmi_log_node(block->log, logWarning, "MinScalingExceeded", "Poor scaling in <block: %s> but will allow it based on structure of block. "
     860                        "Consider changes in the model. Partial derivative of <equation: %I> with respect to <Iter: #r%d#> is <dRes_dIter: %g> (<dRes_dIter_scaled: %g>).",
     861                        block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(block->J, i, maxInJacIndex) ,maxInJac);
     862                }
     863            }
     864            else
     865                scale_ptr[i] = 1/scale_ptr[i];
     866        }
     867
     868        if (block->callbacks->log_options.log_level >= 4) {
     869            jmi_log_node_t outer = jmi_log_enter_fmt(block->log, logInfo, "ResidualScalingUpdated", "<block:%s>", block->label);
     870            if (block->callbacks->log_options.log_level >= 5) {
     871                jmi_log_node_t inner = jmi_log_enter_vector_(block->log, outer, logInfo, "scaling");
     872                realtype* res = scale_ptr;
     873                for (i=0;i<N;i++) jmi_log_real_(block->log, 1/res[i]);
     874                jmi_log_leave(block->log, inner);
     875            }
     876            jmi_log_leave(block->log, outer);
     877        }
     878    } else if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) {
     879        for(i = 0; i < N; i++) {
     880            scale_ptr[i] = 1/scale_ptr[i];
     881        }
     882        /* Scaling is not updated */
     883        if (block->callbacks->log_options.log_level >= 4) {
     884            jmi_log_node_t outer = jmi_log_enter_fmt(block->log, logInfo, "ResidualScalingUpdated", "<block:%s>", block->label);
     885            if (block->callbacks->log_options.log_level >= 5) {
     886                jmi_log_node_t inner = jmi_log_enter_vector_(block->log, outer, logInfo, "scaling");
     887                realtype* res = scale_ptr;
     888                for (i=0;i<N;i++) jmi_log_real_(block->log, 1/res[i]);
     889                jmi_log_leave(block->log, inner);
     890            }
     891            jmi_log_leave(block->log, outer);
     892        }
     893    }     
     894    return;
     895}
     896
     897void jmi_setup_f_residual_scaling(jmi_block_solver_t *block) {
     898    jmi_block_solver_options_t* bsop = block->options;
     899    int i, N = block->n;
     900    realtype* dummy = 0;
     901   
     902    /* Read manual scaling from annotations and propagated nominal values for
     903     * residuals and put them in residual_nominal & scale_ptr*/
     904    if (bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual ||
     905        bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid)
     906    {
     907        block->F(block->problem_data,dummy, block->residual_nominal, JMI_BLOCK_EQUATION_NOMINAL);
     908        for (i = 0; i < N; i++) {
     909            if(block->residual_nominal[i] != 0.0) {
     910                if(block->residual_nominal[i] < 1/bsop->max_residual_scaling_factor) {
     911                    block->residual_nominal[i] = 1/bsop->max_residual_scaling_factor;
     912                    jmi_log_node(block->log, logWarning, "MaxScalingUsed", "Using maximum scaling factor in <block: %s>, "
     913                        "<equation: %I> since specified residual nominal is too small.", block->label, i);
     914                }
     915                else if(block->residual_nominal[i] > 1/bsop->min_residual_scaling_factor) {
     916                    block->residual_nominal[i] = 1/bsop->min_residual_scaling_factor;
     917                    jmi_log_node(block->log, logWarning, "MinScalingUsed", "Using minimal scaling factor in <block: %s>, "
     918                        "<equation: %I> since specified residual nominal is too big.", block->label, i);
     919                }
     920            }
     921            else if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) {
     922                block->residual_nominal[i] = 1.0;
     923            }
     924        }
     925    }
     926   
     927    block->F(block->problem_data,dummy, block->residual_heuristic_nominal, JMI_BLOCK_EQUATION_NOMINAL_AUTO);
     928    block->F(block->problem_data,dummy, block->residual_heuristic_nominal, JMI_BLOCK_EQUATION_NOMINAL);
     929    for (i = 0; i < N; i++) {
     930        if (block->residual_heuristic_nominal[i] == 0.0) {
     931            /* TODO: This information could be used to see that the system is structually singular */
     932            block->residual_heuristic_nominal[i] = 1.0;
     933        }
     934    }
     935}
     936
     937double* jmi_solver_get_f_scales(jmi_block_solver_t* block) {
     938    return N_VGetArrayPointer(block->f_scale);
     939}
     940
     941double calculate_scaled_residual_norm(double* residual, double *f_scale, int n) {
     942    double scaledRes, norm = 0;
     943    int i;
     944    for(i=0; i<n; i++) {
     945        scaledRes = residual[i]*f_scale[i];
     946        norm = norm>= RAbs(scaledRes)?norm:RAbs(scaledRes);
     947    }
     948    return norm;
     949}
     950
  • trunk/RuntimeLibrary/src/jmi/jmi_block_solver.h

    r8507 r8512  
    2828#include "jmi_log.h"
    2929#include <time.h>
    30 
     30#include <sundials/sundials_math.h>
     31#include <sundials/sundials_direct.h>
     32#include <nvector/nvector_serial.h>
     33#include <kinsol/kinsol_direct.h>
     34#include <kinsol/kinsol_impl.h>
     35#include <sundials/sundials_dense.h>
    3136
    3237/** \brief Evaluation modes for the residual function.*/
     
    350355void jmi_block_solver_init_default_options(jmi_block_solver_options_t* op);
    351356
     357/** \brief Update function scaling based on Jacobian information */
     358void jmi_update_f_scale(jmi_block_solver_t *block);
     359
     360/**< \brief Retrieve residual scales used in solver */
     361double* jmi_solver_get_f_scales(jmi_block_solver_t* block);
     362
     363/**< \brief Setup residual scaling */
     364void jmi_setup_f_residual_scaling(jmi_block_solver_t *block);
     365
     366/** \brief Calculate scaled residual max norm */
     367double calculate_scaled_residual_norm(double* residual, double *f_scale, int n);
     368
    352369/** \brief Check and log illegal iv inputs */
    353370int jmi_check_and_log_illegal_iv_input(jmi_block_solver_t* block, double* ivs, int N);
  • trunk/RuntimeLibrary/src/jmi/jmi_block_solver_impl.h

    r8313 r8512  
    3737    jmi_string_t label;
    3838
     39    N_Vector f_scale;          /**< \brief Work vector for scaling of f */
     40    realtype scale_update_time; /**< \brief The last time when f scale was updated */
    3941    int n;                         /**< \brief The number of iteration variables */
    4042    jmi_real_t* x;                 /**< \brief Work vector for the real iteration variables */
    4143    jmi_real_t* last_accepted_x;   /**< \brief Work vector for the real iteration variables holding the last accepted vales by the integrator */
     44    DlsMat J;                       /**< \brief The Jacobian matrix  */
     45    DlsMat J_scale;                 /**< \brief Jacobian matrix scaled with xnorm for used for fnorm calculation */
     46    int using_max_min_scaling_flag; /**< \brief A flag indicating if either the maximum scaling is used of the minimum */
    4247
    4348    jmi_real_t* dx;                /**< \brief Work vector for the seed vector */
     
    6671    jmi_real_t* max;               /**< \brief Max values for iteration variables */
    6772    jmi_real_t* nominal;           /**< \brief Nominal values for iteration variables */
     73    jmi_real_t* residual_nominal;   /**< \brief Nominals values for residual variables */
     74    jmi_real_t* residual_heuristic_nominal;   /**< \brief Heuristic nominals values for residual variables */
    6875    jmi_real_t* initial;           /**< \brief Initial values for iteration variables */
    6976    jmi_real_t* start_set;         /**< \brief If the start value is specified for the iteration variables */
    70    
     77
    7178    int jacobian_variability;      /**< \brief Variability of Jacobian coefficients: JMI_CONSTANT_VARIABILITY
    7279                                         JMI_PARAMETER_VARIABILITY, JMI_DISCRETE_VARIABILITY, JMI_CONTINUOUS_VARIABILITY */
  • trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.c

    r8435 r8512  
    4848static int jmi_kin_lsetup(struct KINMemRec * kin_mem);
    4949
    50 /* Update residual scales based on jacobian information */
    51 static void jmi_update_f_scale(jmi_block_solver_t *block);
     50/* Only call directly after jmi_update_f_scale */
     51static void jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(jmi_block_solver_t *block);
     52
     53static realtype jmi_calculate_jacobian_condition_number(jmi_block_solver_t * block);
    5254
    5355/* Calculate matrix vector product vv = M*v or vv = M^T*v */
     
    174176    /* Test if output is OK (no -1.#IND) */
    175177    n = NV_LENGTH_S(ff);
    176     ret = jmi_check_and_log_illegal_residual_output(block, f, y, solver->residual_heuristic_nominal ,n);
     178    ret = jmi_check_and_log_illegal_residual_output(block, f, y, block->residual_heuristic_nominal ,n);
    177179   
    178180    /* record information for Brent search */
     
    702704            int i, info, N = block->n;
    703705            realtype tol = solver->kin_stol;
    704             realtype *scale_ptr = N_VGetArrayPointer(solver->kin_f_scale);
     706            realtype *scale_ptr = N_VGetArrayPointer(block->f_scale);
    705707            if(block->options->residual_equation_scaling_mode != 0) {
    706708                for(i = 0; i < N; i++){
    707709                    int j;
    708                     realtype* scaled_col_ptr = DENSE_COL(solver->J_scale, i);
    709                     realtype* col_ptr = DENSE_COL(solver->J, i);
     710                    realtype* scaled_col_ptr = DENSE_COL(block->J_scale, i);
     711                    realtype* col_ptr = DENSE_COL(block->J, i);
    710712                    realtype xscale = RAbs(block->nominal[i]);
    711713                    realtype x = RAbs(block->x[i]);
     
    717719                }
    718720            } else {
    719                 DenseCopy(solver->J, solver->J_scale);
    720             }
    721             dgetrf_(  &N, &N, solver->J_scale->data, &N, solver->lapack_iwork, &info);
     721                DenseCopy(block->J, block->J_scale);
     722            }
     723            dgetrf_(  &N, &N, block->J_scale->data, &N, solver->lapack_iwork, &info);
    722724            if(info > 0) {
    723725                jmi_log_node(block->log, logWarning, "SingularJacobian",
     
    725727            }
    726728            else {
    727                 dgecon_(&norm, &N, solver->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
     729                dgecon_(&norm, &N, block->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
    728730
    729731                if(tol * Jcond < UNIT_ROUNDOFF) {
     
    809811    jmi_kinsol_solver_t* solver = block->solver;
    810812    struct KINMemRec* kin_mem = solver->kin_mem;
    811     realtype* residual_scaling_factors = N_VGetArrayPointer(solver->kin_f_scale);
     813    realtype* residual_scaling_factors = N_VGetArrayPointer(block->f_scale);
    812814    jmi_log_t *log = block->log;
    813815    clock_t t = jmi_block_solver_start_clock(block);
     
    10371039    jmi_log_fmt(block->log, node, logInfo, " <max_iter_no_jacobian: %d>", block->options->max_iter_no_jacobian);
    10381040    jmi_log_fmt(block->log, node, logInfo, " <active_bounds_mode: %d>", block->options->active_bounds_mode);
     1041    jmi_log_fmt(block->log, node, logInfo, " <start_from_last_integrator_step: %d>", block->options->start_from_last_integrator_step);
    10391042    jmi_log_leave(block->log, node);
    10401043
     
    14901493    int i,j,k;
    14911494    realtype **JTJ_c =  solver->JTJ->cols;
    1492     realtype **jac = solver->J->cols;
     1495    realtype **jac = block->J->cols;
    14931496    int N = block->n;
    14941497    realtype * uscale_data = N_VGetArrayPointer(solver->kin_y_scale);
    1495     realtype * fscale_data = N_VGetArrayPointer(solver->kin_f_scale);   
     1498    realtype * fscale_data = N_VGetArrayPointer(block->f_scale);   
    14961499
    14971500    for (i=0;i<N;i++) {
     
    15471550}
    15481551
     1552/* Call this function directly after jmi_update_f_scale */
     1553static void jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(jmi_block_solver_t *block) {
     1554
     1555    jmi_kinsol_solver_t* solver = block->solver;
     1556    int i, N = block->n;
     1557    jmi_block_solver_options_t* bsop = block->options;
     1558    int use_scaling_flag = bsop->residual_equation_scaling_mode;
     1559    realtype tol = solver->kin_stol;
     1560    realtype* scale_ptr = N_VGetArrayPointer(block->f_scale);
     1561
     1562    if (block->using_max_min_scaling_flag && !solver->J_is_singular_flag) {
     1563        realtype cond = jmi_calculate_jacobian_condition_number(block);
     1564
     1565        jmi_log_node(block->log, logInfo, "Regularization",
     1566            "Calculated condition number in <block: %s>. Regularizing if <cond: %E> is greater than <regtol: %E>", block->label, cond, solver->kin_reg_tol);
     1567        if (cond > solver->kin_reg_tol) {
     1568            if(N > 1 && solver->handling_of_singular_jacobian_flag == JMI_REGULARIZATION) {
     1569                int info;
     1570                jmi_kinsol_reg_matrix(block);
     1571                dgetrf_(  &block->n, &block->n, solver->JTJ->data, &block->n, solver->lapack_ipiv, &info);
     1572            }
     1573            solver->J_is_singular_flag = 1;
     1574        }
     1575    }
     1576
     1577
     1578    /* estimate condition number of the scaled jacobian
     1579    and scale function tolerance with it. */
     1580    if((N > 1) && bsop->check_jac_cond_flag){
     1581        realtype* scaled_col_ptr;
     1582        char norm = 'I';
     1583        double Jnorm = 1.0, Jcond = 1.0;
     1584        int info;
     1585        if(use_scaling_flag) {
     1586            for(i = 0; i < N; i++){
     1587                int j;
     1588                scaled_col_ptr = DENSE_COL(block->J_scale, i);
     1589                for(j = 0; j < N; j++){
     1590                    scaled_col_ptr[j] = scaled_col_ptr[j] * scale_ptr[j];
     1591                }
     1592            }
     1593        } else {
     1594            DenseCopy(block->J, block->J_scale);
     1595        }
     1596
     1597        dgetrf_(  &N, &N, block->J_scale->data, &N, solver->lapack_iwork, &info);
     1598        if(info > 0) {
     1599            jmi_log_node(block->log, logWarning, "SingularJacobian",
     1600                "Singular Jacobian detected when checking condition number in <block:%s>. Solver may fail to converge.", block->label);
     1601        }
     1602        else {
     1603            dgecon_(&norm, &N, block->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
     1604
     1605            if(tol * Jcond < UNIT_ROUNDOFF) {
     1606                jmi_log_node_t node = jmi_log_enter_fmt(block->log, logWarning, "IllConditionedJacobian",
     1607                    "<JacobianInverseConditionEstimate:%E> Solver may fail to converge in <block: %s>.", Jcond, block->label);
     1608                if (block->callbacks->log_options.log_level >= 4) {
     1609                    jmi_log_reals(block->log, node, logWarning, "ivs", N_VGetArrayPointer(solver->kin_y), block->n);
     1610                }
     1611                jmi_log_leave(block->log, node);
     1612
     1613            }
     1614            else {
     1615                jmi_log_node(block->log, logInfo, "JacobianCondition",
     1616                    "<JacobianInverseConditionEstimate:%E>", Jcond);
     1617            }
     1618        }
     1619    }
     1620}
    15491621
    15501622/* Estimate condition number utilizing dgecon from LAPACK*/
     
    15581630
    15591631    /* Copy Jacobian to factorization matrix */
    1560     DenseCopy(solver->J, solver->J_LU);
     1632    DenseCopy(block->J, solver->J_LU);
    15611633    /* Perform LU factorization to be used with dgecon */
    15621634    dgetrf_(&N, &N, solver->J_LU->data, &N, solver->lapack_ipiv, &info);
     
    15671639
    15681640    /* Compute infinity norm of J to be used with dgecon */
    1569     J_norm = dlange_(&norm, &N, &N, solver->J->data, &N, solver->lapack_work);
     1641    J_norm = dlange_(&norm, &N, &N, block->J->data, &N, solver->lapack_work);
    15701642
    15711643    /* Compute reciprocal condition number */
     
    15971669        return 0;
    15981670    }
    1599     SetToZero(solver->J);
     1671    SetToZero(block->J);
    16001672
    16011673    /* Evaluate Jacobian */
    1602     ret = kin_dF(N, solver->kin_y, kin_mem->kin_fval, solver->J, block, kin_mem->kin_vtemp1, kin_mem->kin_vtemp2);
     1674    ret = kin_dF(N, solver->kin_y, kin_mem->kin_fval, block->J, block, kin_mem->kin_vtemp1, kin_mem->kin_vtemp2);
    16031675    solver->updated_jacobian_flag = 1; /* The Jacobian is current */
    16041676   
     
    16071679        jmi_log_node(block->log, logWarning, "CompressedJacobianEvaluation", "Failed to evaluate Jacobian using compression. Will try full finite differences.");
    16081680        /* Try if it helps with full finite differences with possibility to step in other direction */
    1609         ret = kin_dF(N, solver->kin_y, kin_mem->kin_fval, solver->J, block, kin_mem->kin_vtemp1, kin_mem->kin_vtemp2);
     1681        ret = kin_dF(N, solver->kin_y, kin_mem->kin_fval, block->J, block, kin_mem->kin_vtemp1, kin_mem->kin_vtemp2);
    16101682        solver->has_compression_setup_flag = TRUE;
    16111683    }
     
    16261698    }
    16271699
    1628     if(block->force_rescaling)
     1700    if(block->force_rescaling) {
    16291701        jmi_update_f_scale(block);
     1702        jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block);
     1703    }
    16301704   
    16311705    return 0;
     
    16461720    denom = jmi_kinsol_calc_v1twwv2(kin_mem->kin_pp,kin_mem->kin_pp,solver->kin_y_scale);
    16471721    /* work_vector = Jac * step */
    1648     jmi_kinsol_calc_Mv(solver->J, 0, kin_mem->kin_pp, solver->work_vector);
     1722    jmi_kinsol_calc_Mv(block->J, 0, kin_mem->kin_pp, solver->work_vector);
    16491723
    16501724    /* work_vector = (ResidualDelta  - Jac * step) unless below update tolerance */
     
    16591733
    16601734    for(j=0; j < N; j++) {
    1661         realtype *jacCol_j = solver->J->cols[j];
     1735        realtype *jacCol_j = block->J->cols[j];
    16621736        double tempj = Ith(kin_mem->kin_pp, j)*Ith(solver->kin_y_scale,j)*Ith(solver->kin_y_scale,j);
    16631737        for(i = 0; i < N; i++) {
     
    16801754            jmi_log_fmt(block->log, node,logInfo, "Update denominator <step_norm: %g>", denom);
    16811755#endif
    1682             jmi_log_real_matrix(block->log, node, logInfo, "jacobian", solver->J->data, N, N);
     1756            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", block->J->data, N, N);
    16831757        }
    16841758        jmi_log_leave(block->log, node);
     
    16991773
    17001774    /* work_vector = Jac * step */
    1701     jmi_kinsol_calc_Mv(solver->J, 0, kin_mem->kin_pp, solver->work_vector);
     1775    jmi_kinsol_calc_Mv(block->J, 0, kin_mem->kin_pp, solver->work_vector);
    17021776
    17031777    /* work_vector = (ResidualDelta  - Jac * step) unless below update tolerance */
     
    17151789        double denom = 0;
    17161790        for(j=0; j < N; j++) {
    1717             denom += Ith(kin_mem->kin_pp, j)*Ith(kin_mem->kin_pp, j)*solver->J->cols[j][i]*solver->J->cols[j][i];
     1791            denom += Ith(kin_mem->kin_pp, j)*Ith(kin_mem->kin_pp, j)*block->J->cols[j][i]*block->J->cols[j][i];
    17181792        }
    17191793        Ith(solver->work_vector2,i) = denom;
     
    17211795
    17221796    for(j=0; j < N; j++) {
    1723         realtype *jacCol_j = solver->J->cols[j];
     1797        realtype *jacCol_j = block->J->cols[j];
    17241798        for(i = 0; i < N; i++) {
    17251799            if( Ith(solver->work_vector2,i) >= UNIT_ROUNDOFF) {
     
    17391813            jmi_log_reals(block->log, node,logInfo, "res_JacStep", N_VGetArrayPointer(solver->work_vector), N);
    17401814#endif
    1741             jmi_log_real_matrix(block->log, node, logInfo, "jacobian", solver->J->data, N, N);
     1815            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", block->J->data, N, N);
    17421816        }
    17431817        jmi_log_leave(block->log, node);
     
    17591833    /* work_vector2 = res_diff */
    17601834    N_VLinearSum(1, solver->last_residual, -1, b, solver->work_vector);
    1761     N_VProd(solver->work_vector, solver->kin_f_scale, solver->work_vector);
    1762     N_VProd(solver->work_vector, solver->kin_f_scale, solver->work_vector);
     1835    N_VProd(solver->work_vector, block->f_scale, solver->work_vector);
     1836    N_VProd(solver->work_vector, block->f_scale, solver->work_vector);
    17631837    /* work_vector2 = res_diff^T*f_scale^2*Jac */
    1764     jmi_kinsol_calc_Mv(solver->J, TRUE, solver->work_vector, solver->work_vector2);
     1838    jmi_kinsol_calc_Mv(block->J, TRUE, solver->work_vector, solver->work_vector2);
    17651839
    17661840    /* denom1 = step^T.*kin_y_scale^2*step */
    17671841    denom1 = jmi_kinsol_calc_v1twwv2(kin_mem->kin_pp,kin_mem->kin_pp,solver->kin_y_scale);
    17681842    /* work_vector = Jac * step */
    1769     jmi_kinsol_calc_Mv(solver->J, 0, kin_mem->kin_pp, solver->work_vector);
     1843    jmi_kinsol_calc_Mv(block->J, 0, kin_mem->kin_pp, solver->work_vector);
    17701844    /* denom2 = res_diff^T.*kin_f_scale^2*Jac*step */
    17711845    denom2 = 0;
    17721846    for(i = 0; i < N; i++) {
    1773         denom2 += -(Ith(b,i) - Ith(solver->last_residual,i))*Ith(solver->kin_f_scale,i)*Ith(solver->kin_f_scale,i)*Ith(solver->work_vector,i);
     1847        denom2 += -(Ith(b,i) - Ith(solver->last_residual,i))*Ith(block->f_scale,i)*Ith(block->f_scale,i)*Ith(solver->work_vector,i);
    17741848    }
    17751849   
    17761850    for(j=0; j < N; j++) {
    1777         realtype *jacCol_j = solver->J->cols[j];
     1851        realtype *jacCol_j = block->J->cols[j];
    17781852        double tempj1 = Ith(kin_mem->kin_pp,j)*Ith(solver->kin_y_scale,j)*Ith(solver->kin_y_scale,j)/denom1;
    17791853        double tempj2 = Ith(solver->work_vector2,j)/denom2;
     
    17971871            jmi_log_fmt(block->log, node,logInfo, "Update denominator <denom2: %g>", denom2);
    17981872#endif
    1799             jmi_log_real_matrix(block->log, node, logInfo, "jacobian", solver->J->data, N, N);
     1873            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", block->J->data, N, N);
    18001874        }
    18011875        jmi_log_leave(block->log, node);
     
    18151889    clock_t t = jmi_block_solver_start_clock(block);
    18161890     
    1817     DenseCopy(solver->J, solver->J_LU); /* make a copy of the Jacobian that will be used for LU factorization */
     1891    DenseCopy(block->J, solver->J_LU); /* make a copy of the Jacobian that will be used for LU factorization */
    18181892
    18191893    /* Equillibrate if corresponding option is set */
     
    18751949                             "Will try to find the minimum norm solution in <block: %s>", block->label);
    18761950                SetToZero(solver->J_sing);
    1877                 DenseCopy(solver->J, solver->J_sing);
     1951                DenseCopy(block->J, solver->J_sing);
    18781952            } else {
    18791953                /* Error */
     
    19502024        int i;
    19512025        jmi_log_node_t node = jmi_log_enter_fmt(block->log, logInfo, "AggressiveResidualScalingUpdate", "Updating f_scale aggressively");
    1952         N_VScale(1.0, solver->kin_f_scale, solver->work_vector);
     2026        N_VScale(1.0, block->f_scale, solver->work_vector);
    19532027        jmi_update_f_scale(block);
     2028        jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block);
    19542029        for(i=0; i<block->n; i++) {
    1955             Ith(solver->kin_f_scale, i) = Ith(solver->kin_f_scale, i) > Ith( solver->work_vector, i)? Ith( solver->work_vector, i):Ith(solver->kin_f_scale, i);
     2030            Ith(block->f_scale, i) = Ith(block->f_scale, i) > Ith( solver->work_vector, i)? Ith( solver->work_vector, i):Ith(block->f_scale, i);
    19562031        }
    19572032        if (block->callbacks->log_options.log_level >= 4) {
     
    19592034            if (block->callbacks->log_options.log_level >= 5) {
    19602035                jmi_log_node_t inner = jmi_log_enter_vector_(block->log, outer, logInfo, "scaling");
    1961                 realtype* res = N_VGetArrayPointer(solver->kin_f_scale);
     2036                realtype* res = N_VGetArrayPointer(block->f_scale);
    19622037                for (i=0;i<N;i++) jmi_log_real_(block->log, 1/res[i]);
    19632038                jmi_log_leave(block->log, inner);
     
    19932068    }
    19942069 
    1995     kin_mem->kin_sJpnorm = N_VWL2Norm(b,solver->kin_f_scale);
     2070    kin_mem->kin_sJpnorm = N_VWL2Norm(b,block->f_scale);
    19962071    solver->sJpnorm = kin_mem->kin_sJpnorm;
    1997     N_VProd(b, solver->kin_f_scale, x);
    1998     N_VProd(x, solver->kin_f_scale, x);
     2072    N_VProd(b, block->f_scale, x);
     2073    N_VProd(x, block->f_scale, x);
    19992074   
    20002075    kin_mem->kin_sfdotJp = N_VDotProd(kin_mem->kin_fval, x);
     
    20142089         */
    20152090
    2016         jmi_kinsol_calc_Mv(solver->J, TRUE, x, solver->gradient);
     2091        jmi_kinsol_calc_Mv(block->J, TRUE, x, solver->gradient);
    20172092    }
    20182093    if(solver->use_steepest_descent_flag) {
     
    20232098    }
    20242099    else if(solver->J_is_singular_flag) {
    2025         realtype** jac = solver->J->cols;
     2100        realtype** jac = block->J->cols;
    20262101        if (N == 1) {
    20272102            xd[0] = block->nominal[0] * 0.1 *((bd[0] > 0)?1:-1) * ((jac[0][0] > 0)?1:-1);
     
    20302105            /* solve the regularized problem */
    20312106            int i,j;
    2032             realtype * fscale_data = N_VGetArrayPointer(solver->kin_f_scale);
     2107            realtype * fscale_data = N_VGetArrayPointer(block->f_scale);
    20332108            realtype gnorm;/*gradient norm*/
    20342109
     
    21202195       
    21212196        if((block->callbacks->log_options.log_level >= 6)) {
    2122             jmi_log_real_matrix(block->log, node, logInfo, "jacobian", solver->J->data, N, N);
     2197            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", block->J->data, N, N);
    21232198            jmi_log_reals(block->log, node, logInfo, "rhs", xd, N);
    21242199        }
     
    22032278           
    22042279            /* recalculate sJpnorm  */
    2205             jmi_kinsol_calc_Mv(solver->J, FALSE, x, b);
    2206             sJpnorm = N_VWL2Norm(b,solver->kin_f_scale);
     2280            jmi_kinsol_calc_Mv(block->J, FALSE, x, b);
     2281            sJpnorm = N_VWL2Norm(b,block->f_scale);
    22072282            /* update sJpnorm and sfdotJp for Kinsol */
    22082283            kin_mem->kin_sJpnorm = sJpnorm;
     
    22192294    block->step_calc_time += jmi_block_solver_elapsed_time(block, t);
    22202295    return 0;
    2221 }
    2222 
    2223 double* jmi_kinsol_solver_get_f_scales(jmi_kinsol_solver_t* solver) {
    2224     return N_VGetArrayPointer(solver->kin_f_scale);
    2225 }
    2226 
    2227 static void jmi_setup_f_residual_scaling(jmi_block_solver_t *block) {
    2228     jmi_block_solver_options_t* bsop = block->options;
    2229     jmi_kinsol_solver_t* solver = block->solver;
    2230     int i, N = block->n;
    2231     realtype* dummy = 0;
    2232    
    2233     /* Read manual scaling from annotations and propagated nominal values for
    2234      * residuals and put them in residual_nominal & scale_ptr*/
    2235     if (bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual ||
    2236         bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid)
    2237     {
    2238         block->F(block->problem_data,dummy, solver->residual_nominal, JMI_BLOCK_EQUATION_NOMINAL);
    2239         for (i = 0; i < N; i++) {
    2240             if(solver->residual_nominal[i] != 0.0) {
    2241                 if(solver->residual_nominal[i] < 1/bsop->max_residual_scaling_factor) {
    2242                     solver->residual_nominal[i] = 1/bsop->max_residual_scaling_factor;
    2243                     jmi_log_node(block->log, logWarning, "MaxScalingUsed", "Using maximum scaling factor in <block: %s>, "
    2244                         "<equation: %I> since specified residual nominal is too small.", block->label, i);
    2245                 }
    2246                 else if(solver->residual_nominal[i] > 1/bsop->min_residual_scaling_factor) {
    2247                     solver->residual_nominal[i] = 1/bsop->min_residual_scaling_factor;
    2248                     jmi_log_node(block->log, logWarning, "MinScalingUsed", "Using minimal scaling factor in <block: %s>, "
    2249                         "<equation: %I> since specified residual nominal is too big.", block->label, i);
    2250                 }
    2251             }
    2252             else if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) {
    2253                 solver->residual_nominal[i] = 1.0;
    2254             }
    2255         }
    2256     }
    2257    
    2258     block->F(block->problem_data,dummy, solver->residual_heuristic_nominal, JMI_BLOCK_EQUATION_NOMINAL_AUTO);
    2259     block->F(block->problem_data,dummy, solver->residual_heuristic_nominal, JMI_BLOCK_EQUATION_NOMINAL);
    2260     for (i = 0; i < N; i++) {
    2261         if (solver->residual_heuristic_nominal[i] == 0.0) {
    2262             /* TODO: This information could be used to see that the system is structually singular */
    2263             solver->residual_heuristic_nominal[i] = 1.0;
    2264         }
    2265     }
    2266 }
    2267 
    2268 /*
    2269     Compute appropriate equation scaling and function tolerance based on Jacobian J,
    2270     nominal values (block->nominal) and current point (block->x).
    2271     Store result in solver->kin_f_scale.
    2272 */
    2273 static void jmi_update_f_scale(jmi_block_solver_t *block) {
    2274     jmi_kinsol_solver_t* solver = block->solver;
    2275     int i, N = block->n;
    2276     realtype tol = solver->kin_stol;
    2277     realtype curtime = block->cur_time;
    2278     realtype* scale_ptr = 0;
    2279     jmi_block_solver_options_t* bsop = block->options;
    2280     int use_scaling_flag = bsop->residual_equation_scaling_mode;
    2281 
    2282 
    2283     solver->kin_scale_update_time = curtime;
    2284     block->force_rescaling = 0;
    2285    
    2286     if(bsop->residual_equation_scaling_mode != jmi_residual_scaling_none) {
    2287         /* Zero out the scales initially if we're modify this. */
    2288         N_VConst_Serial(0,solver->kin_f_scale);
    2289     }
    2290 
    2291     scale_ptr = N_VGetArrayPointer(solver->kin_f_scale);
    2292 
    2293     /* Form scaled Jacobian as needed for automatic scaling and condition number checking*/
    2294     if (bsop->residual_equation_scaling_mode == jmi_residual_scaling_auto ||
    2295         bsop->residual_equation_scaling_mode == jmi_residual_scaling_aggressive_auto ||
    2296         bsop->residual_equation_scaling_mode == jmi_residual_scaling_full_jacobian_auto ||
    2297         bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid ||
    2298         bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual)
    2299     {
    2300         kin_char_log(solver, 's');
    2301         for (i = 0; i < N; i++) {
    2302             int j;
    2303             /* column scaling is formed by max(abs(nominal), abs(actual_value)) */
    2304             realtype xscale = RAbs(block->nominal[i]);
    2305             realtype x = RAbs(block->x[i]);
    2306             realtype* col_ptr;
    2307             realtype* scaled_col_ptr;
    2308             if(x < xscale) x = xscale;
    2309             if(x < tol) x = tol;
    2310             col_ptr = DENSE_COL(solver->J, i);
    2311             scaled_col_ptr = DENSE_COL(solver->J_scale, i);
    2312 
    2313             /* row scaling is product of Jac entry and column scaling */
    2314             for (j = 0; j < N; j++) {
    2315                 realtype dres = col_ptr[j];
    2316                 realtype fscale;
    2317                 fscale = dres * x;
    2318                 scaled_col_ptr[j] = fscale;
    2319                 scale_ptr[j] = MAX(scale_ptr[j], RAbs(fscale));
    2320             }
    2321         }
    2322     }
    2323 
    2324     /* put in manual non-zero scales in hybrid & manual cases */
    2325     if(    (bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual)
    2326         || (bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid))    {
    2327             for(i = 0; i < N; i++) {
    2328                 if(solver->residual_nominal[i] != 0) {
    2329                     scale_ptr[i] = solver->residual_nominal[i];
    2330                 }
    2331             }
    2332     }
    2333 
    2334     if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_auto
    2335         || bsop->residual_equation_scaling_mode == jmi_residual_scaling_aggressive_auto
    2336         || bsop->residual_equation_scaling_mode == jmi_residual_scaling_full_jacobian_auto
    2337         || bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid) {
    2338         solver->using_max_min_scaling_flag = 0; /* NOT using max/min scaling */
    2339         /* check that scaling factors have reasonable magnitude */
    2340         for(i = 0; i < N; i++) {
    2341             if(scale_ptr[i] < 1/bsop->max_residual_scaling_factor) {
    2342                 int j, maxInJacIndex = 0;
    2343                 realtype **jac = solver->J_scale->cols;
    2344                 realtype maxInJac = RAbs(jac[0][i]);
    2345                 solver->using_max_min_scaling_flag = 1; /* Using maximum scaling, singular Jacobian? */
    2346                 for(j = 0; j < N; j++) {
    2347                     if(RAbs(maxInJac) < RAbs(jac[j][i])) {
    2348                         maxInJac = jac[j][i];
    2349                         maxInJacIndex = j;
    2350                     }
    2351                 }
    2352                
    2353                 scale_ptr[i] = bsop->max_residual_scaling_factor;
    2354                 jmi_log_node(block->log, logWarning, "MaxScalingUsed", "Poor scaling in <block: %s>. "
    2355                     "Consider changes in the model. Partial derivative of <equation: %I> with respect to <Iter: #r%d#> is <dRes_dIter: %g> (<dRes_dIter_scaled: %g>).",
    2356                     block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(solver->J, i, maxInJacIndex) ,maxInJac);
    2357             }
    2358             else if(scale_ptr[i] > 1/bsop->min_residual_scaling_factor) {
    2359                 int j, maxInJacIndex = 0;
    2360                 realtype **jac = solver->J_scale->cols;
    2361                 realtype maxInJac = RAbs(jac[0][i]);
    2362                 /* Likely not a problem: solver->using_max_min_scaling_flag = 1; -- Using minimum scaling */
    2363                 for(j = 0; j < N; j++) {
    2364                     if(RAbs(maxInJac) < RAbs(jac[j][i])) {
    2365                         maxInJac = jac[j][i];
    2366                         maxInJacIndex = j;
    2367                     }
    2368                 }
    2369                
    2370                 if (scale_ptr[i] > solver->residual_heuristic_nominal[i]/bsop->min_residual_scaling_factor) {
    2371                     scale_ptr[i] = bsop->min_residual_scaling_factor/solver->residual_heuristic_nominal[i];
    2372                     jmi_log_node(block->log, logWarning, "MinScalingUsed", "Poor scaling in <block: %s>. "
    2373                         "Consider changes in the model. Partial derivative of <equation: %I> with respect to <Iter: #r%d#> is <dRes_dIter: %g> (<dRes_dIter_scaled: %g>).",
    2374                         block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(solver->J, i, maxInJacIndex) ,maxInJac);
    2375                 } else {
    2376                     scale_ptr[i] = 1/scale_ptr[i];
    2377                     jmi_log_node(block->log, logWarning, "MinScalingExceeded", "Poor scaling in <block: %s> but will allow it based on structure of block. "
    2378                         "Consider changes in the model. Partial derivative of <equation: %I> with respect to <Iter: #r%d#> is <dRes_dIter: %g> (<dRes_dIter_scaled: %g>).",
    2379                         block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(solver->J, i, maxInJacIndex) ,maxInJac);
    2380                 }
    2381             }
    2382             else
    2383                 scale_ptr[i] = 1/scale_ptr[i];
    2384         }
    2385 
    2386         if (block->callbacks->log_options.log_level >= 4) {
    2387             jmi_log_node_t outer = jmi_log_enter_fmt(block->log, logInfo, "ResidualScalingUpdated", "<block:%s>", block->label);
    2388             if (block->callbacks->log_options.log_level >= 5) {
    2389                 jmi_log_node_t inner = jmi_log_enter_vector_(block->log, outer, logInfo, "scaling");
    2390                 realtype* res = scale_ptr;
    2391                 for (i=0;i<N;i++) jmi_log_real_(block->log, 1/res[i]);
    2392                 jmi_log_leave(block->log, inner);
    2393             }
    2394             jmi_log_leave(block->log, outer);
    2395         }
    2396     } else if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) {
    2397         for(i = 0; i < N; i++) {
    2398             scale_ptr[i] = 1/scale_ptr[i];
    2399         }
    2400         /* Scaling is not updated */
    2401         if (block->callbacks->log_options.log_level >= 4) {
    2402             jmi_log_node_t outer = jmi_log_enter_fmt(block->log, logInfo, "ResidualScalingUpdated", "<block:%s>", block->label);
    2403             if (block->callbacks->log_options.log_level >= 5) {
    2404                 jmi_log_node_t inner = jmi_log_enter_vector_(block->log, outer, logInfo, "scaling");
    2405                 realtype* res = scale_ptr;
    2406                 for (i=0;i<N;i++) jmi_log_real_(block->log, 1/res[i]);
    2407                 jmi_log_leave(block->log, inner);
    2408             }
    2409             jmi_log_leave(block->log, outer);
    2410         }
    2411     }
    2412 
    2413     if (solver->using_max_min_scaling_flag && !solver->J_is_singular_flag) {
    2414         realtype cond = jmi_calculate_jacobian_condition_number(block);
    2415 
    2416         jmi_log_node(block->log, logInfo, "Regularization",
    2417             "Calculated condition number in <block: %s>. Regularizing if <cond: %E> is greater than <regtol: %E>", block->label, cond, solver->kin_reg_tol);
    2418         if (cond > solver->kin_reg_tol) {
    2419             if(N > 1 && solver->handling_of_singular_jacobian_flag == JMI_REGULARIZATION) {
    2420                 int info;
    2421                 jmi_kinsol_reg_matrix(block);
    2422                 dgetrf_(  &block->n, &block->n, solver->JTJ->data, &block->n, solver->lapack_ipiv, &info);
    2423             }
    2424             solver->J_is_singular_flag = 1;
    2425         }
    2426     }
    2427 
    2428 
    2429     /* estimate condition number of the scaled jacobian
    2430     and scale function tolerance with it. */
    2431     if((N > 1) && bsop->check_jac_cond_flag){
    2432         realtype* scaled_col_ptr;
    2433         char norm = 'I';
    2434         double Jnorm = 1.0, Jcond = 1.0;
    2435         int info;
    2436         if(use_scaling_flag) {
    2437             for(i = 0; i < N; i++){
    2438                 int j;
    2439                 scaled_col_ptr = DENSE_COL(solver->J_scale, i);
    2440                 for(j = 0; j < N; j++){
    2441                     scaled_col_ptr[j] = scaled_col_ptr[j] * scale_ptr[j];
    2442                 }
    2443             }
    2444         } else {
    2445             DenseCopy(solver->J, solver->J_scale);
    2446         }
    2447 
    2448         dgetrf_(  &N, &N, solver->J_scale->data, &N, solver->lapack_iwork, &info);
    2449         if(info > 0) {
    2450             jmi_log_node(block->log, logWarning, "SingularJacobian",
    2451                 "Singular Jacobian detected when checking condition number in <block:%s>. Solver may fail to converge.", block->label);
    2452         }
    2453         else {
    2454             dgecon_(&norm, &N, solver->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
    2455 
    2456             if(tol * Jcond < UNIT_ROUNDOFF) {
    2457                 jmi_log_node_t node = jmi_log_enter_fmt(block->log, logWarning, "IllConditionedJacobian",
    2458                     "<JacobianInverseConditionEstimate:%E> Solver may fail to converge in <block: %s>.", Jcond, block->label);
    2459                 if (block->callbacks->log_options.log_level >= 4) {
    2460                     jmi_log_reals(block->log, node, logWarning, "ivs", N_VGetArrayPointer(solver->kin_y), block->n);
    2461                 }
    2462                 jmi_log_leave(block->log, node);
    2463 
    2464             }
    2465             else {
    2466                 jmi_log_node(block->log, logInfo, "JacobianCondition",
    2467                     "<JacobianInverseConditionEstimate:%E>", Jcond);
    2468             }
    2469         }
    2470     }
    2471     return;
    24722296}
    24732297
     
    24932317    solver->kin_y = N_VMake_Serial(n, block->x);
    24942318    solver->kin_y_scale = N_VNew_Serial(n);
    2495     solver->kin_f_scale = N_VNew_Serial(n);
    24962319    solver->gradient  = N_VNew_Serial(n);   
    24972320    solver->last_residual = N_VNew_Serial(n);
    2498     solver->residual_nominal = (realtype*)calloc(n+1,sizeof(realtype));
    2499     solver->residual_heuristic_nominal = (realtype*)calloc(n+1,sizeof(realtype));
    2500     solver->kin_scale_update_time = -1.0;
    25012321    solver->kin_jac_update_time = -1.0;
    25022322    /*NOTE: it'd be nice to use "jmi->newton_tolerance" here
     
    25052325    solver->kin_ftol = UROUND; /* block->options->min_tol; */
    25062326    solver->kin_stol = block->options->min_tol;
    2507     solver->using_max_min_scaling_flag = 0; /* Not using max/min scaling */
    25082327    solver->has_compression_setup_flag = FALSE;
    25092328    solver->is_first_newton_solve_flag = TRUE;
    25102329    solver->current_nni = 0;
    25112330   
    2512     solver->J = NewDenseMat(n ,n);
    25132331    solver->JTJ = NewDenseMat(n ,n);
    25142332    solver->J_LU = NewDenseMat(n ,n);
    2515     solver->J_scale = NewDenseMat(n ,n);
    25162333    solver->J_sing = NewDenseMat(n, n);
    25172334    solver->J_SVD_U = NewDenseMat(n, n);
     
    25492366    /* Initialize scaling to 1.0 - defaults */
    25502367    N_VConst_Serial(1.0,solver->kin_y_scale);
    2551 
    2552     /* Initial equation scaling is 1.0 */
    2553     N_VConst_Serial(1.0,solver->kin_f_scale);
    25542368               
    25552369    flag = KINInit(solver->kin_mem, kin_f, solver->kin_y); /*Initialize Kinsol*/
     
    26372451    N_VDestroy_Serial(solver->kin_y);
    26382452    N_VDestroy_Serial(solver->kin_y_scale);
    2639     N_VDestroy_Serial(solver->kin_f_scale);
    26402453    N_VDestroy_Serial(solver->gradient);
    26412454    N_VDestroy_Serial(solver->last_residual);
    2642     free(solver->residual_nominal);
    2643     free(solver->residual_heuristic_nominal);
    2644     DestroyMat(solver->J);
    26452455    DestroyMat(solver->JTJ);
    26462456    DestroyMat(solver->J_LU);
    2647     DestroyMat(solver->J_scale);
    26482457    DestroyMat(solver->J_sing);
    26492458    DestroyMat(solver->J_SVD_U);
     
    27712580    }
    27722581    jmi_kinsol_solver_print_solve_start(block, &topnode);
    2773     flag = KINSol(solver->kin_mem, solver->kin_y, strategy, solver->kin_y_scale, solver->kin_f_scale);
     2582    flag = KINSol(solver->kin_mem, solver->kin_y, strategy, solver->kin_y_scale, block->f_scale);
    27742583    if(flag == KIN_INITIAL_GUESS_OK) { /* last_fnorm not up to date in this case */
    27752584        double max_residual = 0.0;
     
    27772586        realtype* f;
    27782587        char msg[256];
    2779         realtype* residual_scaling_factors = NV_DATA_S(solver->kin_f_scale);
     2588        realtype* residual_scaling_factors = NV_DATA_S(block->f_scale);
    27802589
    27812590        /*kin_f(solver->kin_y, solver->work_vector, block); */
     
    27942603        solver->last_max_residual= max_residual;
    27952604        solver->last_max_residual_index = max_index;
    2796         solver->last_fnorm = N_VWL2Norm(kin_mem->kin_fval, solver->kin_f_scale);
     2605        solver->last_fnorm = N_VWL2Norm(kin_mem->kin_fval, block->f_scale);
    27972606        kin_mem->kin_fnorm = solver->last_fnorm;
    27982607        sprintf(msg, "nni = %4ld   nfe = %6ld   fnorm = %26.16g", kin_mem->kin_nni, kin_mem->kin_nfe, solver->last_fnorm);
     
    28102619    } else if (flag == KIN_LINESEARCH_NONCONV || flag == KIN_STEP_LT_STPTOL) {
    28112620        realtype fnorm;
    2812         N_VProd(solver->kin_f_scale, kin_mem->kin_fval, solver->work_vector);
     2621        N_VProd(block->f_scale, kin_mem->kin_fval, solver->work_vector);
    28132622        fnorm =N_VMaxNorm(solver->work_vector);
    28142623        if((fnorm <= solver->kin_stol && !block->options->solver_exit_criterion_mode == jmi_exit_criterion_step_residual)
     
    28262635        if(block->options->check_jac_cond_flag) {
    28272636            int i, info, N = block->n;
    2828             realtype* scale_ptr = N_VGetArrayPointer(solver->kin_f_scale);
     2637            realtype* scale_ptr = N_VGetArrayPointer(block->f_scale);
    28292638            realtype tol = solver->kin_stol;
    28302639            if(block->options->residual_equation_scaling_mode !=0) {
    28312640                for(i = 0; i < N; i++){
    28322641                    int j;
    2833                     realtype* scaled_col_ptr = DENSE_COL(solver->J_scale, i);
    2834                     realtype* col_ptr = DENSE_COL(solver->J, i);
     2642                    realtype* scaled_col_ptr = DENSE_COL(block->J_scale, i);
     2643                    realtype* col_ptr = DENSE_COL(block->J, i);
    28352644                    realtype xscale = RAbs(block->nominal[i]);
    28362645                    realtype x = RAbs(block->x[i]);
     
    28422651                }
    28432652            } else {
    2844                 DenseCopy(solver->J, solver->J_scale);
    2845             }
    2846 
    2847             dgetrf_(  &N, &N, solver->J_scale->data, &N, solver->lapack_iwork, &info);
     2653                DenseCopy(block->J, block->J_scale);
     2654            }
     2655
     2656            dgetrf_(  &N, &N, block->J_scale->data, &N, solver->lapack_iwork, &info);
    28482657            if(info > 0) {
    28492658                jmi_log_node(block->log, logWarning, "SingularJacobian",
     
    28532662                char norm = 'I';
    28542663                double Jnorm = 1.0, Jcond = 1.0;
    2855                 dgecon_(&norm, &N, solver->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
     2664                dgecon_(&norm, &N, block->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
    28562665
    28572666                if(tol * Jcond < UNIT_ROUNDOFF) {
     
    28932702   
    28942703    if(block->init) {
    2895         jmi_setup_f_residual_scaling(block);
    28962704        flag = jmi_kinsol_init(block);
    28972705        if(flag) return flag;
     
    29252733
    29262734    /* update the scaling only once per time step */
    2927     if(block->init || (block->options->rescale_each_step_flag && (curtime > solver->kin_scale_update_time)) || block->force_rescaling) {
     2735    if(block->init || (block->options->rescale_each_step_flag && (curtime > block->scale_update_time)) || block->force_rescaling) {
    29282736        jmi_update_f_scale(block);
     2737        jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block);
    29292738    }
    29302739   
     
    29712780        solver->use_steepest_descent_flag = 1;
    29722781       
    2973         flag = KINSol(solver->kin_mem, solver->kin_y, KIN_LINESEARCH, solver->kin_y_scale, solver->kin_f_scale);
     2782        flag = KINSol(solver->kin_mem, solver->kin_y, KIN_LINESEARCH, solver->kin_y_scale, block->f_scale);
    29742783        if(flag == KIN_INITIAL_GUESS_OK) {
    29752784            flag = KIN_SUCCESS;
     
    29982807        /* Update the scaling  */
    29992808        jmi_update_f_scale(block);
     2809        jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block);
    30002810       
    30012811        /* For the second solve set tight tolerance on the step and specified tolerance on the residual */
     
    31212931            DenseCopy(solver->saved_state->J_modified, solver->J_LU);
    31222932    }
    3123     DenseCopy(solver->saved_state->J, solver->J);
    3124    
    3125     solver->kin_scale_update_time = solver->saved_state->kin_scale_update_time;
     2933    DenseCopy(solver->saved_state->J, block->J);
     2934   
     2935    block->scale_update_time = solver->saved_state->kin_scale_update_time;
    31262936    block->force_rescaling        = solver->saved_state->force_rescaling;
    31272937
     
    31302940    solver->J_is_singular_flag = solver->saved_state->J_is_singular_flag;
    31312941    ((struct KINMemRec *)solver->kin_mem)->kin_nnilset = ((struct KINMemRec *)solver->kin_mem)->kin_nni - solver->saved_state->mbset;
    3132     memcpy(N_VGetArrayPointer(solver->kin_f_scale), N_VGetArrayPointer(solver->saved_state->kin_f_scale), block->n*sizeof(realtype));
     2942    memcpy(N_VGetArrayPointer(block->f_scale), N_VGetArrayPointer(solver->saved_state->kin_f_scale), block->n*sizeof(realtype));
    31332943    memcpy(N_VGetArrayPointer(solver->kin_y_scale), N_VGetArrayPointer(solver->saved_state->kin_y_scale), block->n*sizeof(realtype));
    31342944    memcpy(solver->lapack_ipiv, solver->saved_state->lapack_ipiv, (block->n+2)*sizeof(int));
     
    31392949        jmi_log_ints(block->log, node, logInfo, "is_singular", &(block->force_rescaling), 1);
    31402950        jmi_log_ints(block->log, node, logInfo, "force_rescaling", &(solver->J_is_singular_flag), 1);
    3141         jmi_log_reals(block->log, node, logInfo, "residual_scaling_update_time", &(solver->kin_scale_update_time), 1);
     2951        jmi_log_reals(block->log, node, logInfo, "residual_scaling_update_time", &(block->scale_update_time), 1);
    31422952        jmi_log_reals(block->log, node, logInfo, "iv_scaling", N_VGetArrayPointer(solver->kin_y_scale), block->n);
    3143         jmi_log_reals(block->log, node, logInfo, "residual_scaling", N_VGetArrayPointer(solver->kin_f_scale), block->n);
     2953        jmi_log_reals(block->log, node, logInfo, "residual_scaling", N_VGetArrayPointer(block->f_scale), block->n);
    31442954        jmi_log_leave(block->log, node);
    31452955    }
     
    31802990            DenseCopy(solver->J_LU, solver->saved_state->J_modified);
    31812991        }
    3182         DenseCopy(solver->J, solver->saved_state->J);
     2992        DenseCopy(block->J, solver->saved_state->J);
    31832993       
    3184         solver->saved_state->kin_scale_update_time = solver->kin_scale_update_time;
     2994        solver->saved_state->kin_scale_update_time = block->scale_update_time;
    31852995        solver->saved_state->force_rescaling       = block->force_rescaling;
    31862996        solver->saved_state->force_new_J_flag      = solver->force_new_J_flag;
     
    31892999        solver->saved_state->J_is_singular_flag = solver->J_is_singular_flag;
    31903000        solver->saved_state->mbset = ((struct KINMemRec *)solver->kin_mem)->kin_nni - ((struct KINMemRec *)solver->kin_mem)->kin_nnilset;
    3191         memcpy(N_VGetArrayPointer(solver->saved_state->kin_f_scale),N_VGetArrayPointer(solver->kin_f_scale), block->n*sizeof(realtype));
     3001        memcpy(N_VGetArrayPointer(solver->saved_state->kin_f_scale),N_VGetArrayPointer(block->f_scale), block->n*sizeof(realtype));
    31923002        memcpy(N_VGetArrayPointer(solver->saved_state->kin_y_scale),N_VGetArrayPointer(solver->kin_y_scale), block->n*sizeof(realtype));
    31933003        memcpy(solver->saved_state->lapack_ipiv, solver->lapack_ipiv, (block->n+2)*sizeof(int));
     
    31983008            jmi_log_ints(block->log, node, logInfo, "is_singular", &(solver->J_is_singular_flag), 1);
    31993009            jmi_log_ints(block->log, node, logInfo, "force_rescaling", &(solver->J_is_singular_flag), 1);
    3200             jmi_log_reals(block->log, node, logInfo, "residual_scaling_update_time", &(solver->kin_scale_update_time), 1);
     3010            jmi_log_reals(block->log, node, logInfo, "residual_scaling_update_time", &(block->scale_update_time), 1);
    32013011            jmi_log_reals(block->log, node, logInfo, "iv_scaling", N_VGetArrayPointer(solver->kin_y_scale), block->n);
    3202             jmi_log_reals(block->log, node, logInfo, "residual_scaling", N_VGetArrayPointer(solver->kin_f_scale), block->n);
     3012            jmi_log_reals(block->log, node, logInfo, "residual_scaling", N_VGetArrayPointer(block->f_scale), block->n);
    32033013            jmi_log_leave(block->log, node);
    32043014        }
     
    32063016    return 0;
    32073017}
    3208 
  • trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.h

    r8384 r8512  
    5858int jmi_kinsol_restore_state(jmi_block_solver_t* block);
    5959
    60 /**< \brief Retrieve residual scales used in Kinsol solver */
    61 double* jmi_kinsol_solver_get_f_scales(jmi_kinsol_solver_t* solver);
    62 
    6360/**< \brief Convert Kinsol return flag to readable name */
    6461const char *jmi_kinsol_flag_to_name(int flag);
     
    8784
    8885    N_Vector kin_y_scale;          /**< \brief Work vector for Kinsol scaling of y */
    89     N_Vector kin_f_scale;          /**< \brief Work vector for Kinsol scaling of f */
    9086    N_Vector gradient;              /**< \brief Steepest descent direction */
    91     realtype* residual_nominal;      /**< \brief Vector for reading in manual scaling factors for f  */
    92     realtype* residual_heuristic_nominal; /**< \brief Vector for reading in heuristic scaling factors for f  */
    93     realtype kin_scale_update_time; /**< \brief The last time when Kinsol scale was updated */
    9487    realtype kin_jac_update_time; /**< \brief The last time when Jacobian was updated */
    9588    realtype kin_ftol;             /**< \brief Tolerance for F */
     
    9790    realtype kin_reg_tol;          /**< \brief Regularization tolerance */
    9891   
    99     DlsMat J;                       /**< \brief The Jacobian matrix  */   
    10092    DlsMat JTJ;                     /**< \brief The Transpose(J).J used if J is singular */
    10193    int J_is_singular_flag;         /**< \brief A flag indicating that J is singular. Regularized JTJ is setup */
    10294    int use_steepest_descent_flag;  /**< \brief A flag indicating that steepest descent and not Newton direction should be used */
    10395    int force_new_J_flag;           /**< \brief A flag indicating that J needs to be recalculated */
    104     int using_max_min_scaling_flag; /**< \brief A flag indicating if either the maximum scaling is used of the minimum */
    10596    int updated_jacobian_flag;      /**< \brief A flag indicating if an updated Jacobian is used to solve the system */
    10697    int handling_of_singular_jacobian_flag; /**< \brief A flag for determining how singular systems should be treated */
    10798    DlsMat J_LU;                    /**< \brief Jacobian matrix/it's LU decomposition */
    108     DlsMat J_scale;                 /**< \brief Jacobian matrix scaled with xnorm for used for fnorm calculation */
    10999    DlsMat J_sing;                  /**< \brief Jacobian matrix/it's right singular vectors */
    110100    DlsMat J_SVD_U;                 /**< \brief The left singular vectors */
  • trunk/RuntimeLibrary/src/jmi/jmi_linear_solver.c

    r8052 r8512  
    4040    /* Initialize work vectors.*/
    4141    solver->factorization = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t));
    42     solver->jacobian = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t));
    4342    solver->dependent_set = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t));
    4443    solver->jacobian_temp = (jmi_real_t*)calloc(2*n_x*n_x,sizeof(jmi_real_t));
     
    8281    for (i = 0; i < n_x; i++) { solver->dependent_set[(n_x-1)*n_x+i] = 0; }
    8382   
    84     memcpy(solver->jacobian_temp, solver->jacobian, n_x*n_x*sizeof(jmi_real_t));
     83    memcpy(solver->jacobian_temp, block->J->data, n_x*n_x*sizeof(jmi_real_t));
    8584    /* Compute the null space of A */
    8685    /*
     
    197196    if (solver->n_extra_rows > 0) {
    198197        for (j = 0; j < n_x; j++) {
    199             memcpy(&solver->jacobian_temp[j*(n_x+solver->n_extra_rows)], &solver->jacobian[j*n_x], n_x*sizeof(jmi_real_t));
     198            memcpy(&solver->jacobian_temp[j*(n_x+solver->n_extra_rows)], &block->J->data[j*n_x], n_x*sizeof(jmi_real_t));
    200199            memcpy(&solver->jacobian_temp[j*(n_x+solver->n_extra_rows)+n_x], &solver->jacobian_extension[j*n_x], solver->n_extra_rows*sizeof(jmi_real_t));
    201200        }
     
    246245             A regularization strategy for simple cases singular jac should be introduced.
    247246          */
    248         info = block->F(block->problem_data,NULL,solver->jacobian,JMI_BLOCK_EVALUATE_JACOBIAN);
    249         memcpy(solver->factorization, solver->jacobian, n_x*n_x*sizeof(jmi_real_t));
     247        info = block->F(block->problem_data,NULL,block->J->data,JMI_BLOCK_EVALUATE_JACOBIAN);
     248        memcpy(solver->factorization, block->J->data, n_x*n_x*sizeof(jmi_real_t));
    250249        if(info) {
    251250            if(block->init) {
     
    275274            for (j = 0; j < n_x; j++) {
    276275                /* Unrecoverable error*/
    277                 if ( solver->jacobian[i*n_x+j] - solver->jacobian[i*n_x+j] != 0) {
     276                if ( block->J->data[i*n_x+j] - block->J->data[i*n_x+j] != 0) {
    278277                    jmi_log_node(block->log, logError, "NaNOutput", "Not a number in the Jacobian <row: %I> <col: %I> from <block: %s>",
    279278                            i,j, block->label);
     
    289288                                     "Linear solver invoked for <block:%s>", block->label);
    290289        jmi_log_reals(block->log, destnode, logInfo, "ivs", block->x, block->n);
    291         jmi_log_real_matrix(block->log, destnode, logInfo, "A", solver->jacobian, block->n, block->n);
     290        jmi_log_real_matrix(block->log, destnode, logInfo, "A", block->J->data, block->n, block->n);
    292291    }
    293292
     
    386385            dgelss_(&n_rows, &n_x, &i, solver->jacobian_temp, &n_rows, solver->rhs, &n_rows ,solver->singular_values, &rcond, &rank, solver->rwork, &iwork, &info);
    387386        } else {
    388             memcpy(solver->jacobian_temp, solver->jacobian, n_x*n_x*sizeof(jmi_real_t));
     387            memcpy(solver->jacobian_temp, block->J->data, n_x*n_x*sizeof(jmi_real_t));
    389388            dgelss_(&n_x, &n_x, &i, solver->jacobian_temp, &n_x, solver->rhs, &n_x ,solver->singular_values, &rcond, &rank, solver->rwork, &iwork, &info);
    390389        }
     
    395394        }
    396395       
    397         if(block->callbacks->log_options.log_level >= 5){
    398             jmi_log_node(block->log, logInfo, "Info", "Successfully calculated the minimum norm solution to the linear system in <block: %s>", block->label);
    399         }
    400396    }else{
    401397        /*
     
    449445    block->F(block->problem_data,block->x, solver->rhs, JMI_BLOCK_EVALUATE);
    450446
     447    /* Check if the calculated solution from minimum norm is a valid solution to the original problem */
     448    if(solver->singular_jacobian==1) {
     449        double scaledMaxNorm;
     450        jmi_update_f_scale(block);
     451        scaledMaxNorm = calculate_scaled_residual_norm(solver->rhs, N_VGetArrayPointer(block->f_scale), block->n);
     452        if(scaledMaxNorm <= block->options->res_tol) {
     453            if(block->callbacks->log_options.log_level >= 5){
     454                jmi_log_node(block->log, logInfo, "Info", "Successfully calculated the minimum norm solution to the linear system in <block: %s>", block->label);
     455            }
     456        } else {
     457            info = -1;
     458            destnode = jmi_log_enter_fmt(block->log, logError, "UnsolveableLinearSystem", "Failed to calculate a valid minimum norm solution to the linear system in <block: %s> at <t: %f>", block->label, block->cur_time);
     459            jmi_log_reals(block->log, destnode, logError, "residual", solver->rhs, block->n);
     460            jmi_log_reals(block->log, destnode, logError, "scaled_max_norm", &(scaledMaxNorm), 1);
     461            jmi_log_leave(block->log, destnode);
     462        }
     463
     464    }
     465
    451466    return info==0? 0: -1;
    452467}
     
    456471    free(solver->ipiv);
    457472    free(solver->factorization);
    458     free(solver->jacobian);
    459473    free(solver->singular_values);
    460474    free(solver->singular_vectors);
  • trunk/RuntimeLibrary/src/jmi/jmi_linear_solver.h

    r7216 r8512  
    5858    int* ipiv;                     /**< \brief Work vector needed for dgesv */
    5959    jmi_real_t* factorization;      /**< \brief Matrix for storing the Jacobian factorization */
    60     jmi_real_t* jacobian;         /**< \brief Matrix for storing the Jacobian */
    6160    jmi_real_t* jacobian_temp;         /**< \brief Matrix for storing the Jacobian */
    6261    jmi_real_t* singular_values;  /**< \brief Vector for the singular values of the Jacobian */
Note: See TracChangeset for help on using the changeset viewer.