Changeset 8512

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:

• Python/src/tests_jmodelica/files/Modelica/Singular.mo (modified) (1 diff)
• Python/src/tests_jmodelica/simulation/test_assimulo_interface_fmi.py (modified) (2 diffs)
• RuntimeLibrary/src/jmi/jmi_block_solver.c (modified) (9 diffs)
• RuntimeLibrary/src/jmi/jmi_block_solver.h (modified) (2 diffs)
• RuntimeLibrary/src/jmi/jmi_block_solver_impl.h (modified) (2 diffs)
• RuntimeLibrary/src/jmi/jmi_kinsol_solver.c (modified) (56 diffs)
• RuntimeLibrary/src/jmi/jmi_kinsol_solver.h (modified) (3 diffs)
• RuntimeLibrary/src/jmi/jmi_linear_solver.c (modified) (10 diffs)
• RuntimeLibrary/src/jmi/jmi_linear_solver.h (modified) (1 diff)

trunk/Python/src/tests_jmodelica/files/Modelica/Singular.mo

@@ -205 +205 @@
         (sin(x*y*z)+z) = 0;
     end NonLinear5;
+    
+    model NoMinimumNormSolution
+        Real x,y,z,v;
+        parameter Real a11 = 1;
+        parameter Real a12 = 2;
+        parameter Real a13 = 3;
+        parameter Real a21 = 1;
+        parameter Real a22 = 1;
+        parameter Real a23 = 1;
+        Real a31(start=0);
+        Real a32(start=0);
+        Real a33(start=0);
+        parameter Real b[3] = {1,2,3};
+    equation
+        a33=a31+a32;
+        a33=a31^2+a32^2;
+        a33=a31-a32;
+        a11*x+a12*y+a13*z = b[1];
+        a21*x+a22*y+a23*z = b[2];
+        a31*x+a32*y+a33*z = b[3];
+        der(v) = time;
+    end NoMinimumNormSolution;
 
 end Singular;

trunk/Python/src/tests_jmodelica/simulation/test_assimulo_interface_fmi.py

@@ -621 +621 @@
         compile_fmu("Singular.NonLinear4", file_name)
         compile_fmu("Singular.NonLinear5", file_name)
+        compile_fmu("Singular.NoMinimumNormSolution", file_name)
     
     @testattr(stddist = True)
@@ -749 +750 @@
         nose.tools.assert_almost_equal(res["z"][0] ,0.000000000)
         nose.tools.assert_almost_equal(res["z"][-1] ,-1.000000000)
+    
+    @testattr(stddist = True)
+    def test_no_valid_minimum_norm_sol(self):
+        model = load_fmu("NoMinimumNormSolution.fmu", log_level=3)
+        model.set("_log_level", 3)
+        model.set_log_level(log_level)
+        nose.tools.assert_raises(FMUException, model.initialize)
 
 class Test_FMI_ODE_CS_2:

trunk/RuntimeLibrary/src/jmi/jmi_block_solver.c

@@ -28 +28 @@
 #include <string.h>
 #include <math.h>
+
+#include <sundials/sundials_math.h>
+#include <sundials/sundials_direct.h>
+#include <nvector/nvector_serial.h>
+#include <kinsol/kinsol_direct.h>
+#include <kinsol/kinsol_impl.h>
+#include <sundials/sundials_dense.h>
 
 #include "jmi_log.h"
@@ -56 +63 @@
                            jmi_block_solver_options_t* options,
                            void* problem_data){
+    int i;
     jmi_block_solver_t* block_solver = (jmi_block_solver_t*)calloc(1, sizeof(jmi_block_solver_t));
     if(!block_solver) {
@@ -74 +82 @@
 
     block_solver->dx=(jmi_real_t*)calloc(n,sizeof(jmi_real_t));;                /**< \brief Work vector for the seed vector */
+
+    block_solver->f_scale = N_VNew_Serial(n);
+    block_solver->scale_update_time = -1.0;
+    if(n>0) {
+        block_solver->J = NewDenseMat(n ,n);
+        block_solver->J_scale = NewDenseMat(n ,n);
+    }
 
     block_solver->res = (jmi_real_t*)calloc(n,sizeof(jmi_real_t));
@@ -84 +99 @@
     block_solver->max = (jmi_real_t*)calloc(n,sizeof(jmi_real_t));
     block_solver->nominal = (jmi_real_t*)calloc(n,sizeof(jmi_real_t));
+    block_solver->residual_nominal = (jmi_real_t*)calloc(n+1,sizeof(jmi_real_t));
+    block_solver->residual_heuristic_nominal = (jmi_real_t*)calloc(n+1,sizeof(jmi_real_t));
     block_solver->initial = (jmi_real_t*)calloc(n,sizeof(jmi_real_t));
     block_solver->start_set = (jmi_real_t*)calloc(n,sizeof(jmi_real_t));
@@ -97 +114 @@
 
     block_solver->cur_time = 0;
+    block_solver->using_max_min_scaling_flag = 0; /* Not using max/min scaling */
+
+        /* Initial equation scaling is 1.0 */
+    N_VConst_Serial(1.0,block_solver->f_scale);
+
+    /* Make sure that residual_nominals and heuristic nominals always have values */
+    for(i=0; i<n; i++) {
+        block_solver->residual_nominal[i] = 1.0;
+        block_solver->residual_heuristic_nominal[i] = 1.0;
+    }
+
     switch(options->solver) {
         case JMI_KINSOL_SOLVER: {
@@ -196 +224 @@
     free(block_solver->dx);
 
+    N_VDestroy_Serial(block_solver->f_scale);
+    if(block_solver->n > 0) {
+        DestroyMat(block_solver->J);
+        DestroyMat(block_solver->J_scale);
+    }
+
     free(block_solver->res);
     free(block_solver->dres);
@@ -205 +239 @@
     free(block_solver->max);
     free(block_solver->nominal);
+    free(block_solver->residual_nominal);
+    free(block_solver->residual_heuristic_nominal);
     free(block_solver->initial);
     free(block_solver->value_references);
@@ -297 +333 @@
         block_solver->F(block_solver->problem_data,real_vrs,block_solver->res,JMI_BLOCK_VALUE_REFERENCE);
         block_solver->F(block_solver->problem_data,block_solver->start_set,block_solver->res,JMI_BLOCK_START_SET);
+        jmi_setup_f_residual_scaling(block_solver);
 
         for (i=0;i<block_solver->n;i++) {
@@ -710 +747 @@
 }
 
+/* 
+    Compute appropriate equation scaling and function tolerance based on Jacobian J,
+    nominal values (block->nominal) and current point (block->x).
+    Store result in block->f_scale.
+*/
+void jmi_update_f_scale(jmi_block_solver_t *block) { 
+    int i, N = block->n;
+    realtype curtime = block->cur_time;
+    realtype* scale_ptr = 0;
+    jmi_block_solver_options_t* bsop = block->options;
+    int use_scaling_flag = bsop->residual_equation_scaling_mode;
+    realtype tol = bsop->res_tol;
+
+    block->scale_update_time = curtime;
+    block->force_rescaling = 0;
+    
+    if(bsop->residual_equation_scaling_mode != jmi_residual_scaling_none) {
+        /* Zero out the scales initially if we're modify this. */
+        N_VConst_Serial(0,block->f_scale);
+    }
+
+    scale_ptr = N_VGetArrayPointer(block->f_scale);
+
+    /* Form scaled Jacobian as needed for automatic scaling and condition number checking*/
+    if (bsop->residual_equation_scaling_mode == jmi_residual_scaling_auto || 
+        bsop->residual_equation_scaling_mode == jmi_residual_scaling_aggressive_auto ||
+        bsop->residual_equation_scaling_mode == jmi_residual_scaling_full_jacobian_auto ||
+        bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid ||
+        bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual)
+    {
+        for (i = 0; i < N; i++) {
+            int j;
+            /* column scaling is formed by max(abs(nominal), abs(actual_value)) */
+            realtype xscale = RAbs(block->nominal[i]);
+            realtype x = RAbs(block->x[i]);
+            realtype* col_ptr;
+            realtype* scaled_col_ptr;
+            if(x < xscale) x = xscale;
+            if(x < tol) x = tol;
+            col_ptr = DENSE_COL(block->J, i);
+            scaled_col_ptr = DENSE_COL(block->J_scale, i);
+
+            /* row scaling is product of Jac entry and column scaling */
+            for (j = 0; j < N; j++) {
+                realtype dres = col_ptr[j];
+                realtype fscale;
+                fscale = dres * x;
+                scaled_col_ptr[j] = fscale;
+                scale_ptr[j] = MAX(scale_ptr[j], RAbs(fscale));
+            }
+        }
+    }
+
+    /* put in manual non-zero scales in hybrid & manual cases */
+    if(    (bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) 
+        || (bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid))    {
+            for(i = 0; i < N; i++) {
+                if(block->residual_nominal[i] != 0) {
+                    scale_ptr[i] = block->residual_nominal[i];
+                }
+            }
+    }
+
+    if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_auto 
+        || bsop->residual_equation_scaling_mode == jmi_residual_scaling_aggressive_auto
+        || bsop->residual_equation_scaling_mode == jmi_residual_scaling_full_jacobian_auto
+        || bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid) {
+        block->using_max_min_scaling_flag = 0; /* NOT using max/min scaling */
+        /* check that scaling factors have reasonable magnitude */
+        for(i = 0; i < N; i++) {
+            if(scale_ptr[i] == 0.0) {
+                jmi_log_node(block->log, logInfo, "JacobianZeroRow", 
+                    "The derivatives for <residual: %d> are all equal 0. Using heuristic residual nominals for scaling in <block: %s>", i, block->label);
+                scale_ptr[i] = block->residual_heuristic_nominal[i];
+            } else if(scale_ptr[i] < 1/bsop->max_residual_scaling_factor) {
+                int j, maxInJacIndex = 0;
+                realtype **jac = block->J_scale->cols;
+                realtype maxInJac = RAbs(jac[0][i]); 
+                block->using_max_min_scaling_flag = 1; /* Using maximum scaling, singular Jacobian? */
+                for(j = 0; j < N; j++) {
+                    if(RAbs(maxInJac) < RAbs(jac[j][i])) {
+                        maxInJac = jac[j][i];
+                        maxInJacIndex = j;
+                    }
+                }
+                
+                scale_ptr[i] = bsop->max_residual_scaling_factor;
+                jmi_log_node(block->log, logWarning, "MaxScalingUsed", "Poor scaling in <block: %s>. "
+                    "Consider changes in the model. Partial derivative of <equation: %I> with respect to <Iter: #r%d#> is <dRes_dIter: %g> (<dRes_dIter_scaled: %g>).",
+                    block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(block->J, i, maxInJacIndex) ,maxInJac);
+            }
+            else if(scale_ptr[i] > 1/bsop->min_residual_scaling_factor) {
+                int j, maxInJacIndex = 0;
+                realtype **jac = block->J_scale->cols;
+                realtype maxInJac = RAbs(jac[0][i]);
+                /* Likely not a problem: solver->using_max_min_scaling_flag = 1; -- Using minimum scaling */
+                for(j = 0; j < N; j++) {
+                    if(RAbs(maxInJac) < RAbs(jac[j][i])) {
+                        maxInJac = jac[j][i];
+                        maxInJacIndex = j;
+                    }
+                }
+                
+                if (scale_ptr[i] > block->residual_heuristic_nominal[i]/bsop->min_residual_scaling_factor) {
+                    scale_ptr[i] = bsop->min_residual_scaling_factor/block->residual_heuristic_nominal[i];
+                    jmi_log_node(block->log, logWarning, "MinScalingUsed", "Poor scaling in <block: %s>. "
+                        "Consider changes in the model. Partial derivative of <equation: %I> with respect to <Iter: #r%d#> is <dRes_dIter: %g> (<dRes_dIter_scaled: %g>).",
+                        block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(block->J, i, maxInJacIndex) ,maxInJac);
+                } else {
+                    scale_ptr[i] = 1/scale_ptr[i];
+                    jmi_log_node(block->log, logWarning, "MinScalingExceeded", "Poor scaling in <block: %s> but will allow it based on structure of block. "
+                        "Consider changes in the model. Partial derivative of <equation: %I> with respect to <Iter: #r%d#> is <dRes_dIter: %g> (<dRes_dIter_scaled: %g>).",
+                        block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(block->J, i, maxInJacIndex) ,maxInJac);
+                }
+            }
+            else
+                scale_ptr[i] = 1/scale_ptr[i];
+        }
+
+        if (block->callbacks->log_options.log_level >= 4) {
+            jmi_log_node_t outer = jmi_log_enter_fmt(block->log, logInfo, "ResidualScalingUpdated", "<block:%s>", block->label);
+            if (block->callbacks->log_options.log_level >= 5) {
+                jmi_log_node_t inner = jmi_log_enter_vector_(block->log, outer, logInfo, "scaling");
+                realtype* res = scale_ptr;
+                for (i=0;i<N;i++) jmi_log_real_(block->log, 1/res[i]);
+                jmi_log_leave(block->log, inner);
+            }
+            jmi_log_leave(block->log, outer);
+        }
+    } else if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) {
+        for(i = 0; i < N; i++) {
+            scale_ptr[i] = 1/scale_ptr[i];
+        }
+        /* Scaling is not updated */
+        if (block->callbacks->log_options.log_level >= 4) {
+            jmi_log_node_t outer = jmi_log_enter_fmt(block->log, logInfo, "ResidualScalingUpdated", "<block:%s>", block->label);
+            if (block->callbacks->log_options.log_level >= 5) {
+                jmi_log_node_t inner = jmi_log_enter_vector_(block->log, outer, logInfo, "scaling");
+                realtype* res = scale_ptr;
+                for (i=0;i<N;i++) jmi_log_real_(block->log, 1/res[i]);
+                jmi_log_leave(block->log, inner);
+            }
+            jmi_log_leave(block->log, outer);
+        }
+    }     
+    return;
+}
+
+void jmi_setup_f_residual_scaling(jmi_block_solver_t *block) {
+    jmi_block_solver_options_t* bsop = block->options;
+    int i, N = block->n;
+    realtype* dummy = 0;
+    
+    /* Read manual scaling from annotations and propagated nominal values for
+     * residuals and put them in residual_nominal & scale_ptr*/
+    if (bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual ||
+        bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid)
+    {
+        block->F(block->problem_data,dummy, block->residual_nominal, JMI_BLOCK_EQUATION_NOMINAL);
+        for (i = 0; i < N; i++) {
+            if(block->residual_nominal[i] != 0.0) {
+                if(block->residual_nominal[i] < 1/bsop->max_residual_scaling_factor) {
+                    block->residual_nominal[i] = 1/bsop->max_residual_scaling_factor;
+                    jmi_log_node(block->log, logWarning, "MaxScalingUsed", "Using maximum scaling factor in <block: %s>, "
+                        "<equation: %I> since specified residual nominal is too small.", block->label, i);
+                }
+                else if(block->residual_nominal[i] > 1/bsop->min_residual_scaling_factor) {
+                    block->residual_nominal[i] = 1/bsop->min_residual_scaling_factor;
+                    jmi_log_node(block->log, logWarning, "MinScalingUsed", "Using minimal scaling factor in <block: %s>, "
+                        "<equation: %I> since specified residual nominal is too big.", block->label, i);
+                }
+            }
+            else if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) {
+                block->residual_nominal[i] = 1.0;
+            }
+        }
+    }
+    
+    block->F(block->problem_data,dummy, block->residual_heuristic_nominal, JMI_BLOCK_EQUATION_NOMINAL_AUTO);
+    block->F(block->problem_data,dummy, block->residual_heuristic_nominal, JMI_BLOCK_EQUATION_NOMINAL);
+    for (i = 0; i < N; i++) {
+        if (block->residual_heuristic_nominal[i] == 0.0) {
+            /* TODO: This information could be used to see that the system is structually singular */
+            block->residual_heuristic_nominal[i] = 1.0;
+        }
+    }
+}
+
+double* jmi_solver_get_f_scales(jmi_block_solver_t* block) {
+    return N_VGetArrayPointer(block->f_scale);
+}
+
+double calculate_scaled_residual_norm(double* residual, double *f_scale, int n) {
+    double scaledRes, norm = 0;
+    int i;
+    for(i=0; i<n; i++) {
+        scaledRes = residual[i]*f_scale[i];
+        norm = norm>= RAbs(scaledRes)?norm:RAbs(scaledRes);
+    }
+    return norm;
+}
+

trunk/RuntimeLibrary/src/jmi/jmi_block_solver.h

@@ -28 +28 @@
 #include "jmi_log.h"
 #include <time.h>
-
+#include <sundials/sundials_math.h>
+#include <sundials/sundials_direct.h>
+#include <nvector/nvector_serial.h>
+#include <kinsol/kinsol_direct.h>
+#include <kinsol/kinsol_impl.h>
+#include <sundials/sundials_dense.h>
 
 /** \brief Evaluation modes for the residual function.*/
@@ -350 +355 @@
 void jmi_block_solver_init_default_options(jmi_block_solver_options_t* op);
 
+/** \brief Update function scaling based on Jacobian information */
+void jmi_update_f_scale(jmi_block_solver_t *block);
+
+/**< \brief Retrieve residual scales used in solver */
+double* jmi_solver_get_f_scales(jmi_block_solver_t* block);
+
+/**< \brief Setup residual scaling */
+void jmi_setup_f_residual_scaling(jmi_block_solver_t *block);
+
+/** \brief Calculate scaled residual max norm */
+double calculate_scaled_residual_norm(double* residual, double *f_scale, int n);
+
 /** \brief Check and log illegal iv inputs */
 int 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

@@ -37 +37 @@
     jmi_string_t label;
 
+    N_Vector f_scale;          /**< \brief Work vector for scaling of f */
+    realtype scale_update_time; /**< \brief The last time when f scale was updated */
     int n;                         /**< \brief The number of iteration variables */
     jmi_real_t* x;                 /**< \brief Work vector for the real iteration variables */
     jmi_real_t* last_accepted_x;   /**< \brief Work vector for the real iteration variables holding the last accepted vales by the integrator */
+    DlsMat J;                       /**< \brief The Jacobian matrix  */
+    DlsMat J_scale;                 /**< \brief Jacobian matrix scaled with xnorm for used for fnorm calculation */
+    int using_max_min_scaling_flag; /**< \brief A flag indicating if either the maximum scaling is used of the minimum */
 
     jmi_real_t* dx;                /**< \brief Work vector for the seed vector */
@@ -66 +71 @@
     jmi_real_t* max;               /**< \brief Max values for iteration variables */
     jmi_real_t* nominal;           /**< \brief Nominal values for iteration variables */
+    jmi_real_t* residual_nominal;   /**< \brief Nominals values for residual variables */
+    jmi_real_t* residual_heuristic_nominal;   /**< \brief Heuristic nominals values for residual variables */
     jmi_real_t* initial;           /**< \brief Initial values for iteration variables */
     jmi_real_t* start_set;         /**< \brief If the start value is specified for the iteration variables */
-    
+
     int jacobian_variability;      /**< \brief Variability of Jacobian coefficients: JMI_CONSTANT_VARIABILITY
                                          JMI_PARAMETER_VARIABILITY, JMI_DISCRETE_VARIABILITY, JMI_CONTINUOUS_VARIABILITY */

trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.c

@@ -48 +48 @@
 static int jmi_kin_lsetup(struct KINMemRec * kin_mem);
 
-/* Update residual scales based on jacobian information */
-static void jmi_update_f_scale(jmi_block_solver_t *block);
+/* Only call directly after jmi_update_f_scale */
+static void jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(jmi_block_solver_t *block);
+
+static realtype jmi_calculate_jacobian_condition_number(jmi_block_solver_t * block);
 
 /* Calculate matrix vector product vv = M*v or vv = M^T*v */
@@ -174 +176 @@
     /* Test if output is OK (no -1.#IND) */
     n = NV_LENGTH_S(ff);
-    ret = jmi_check_and_log_illegal_residual_output(block, f, y, solver->residual_heuristic_nominal ,n);
+    ret = jmi_check_and_log_illegal_residual_output(block, f, y, block->residual_heuristic_nominal ,n);
     
     /* record information for Brent search */
@@ -702 +704 @@
             int i, info, N = block->n;
             realtype tol = solver->kin_stol;
-            realtype *scale_ptr = N_VGetArrayPointer(solver->kin_f_scale);
+            realtype *scale_ptr = N_VGetArrayPointer(block->f_scale);
             if(block->options->residual_equation_scaling_mode != 0) {
                 for(i = 0; i < N; i++){
                     int j;
-                    realtype* scaled_col_ptr = DENSE_COL(solver->J_scale, i);
-                    realtype* col_ptr = DENSE_COL(solver->J, i);
+                    realtype* scaled_col_ptr = DENSE_COL(block->J_scale, i);
+                    realtype* col_ptr = DENSE_COL(block->J, i);
                     realtype xscale = RAbs(block->nominal[i]);
                     realtype x = RAbs(block->x[i]);
@@ -717 +719 @@
                 }
             } else {
-                DenseCopy(solver->J, solver->J_scale);
-            }
-            dgetrf_(  &N, &N, solver->J_scale->data, &N, solver->lapack_iwork, &info);
+                DenseCopy(block->J, block->J_scale);
+            }
+            dgetrf_(  &N, &N, block->J_scale->data, &N, solver->lapack_iwork, &info);
             if(info > 0) {
                 jmi_log_node(block->log, logWarning, "SingularJacobian",
@@ -725 +727 @@
             }
             else {
-                dgecon_(&norm, &N, solver->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
+                dgecon_(&norm, &N, block->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
 
                 if(tol * Jcond < UNIT_ROUNDOFF) {
@@ -809 +811 @@
     jmi_kinsol_solver_t* solver = block->solver;
     struct KINMemRec* kin_mem = solver->kin_mem;
-    realtype* residual_scaling_factors = N_VGetArrayPointer(solver->kin_f_scale);
+    realtype* residual_scaling_factors = N_VGetArrayPointer(block->f_scale);
     jmi_log_t *log = block->log;
     clock_t t = jmi_block_solver_start_clock(block);
@@ -1037 +1039 @@
     jmi_log_fmt(block->log, node, logInfo, " <max_iter_no_jacobian: %d>", block->options->max_iter_no_jacobian);
     jmi_log_fmt(block->log, node, logInfo, " <active_bounds_mode: %d>", block->options->active_bounds_mode);
+    jmi_log_fmt(block->log, node, logInfo, " <start_from_last_integrator_step: %d>", block->options->start_from_last_integrator_step);
     jmi_log_leave(block->log, node);
 
@@ -1490 +1493 @@
     int i,j,k;
     realtype **JTJ_c =  solver->JTJ->cols;
-    realtype **jac = solver->J->cols;
+    realtype **jac = block->J->cols;
     int N = block->n;
     realtype * uscale_data = N_VGetArrayPointer(solver->kin_y_scale);
-    realtype * fscale_data = N_VGetArrayPointer(solver->kin_f_scale);    
+    realtype * fscale_data = N_VGetArrayPointer(block->f_scale);    
 
     for (i=0;i<N;i++) {
@@ -1547 +1550 @@
 }
 
+/* Call this function directly after jmi_update_f_scale */
+static void jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(jmi_block_solver_t *block) {
+
+    jmi_kinsol_solver_t* solver = block->solver;
+    int i, N = block->n;
+    jmi_block_solver_options_t* bsop = block->options;
+    int use_scaling_flag = bsop->residual_equation_scaling_mode;
+    realtype tol = solver->kin_stol;
+    realtype* scale_ptr = N_VGetArrayPointer(block->f_scale);
+
+    if (block->using_max_min_scaling_flag && !solver->J_is_singular_flag) {
+        realtype cond = jmi_calculate_jacobian_condition_number(block);
+
+        jmi_log_node(block->log, logInfo, "Regularization",
+            "Calculated condition number in <block: %s>. Regularizing if <cond: %E> is greater than <regtol: %E>", block->label, cond, solver->kin_reg_tol);
+        if (cond > solver->kin_reg_tol) {
+            if(N > 1 && solver->handling_of_singular_jacobian_flag == JMI_REGULARIZATION) {
+                int info;
+                jmi_kinsol_reg_matrix(block);
+                dgetrf_(  &block->n, &block->n, solver->JTJ->data, &block->n, solver->lapack_ipiv, &info);
+            }
+            solver->J_is_singular_flag = 1;
+        }
+    }
+
+
+    /* estimate condition number of the scaled jacobian 
+    and scale function tolerance with it. */
+    if((N > 1) && bsop->check_jac_cond_flag){
+        realtype* scaled_col_ptr;
+        char norm = 'I';
+        double Jnorm = 1.0, Jcond = 1.0;
+        int info;
+        if(use_scaling_flag) {
+            for(i = 0; i < N; i++){
+                int j;
+                scaled_col_ptr = DENSE_COL(block->J_scale, i);
+                for(j = 0; j < N; j++){
+                    scaled_col_ptr[j] = scaled_col_ptr[j] * scale_ptr[j];
+                }
+            }
+        } else {
+            DenseCopy(block->J, block->J_scale);
+        }
+
+        dgetrf_(  &N, &N, block->J_scale->data, &N, solver->lapack_iwork, &info);
+        if(info > 0) {
+            jmi_log_node(block->log, logWarning, "SingularJacobian",
+                "Singular Jacobian detected when checking condition number in <block:%s>. Solver may fail to converge.", block->label);
+        }
+        else {
+            dgecon_(&norm, &N, block->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
+
+            if(tol * Jcond < UNIT_ROUNDOFF) {
+                jmi_log_node_t node = jmi_log_enter_fmt(block->log, logWarning, "IllConditionedJacobian",
+                    "<JacobianInverseConditionEstimate:%E> Solver may fail to converge in <block: %s>.", Jcond, block->label);
+                if (block->callbacks->log_options.log_level >= 4) {
+                    jmi_log_reals(block->log, node, logWarning, "ivs", N_VGetArrayPointer(solver->kin_y), block->n);
+                }
+                jmi_log_leave(block->log, node);
+
+            }
+            else {
+                jmi_log_node(block->log, logInfo, "JacobianCondition",
+                    "<JacobianInverseConditionEstimate:%E>", Jcond);
+            }
+        }
+    }
+}
 
 /* Estimate condition number utilizing dgecon from LAPACK*/
@@ -1558 +1630 @@
 
     /* Copy Jacobian to factorization matrix */
-    DenseCopy(solver->J, solver->J_LU);
+    DenseCopy(block->J, solver->J_LU);
     /* Perform LU factorization to be used with dgecon */
     dgetrf_(&N, &N, solver->J_LU->data, &N, solver->lapack_ipiv, &info);
@@ -1567 +1639 @@
 
     /* Compute infinity norm of J to be used with dgecon */
-    J_norm = dlange_(&norm, &N, &N, solver->J->data, &N, solver->lapack_work);
+    J_norm = dlange_(&norm, &N, &N, block->J->data, &N, solver->lapack_work);
 
     /* Compute reciprocal condition number */
@@ -1597 +1669 @@
         return 0;
     }
-    SetToZero(solver->J);
+    SetToZero(block->J);
 
     /* Evaluate Jacobian */
-    ret = kin_dF(N, solver->kin_y, kin_mem->kin_fval, solver->J, block, kin_mem->kin_vtemp1, kin_mem->kin_vtemp2);
+    ret = kin_dF(N, solver->kin_y, kin_mem->kin_fval, block->J, block, kin_mem->kin_vtemp1, kin_mem->kin_vtemp2);
     solver->updated_jacobian_flag = 1; /* The Jacobian is current */
     
@@ -1607 +1679 @@
         jmi_log_node(block->log, logWarning, "CompressedJacobianEvaluation", "Failed to evaluate Jacobian using compression. Will try full finite differences.");
         /* Try if it helps with full finite differences with possibility to step in other direction */
-        ret = kin_dF(N, solver->kin_y, kin_mem->kin_fval, solver->J, block, kin_mem->kin_vtemp1, kin_mem->kin_vtemp2);
+        ret = kin_dF(N, solver->kin_y, kin_mem->kin_fval, block->J, block, kin_mem->kin_vtemp1, kin_mem->kin_vtemp2);
         solver->has_compression_setup_flag = TRUE;
     }
@@ -1626 +1698 @@
     }
 
-    if(block->force_rescaling)
+    if(block->force_rescaling) {
         jmi_update_f_scale(block);
+        jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block);
+    }
     
     return 0;
@@ -1646 +1720 @@
     denom = jmi_kinsol_calc_v1twwv2(kin_mem->kin_pp,kin_mem->kin_pp,solver->kin_y_scale);
     /* work_vector = Jac * step */
-    jmi_kinsol_calc_Mv(solver->J, 0, kin_mem->kin_pp, solver->work_vector);
+    jmi_kinsol_calc_Mv(block->J, 0, kin_mem->kin_pp, solver->work_vector);
 
     /* work_vector = (ResidualDelta  - Jac * step) unless below update tolerance */
@@ -1659 +1733 @@
 
     for(j=0; j < N; j++) {
-        realtype *jacCol_j = solver->J->cols[j];
+        realtype *jacCol_j = block->J->cols[j];
         double tempj = Ith(kin_mem->kin_pp, j)*Ith(solver->kin_y_scale,j)*Ith(solver->kin_y_scale,j);
         for(i = 0; i < N; i++) {
@@ -1680 +1754 @@
             jmi_log_fmt(block->log, node,logInfo, "Update denominator <step_norm: %g>", denom);
 #endif
-            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", solver->J->data, N, N);
+            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", block->J->data, N, N);
         }
         jmi_log_leave(block->log, node);
@@ -1699 +1773 @@
 
     /* work_vector = Jac * step */
-    jmi_kinsol_calc_Mv(solver->J, 0, kin_mem->kin_pp, solver->work_vector);
+    jmi_kinsol_calc_Mv(block->J, 0, kin_mem->kin_pp, solver->work_vector);
 
     /* work_vector = (ResidualDelta  - Jac * step) unless below update tolerance */
@@ -1715 +1789 @@
         double denom = 0;
         for(j=0; j < N; j++) {
-            denom += Ith(kin_mem->kin_pp, j)*Ith(kin_mem->kin_pp, j)*solver->J->cols[j][i]*solver->J->cols[j][i];
+            denom += Ith(kin_mem->kin_pp, j)*Ith(kin_mem->kin_pp, j)*block->J->cols[j][i]*block->J->cols[j][i];
         }
         Ith(solver->work_vector2,i) = denom;
@@ -1721 +1795 @@
 
     for(j=0; j < N; j++) {
-        realtype *jacCol_j = solver->J->cols[j];
+        realtype *jacCol_j = block->J->cols[j];
         for(i = 0; i < N; i++) {
             if( Ith(solver->work_vector2,i) >= UNIT_ROUNDOFF) {
@@ -1739 +1813 @@
             jmi_log_reals(block->log, node,logInfo, "res_JacStep", N_VGetArrayPointer(solver->work_vector), N);
 #endif
-            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", solver->J->data, N, N);
+            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", block->J->data, N, N);
         }
         jmi_log_leave(block->log, node);
@@ -1759 +1833 @@
     /* work_vector2 = res_diff */
     N_VLinearSum(1, solver->last_residual, -1, b, solver->work_vector);
-    N_VProd(solver->work_vector, solver->kin_f_scale, solver->work_vector);
-    N_VProd(solver->work_vector, solver->kin_f_scale, solver->work_vector);
+    N_VProd(solver->work_vector, block->f_scale, solver->work_vector);
+    N_VProd(solver->work_vector, block->f_scale, solver->work_vector);
     /* work_vector2 = res_diff^T*f_scale^2*Jac */
-    jmi_kinsol_calc_Mv(solver->J, TRUE, solver->work_vector, solver->work_vector2);
+    jmi_kinsol_calc_Mv(block->J, TRUE, solver->work_vector, solver->work_vector2);
 
     /* denom1 = step^T.*kin_y_scale^2*step */
     denom1 = jmi_kinsol_calc_v1twwv2(kin_mem->kin_pp,kin_mem->kin_pp,solver->kin_y_scale);
     /* work_vector = Jac * step */
-    jmi_kinsol_calc_Mv(solver->J, 0, kin_mem->kin_pp, solver->work_vector);
+    jmi_kinsol_calc_Mv(block->J, 0, kin_mem->kin_pp, solver->work_vector);
     /* denom2 = res_diff^T.*kin_f_scale^2*Jac*step */
     denom2 = 0;
     for(i = 0; i < N; i++) {
-        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);
+        denom2 += -(Ith(b,i) - Ith(solver->last_residual,i))*Ith(block->f_scale,i)*Ith(block->f_scale,i)*Ith(solver->work_vector,i);
     }
     
     for(j=0; j < N; j++) {
-        realtype *jacCol_j = solver->J->cols[j];
+        realtype *jacCol_j = block->J->cols[j];
         double tempj1 = Ith(kin_mem->kin_pp,j)*Ith(solver->kin_y_scale,j)*Ith(solver->kin_y_scale,j)/denom1;
         double tempj2 = Ith(solver->work_vector2,j)/denom2;
@@ -1797 +1871 @@
             jmi_log_fmt(block->log, node,logInfo, "Update denominator <denom2: %g>", denom2);
 #endif
-            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", solver->J->data, N, N);
+            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", block->J->data, N, N);
         }
         jmi_log_leave(block->log, node);
@@ -1815 +1889 @@
     clock_t t = jmi_block_solver_start_clock(block);
       
-    DenseCopy(solver->J, solver->J_LU); /* make a copy of the Jacobian that will be used for LU factorization */
+    DenseCopy(block->J, solver->J_LU); /* make a copy of the Jacobian that will be used for LU factorization */
 
     /* Equillibrate if corresponding option is set */
@@ -1875 +1949 @@
                              "Will try to find the minimum norm solution in <block: %s>", block->label);
                 SetToZero(solver->J_sing);
-                DenseCopy(solver->J, solver->J_sing);
+                DenseCopy(block->J, solver->J_sing);
             } else {
                 /* Error */
@@ -1950 +2024 @@
         int i;
         jmi_log_node_t node = jmi_log_enter_fmt(block->log, logInfo, "AggressiveResidualScalingUpdate", "Updating f_scale aggressively");
-        N_VScale(1.0, solver->kin_f_scale, solver->work_vector);
+        N_VScale(1.0, block->f_scale, solver->work_vector);
         jmi_update_f_scale(block);
+        jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block);
         for(i=0; i<block->n; i++) {
-            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);
+            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);
         }
         if (block->callbacks->log_options.log_level >= 4) {
@@ -1959 +2034 @@
             if (block->callbacks->log_options.log_level >= 5) {
                 jmi_log_node_t inner = jmi_log_enter_vector_(block->log, outer, logInfo, "scaling");
-                realtype* res = N_VGetArrayPointer(solver->kin_f_scale);
+                realtype* res = N_VGetArrayPointer(block->f_scale);
                 for (i=0;i<N;i++) jmi_log_real_(block->log, 1/res[i]);
                 jmi_log_leave(block->log, inner);
@@ -1993 +2068 @@
     }
   
-    kin_mem->kin_sJpnorm = N_VWL2Norm(b,solver->kin_f_scale);
+    kin_mem->kin_sJpnorm = N_VWL2Norm(b,block->f_scale);
     solver->sJpnorm = kin_mem->kin_sJpnorm;
-    N_VProd(b, solver->kin_f_scale, x);
-    N_VProd(x, solver->kin_f_scale, x);
+    N_VProd(b, block->f_scale, x);
+    N_VProd(x, block->f_scale, x);
     
     kin_mem->kin_sfdotJp = N_VDotProd(kin_mem->kin_fval, x);
@@ -2014 +2089 @@
          */
 
-        jmi_kinsol_calc_Mv(solver->J, TRUE, x, solver->gradient);
+        jmi_kinsol_calc_Mv(block->J, TRUE, x, solver->gradient);
     }
     if(solver->use_steepest_descent_flag) {
@@ -2023 +2098 @@
     }
     else if(solver->J_is_singular_flag) {
-        realtype** jac = solver->J->cols;
+        realtype** jac = block->J->cols;
         if (N == 1) {
             xd[0] = block->nominal[0] * 0.1 *((bd[0] > 0)?1:-1) * ((jac[0][0] > 0)?1:-1);
@@ -2030 +2105 @@
             /* solve the regularized problem */
             int i,j;
-            realtype * fscale_data = N_VGetArrayPointer(solver->kin_f_scale);
+            realtype * fscale_data = N_VGetArrayPointer(block->f_scale);
             realtype gnorm;/*gradient norm*/
 
@@ -2120 +2195 @@
         
         if((block->callbacks->log_options.log_level >= 6)) {
-            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", solver->J->data, N, N);
+            jmi_log_real_matrix(block->log, node, logInfo, "jacobian", block->J->data, N, N);
             jmi_log_reals(block->log, node, logInfo, "rhs", xd, N);
         }
@@ -2203 +2278 @@
             
             /* recalculate sJpnorm  */ 
-            jmi_kinsol_calc_Mv(solver->J, FALSE, x, b);
-            sJpnorm = N_VWL2Norm(b,solver->kin_f_scale);
+            jmi_kinsol_calc_Mv(block->J, FALSE, x, b);
+            sJpnorm = N_VWL2Norm(b,block->f_scale);
             /* update sJpnorm and sfdotJp for Kinsol */ 
             kin_mem->kin_sJpnorm = sJpnorm;
@@ -2219 +2294 @@
     block->step_calc_time += jmi_block_solver_elapsed_time(block, t);
     return 0;
-}
-
-double* jmi_kinsol_solver_get_f_scales(jmi_kinsol_solver_t* solver) {
-    return N_VGetArrayPointer(solver->kin_f_scale);
-}
-
-static void jmi_setup_f_residual_scaling(jmi_block_solver_t *block) {
-    jmi_block_solver_options_t* bsop = block->options;
-    jmi_kinsol_solver_t* solver = block->solver;
-    int i, N = block->n;
-    realtype* dummy = 0;
-    
-    /* Read manual scaling from annotations and propagated nominal values for
-     * residuals and put them in residual_nominal & scale_ptr*/
-    if (bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual ||
-        bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid)
-    {
-        block->F(block->problem_data,dummy, solver->residual_nominal, JMI_BLOCK_EQUATION_NOMINAL);
-        for (i = 0; i < N; i++) {
-            if(solver->residual_nominal[i] != 0.0) {
-                if(solver->residual_nominal[i] < 1/bsop->max_residual_scaling_factor) {
-                    solver->residual_nominal[i] = 1/bsop->max_residual_scaling_factor;
-                    jmi_log_node(block->log, logWarning, "MaxScalingUsed", "Using maximum scaling factor in <block: %s>, "
-                        "<equation: %I> since specified residual nominal is too small.", block->label, i);
-                }
-                else if(solver->residual_nominal[i] > 1/bsop->min_residual_scaling_factor) {
-                    solver->residual_nominal[i] = 1/bsop->min_residual_scaling_factor;
-                    jmi_log_node(block->log, logWarning, "MinScalingUsed", "Using minimal scaling factor in <block: %s>, "
-                        "<equation: %I> since specified residual nominal is too big.", block->label, i);
-                }
-            }
-            else if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) {
-                solver->residual_nominal[i] = 1.0;
-            }
-        }
-    }
-    
-    block->F(block->problem_data,dummy, solver->residual_heuristic_nominal, JMI_BLOCK_EQUATION_NOMINAL_AUTO);
-    block->F(block->problem_data,dummy, solver->residual_heuristic_nominal, JMI_BLOCK_EQUATION_NOMINAL);
-    for (i = 0; i < N; i++) {
-        if (solver->residual_heuristic_nominal[i] == 0.0) {
-            /* TODO: This information could be used to see that the system is structually singular */
-            solver->residual_heuristic_nominal[i] = 1.0;
-        }
-    }
-}
-
-/* 
-    Compute appropriate equation scaling and function tolerance based on Jacobian J,
-    nominal values (block->nominal) and current point (block->x).
-    Store result in solver->kin_f_scale.
-*/
-static void jmi_update_f_scale(jmi_block_solver_t *block) {
-    jmi_kinsol_solver_t* solver = block->solver; 
-    int i, N = block->n;
-    realtype tol = solver->kin_stol;
-    realtype curtime = block->cur_time;
-    realtype* scale_ptr = 0;
-    jmi_block_solver_options_t* bsop = block->options;
-    int use_scaling_flag = bsop->residual_equation_scaling_mode;
-
-
-    solver->kin_scale_update_time = curtime;
-    block->force_rescaling = 0;
-    
-    if(bsop->residual_equation_scaling_mode != jmi_residual_scaling_none) {
-        /* Zero out the scales initially if we're modify this. */
-        N_VConst_Serial(0,solver->kin_f_scale);
-    }
-
-    scale_ptr = N_VGetArrayPointer(solver->kin_f_scale);
-
-    /* Form scaled Jacobian as needed for automatic scaling and condition number checking*/
-    if (bsop->residual_equation_scaling_mode == jmi_residual_scaling_auto || 
-        bsop->residual_equation_scaling_mode == jmi_residual_scaling_aggressive_auto ||
-        bsop->residual_equation_scaling_mode == jmi_residual_scaling_full_jacobian_auto ||
-        bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid ||
-        bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual)
-    {
-        kin_char_log(solver, 's');
-        for (i = 0; i < N; i++) {
-            int j;
-            /* column scaling is formed by max(abs(nominal), abs(actual_value)) */
-            realtype xscale = RAbs(block->nominal[i]);
-            realtype x = RAbs(block->x[i]);
-            realtype* col_ptr;
-            realtype* scaled_col_ptr;
-            if(x < xscale) x = xscale;
-            if(x < tol) x = tol;
-            col_ptr = DENSE_COL(solver->J, i);
-            scaled_col_ptr = DENSE_COL(solver->J_scale, i);
-
-            /* row scaling is product of Jac entry and column scaling */
-            for (j = 0; j < N; j++) {
-                realtype dres = col_ptr[j];
-                realtype fscale;
-                fscale = dres * x;
-                scaled_col_ptr[j] = fscale;
-                scale_ptr[j] = MAX(scale_ptr[j], RAbs(fscale));
-            }
-        }
-    }
-
-    /* put in manual non-zero scales in hybrid & manual cases */
-    if(    (bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) 
-        || (bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid))    {
-            for(i = 0; i < N; i++) {
-                if(solver->residual_nominal[i] != 0) {
-                    scale_ptr[i] = solver->residual_nominal[i];
-                }
-            }
-    }
-
-    if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_auto 
-        || bsop->residual_equation_scaling_mode == jmi_residual_scaling_aggressive_auto
-        || bsop->residual_equation_scaling_mode == jmi_residual_scaling_full_jacobian_auto
-        || bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid) {
-        solver->using_max_min_scaling_flag = 0; /* NOT using max/min scaling */
-        /* check that scaling factors have reasonable magnitude */
-        for(i = 0; i < N; i++) {
-            if(scale_ptr[i] < 1/bsop->max_residual_scaling_factor) {
-                int j, maxInJacIndex = 0;
-                realtype **jac = solver->J_scale->cols;
-                realtype maxInJac = RAbs(jac[0][i]); 
-                solver->using_max_min_scaling_flag = 1; /* Using maximum scaling, singular Jacobian? */
-                for(j = 0; j < N; j++) {
-                    if(RAbs(maxInJac) < RAbs(jac[j][i])) {
-                        maxInJac = jac[j][i];
-                        maxInJacIndex = j;
-                    }
-                }
-                
-                scale_ptr[i] = bsop->max_residual_scaling_factor;
-                jmi_log_node(block->log, logWarning, "MaxScalingUsed", "Poor scaling in <block: %s>. "
-                    "Consider changes in the model. Partial derivative of <equation: %I> with respect to <Iter: #r%d#> is <dRes_dIter: %g> (<dRes_dIter_scaled: %g>).",
-                    block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(solver->J, i, maxInJacIndex) ,maxInJac);
-            }
-            else if(scale_ptr[i] > 1/bsop->min_residual_scaling_factor) {
-                int j, maxInJacIndex = 0;
-                realtype **jac = solver->J_scale->cols;
-                realtype maxInJac = RAbs(jac[0][i]);
-                /* Likely not a problem: solver->using_max_min_scaling_flag = 1; -- Using minimum scaling */
-                for(j = 0; j < N; j++) {
-                    if(RAbs(maxInJac) < RAbs(jac[j][i])) {
-                        maxInJac = jac[j][i];
-                        maxInJacIndex = j;
-                    }
-                }
-                
-                if (scale_ptr[i] > solver->residual_heuristic_nominal[i]/bsop->min_residual_scaling_factor) {
-                    scale_ptr[i] = bsop->min_residual_scaling_factor/solver->residual_heuristic_nominal[i];
-                    jmi_log_node(block->log, logWarning, "MinScalingUsed", "Poor scaling in <block: %s>. "
-                        "Consider changes in the model. Partial derivative of <equation: %I> with respect to <Iter: #r%d#> is <dRes_dIter: %g> (<dRes_dIter_scaled: %g>).",
-                        block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(solver->J, i, maxInJacIndex) ,maxInJac);
-                } else {
-                    scale_ptr[i] = 1/scale_ptr[i];
-                    jmi_log_node(block->log, logWarning, "MinScalingExceeded", "Poor scaling in <block: %s> but will allow it based on structure of block. "
-                        "Consider changes in the model. Partial derivative of <equation: %I> with respect to <Iter: #r%d#> is <dRes_dIter: %g> (<dRes_dIter_scaled: %g>).",
-                        block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(solver->J, i, maxInJacIndex) ,maxInJac);
-                }
-            }
-            else
-                scale_ptr[i] = 1/scale_ptr[i];
-        }
-
-        if (block->callbacks->log_options.log_level >= 4) {
-            jmi_log_node_t outer = jmi_log_enter_fmt(block->log, logInfo, "ResidualScalingUpdated", "<block:%s>", block->label);
-            if (block->callbacks->log_options.log_level >= 5) {
-                jmi_log_node_t inner = jmi_log_enter_vector_(block->log, outer, logInfo, "scaling");
-                realtype* res = scale_ptr;
-                for (i=0;i<N;i++) jmi_log_real_(block->log, 1/res[i]);
-                jmi_log_leave(block->log, inner);
-            }
-            jmi_log_leave(block->log, outer);
-        }
-    } else if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) {
-        for(i = 0; i < N; i++) {
-            scale_ptr[i] = 1/scale_ptr[i];
-        }
-        /* Scaling is not updated */
-        if (block->callbacks->log_options.log_level >= 4) {
-            jmi_log_node_t outer = jmi_log_enter_fmt(block->log, logInfo, "ResidualScalingUpdated", "<block:%s>", block->label);
-            if (block->callbacks->log_options.log_level >= 5) {
-                jmi_log_node_t inner = jmi_log_enter_vector_(block->log, outer, logInfo, "scaling");
-                realtype* res = scale_ptr;
-                for (i=0;i<N;i++) jmi_log_real_(block->log, 1/res[i]);
-                jmi_log_leave(block->log, inner);
-            }
-            jmi_log_leave(block->log, outer);
-        }
-    } 
-
-    if (solver->using_max_min_scaling_flag && !solver->J_is_singular_flag) {
-        realtype cond = jmi_calculate_jacobian_condition_number(block);
-
-        jmi_log_node(block->log, logInfo, "Regularization",
-            "Calculated condition number in <block: %s>. Regularizing if <cond: %E> is greater than <regtol: %E>", block->label, cond, solver->kin_reg_tol);
-        if (cond > solver->kin_reg_tol) {
-            if(N > 1 && solver->handling_of_singular_jacobian_flag == JMI_REGULARIZATION) {
-                int info;
-                jmi_kinsol_reg_matrix(block);
-                dgetrf_(  &block->n, &block->n, solver->JTJ->data, &block->n, solver->lapack_ipiv, &info);
-            }
-            solver->J_is_singular_flag = 1;
-        }
-    }
-
-
-    /* estimate condition number of the scaled jacobian 
-    and scale function tolerance with it. */
-    if((N > 1) && bsop->check_jac_cond_flag){
-        realtype* scaled_col_ptr;
-        char norm = 'I';
-        double Jnorm = 1.0, Jcond = 1.0;
-        int info;
-        if(use_scaling_flag) {
-            for(i = 0; i < N; i++){
-                int j;
-                scaled_col_ptr = DENSE_COL(solver->J_scale, i);
-                for(j = 0; j < N; j++){
-                    scaled_col_ptr[j] = scaled_col_ptr[j] * scale_ptr[j];
-                }
-            }
-        } else {
-            DenseCopy(solver->J, solver->J_scale);
-        }
-
-        dgetrf_(  &N, &N, solver->J_scale->data, &N, solver->lapack_iwork, &info);
-        if(info > 0) {
-            jmi_log_node(block->log, logWarning, "SingularJacobian",
-                "Singular Jacobian detected when checking condition number in <block:%s>. Solver may fail to converge.", block->label);
-        }
-        else {
-            dgecon_(&norm, &N, solver->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
-
-            if(tol * Jcond < UNIT_ROUNDOFF) {
-                jmi_log_node_t node = jmi_log_enter_fmt(block->log, logWarning, "IllConditionedJacobian",
-                    "<JacobianInverseConditionEstimate:%E> Solver may fail to converge in <block: %s>.", Jcond, block->label);
-                if (block->callbacks->log_options.log_level >= 4) {
-                    jmi_log_reals(block->log, node, logWarning, "ivs", N_VGetArrayPointer(solver->kin_y), block->n);
-                }
-                jmi_log_leave(block->log, node);
-
-            }
-            else {
-                jmi_log_node(block->log, logInfo, "JacobianCondition",
-                    "<JacobianInverseConditionEstimate:%E>", Jcond);
-            }
-        }
-    }
-    return;
 }
 
@@ -2493 +2317 @@
     solver->kin_y = N_VMake_Serial(n, block->x);
     solver->kin_y_scale = N_VNew_Serial(n);
-    solver->kin_f_scale = N_VNew_Serial(n);
     solver->gradient  = N_VNew_Serial(n);    
     solver->last_residual = N_VNew_Serial(n);
-    solver->residual_nominal = (realtype*)calloc(n+1,sizeof(realtype));
-    solver->residual_heuristic_nominal = (realtype*)calloc(n+1,sizeof(realtype));
-    solver->kin_scale_update_time = -1.0;
     solver->kin_jac_update_time = -1.0;
     /*NOTE: it'd be nice to use "jmi->newton_tolerance" here
@@ -2505 +2325 @@
     solver->kin_ftol = UROUND; /* block->options->min_tol; */
     solver->kin_stol = block->options->min_tol;
-    solver->using_max_min_scaling_flag = 0; /* Not using max/min scaling */
     solver->has_compression_setup_flag = FALSE;
     solver->is_first_newton_solve_flag = TRUE;
     solver->current_nni = 0;
     
-    solver->J = NewDenseMat(n ,n);
     solver->JTJ = NewDenseMat(n ,n);
     solver->J_LU = NewDenseMat(n ,n);
-    solver->J_scale = NewDenseMat(n ,n);
     solver->J_sing = NewDenseMat(n, n);
     solver->J_SVD_U = NewDenseMat(n, n);
@@ -2549 +2366 @@
     /* Initialize scaling to 1.0 - defaults */
     N_VConst_Serial(1.0,solver->kin_y_scale);
-
-    /* Initial equation scaling is 1.0 */
-    N_VConst_Serial(1.0,solver->kin_f_scale);
                 
     flag = KINInit(solver->kin_mem, kin_f, solver->kin_y); /*Initialize Kinsol*/
@@ -2637 +2451 @@
     N_VDestroy_Serial(solver->kin_y);
     N_VDestroy_Serial(solver->kin_y_scale);
-    N_VDestroy_Serial(solver->kin_f_scale);
     N_VDestroy_Serial(solver->gradient);
     N_VDestroy_Serial(solver->last_residual);
-    free(solver->residual_nominal);
-    free(solver->residual_heuristic_nominal);
-    DestroyMat(solver->J);
     DestroyMat(solver->JTJ);
     DestroyMat(solver->J_LU);
-    DestroyMat(solver->J_scale);
     DestroyMat(solver->J_sing);
     DestroyMat(solver->J_SVD_U);
@@ -2771 +2580 @@
     } 
     jmi_kinsol_solver_print_solve_start(block, &topnode);
-    flag = KINSol(solver->kin_mem, solver->kin_y, strategy, solver->kin_y_scale, solver->kin_f_scale);
+    flag = KINSol(solver->kin_mem, solver->kin_y, strategy, solver->kin_y_scale, block->f_scale);
     if(flag == KIN_INITIAL_GUESS_OK) { /* last_fnorm not up to date in this case */
         double max_residual = 0.0;
@@ -2777 +2586 @@
         realtype* f;
         char msg[256];
-        realtype* residual_scaling_factors = NV_DATA_S(solver->kin_f_scale);
+        realtype* residual_scaling_factors = NV_DATA_S(block->f_scale);
 
         /*kin_f(solver->kin_y, solver->work_vector, block); */
@@ -2794 +2603 @@
         solver->last_max_residual= max_residual;
         solver->last_max_residual_index = max_index;
-        solver->last_fnorm = N_VWL2Norm(kin_mem->kin_fval, solver->kin_f_scale);
+        solver->last_fnorm = N_VWL2Norm(kin_mem->kin_fval, block->f_scale);
         kin_mem->kin_fnorm = solver->last_fnorm;
         sprintf(msg, "nni = %4ld   nfe = %6ld   fnorm = %26.16g", kin_mem->kin_nni, kin_mem->kin_nfe, solver->last_fnorm);
@@ -2810 +2619 @@
     } else if (flag == KIN_LINESEARCH_NONCONV || flag == KIN_STEP_LT_STPTOL) {
         realtype fnorm;
-        N_VProd(solver->kin_f_scale, kin_mem->kin_fval, solver->work_vector);
+        N_VProd(block->f_scale, kin_mem->kin_fval, solver->work_vector);
         fnorm =N_VMaxNorm(solver->work_vector);
         if((fnorm <= solver->kin_stol && !block->options->solver_exit_criterion_mode == jmi_exit_criterion_step_residual)
@@ -2826 +2635 @@
         if(block->options->check_jac_cond_flag) {
             int i, info, N = block->n;
-            realtype* scale_ptr = N_VGetArrayPointer(solver->kin_f_scale);
+            realtype* scale_ptr = N_VGetArrayPointer(block->f_scale);
             realtype tol = solver->kin_stol;
             if(block->options->residual_equation_scaling_mode !=0) {
                 for(i = 0; i < N; i++){
                     int j;
-                    realtype* scaled_col_ptr = DENSE_COL(solver->J_scale, i);
-                    realtype* col_ptr = DENSE_COL(solver->J, i);
+                    realtype* scaled_col_ptr = DENSE_COL(block->J_scale, i);
+                    realtype* col_ptr = DENSE_COL(block->J, i);
                     realtype xscale = RAbs(block->nominal[i]);
                     realtype x = RAbs(block->x[i]);
@@ -2842 +2651 @@
                 }
             } else {
-                DenseCopy(solver->J, solver->J_scale);
-            }
-
-            dgetrf_(  &N, &N, solver->J_scale->data, &N, solver->lapack_iwork, &info);
+                DenseCopy(block->J, block->J_scale);
+            }
+
+            dgetrf_(  &N, &N, block->J_scale->data, &N, solver->lapack_iwork, &info);
             if(info > 0) {
                 jmi_log_node(block->log, logWarning, "SingularJacobian",
@@ -2853 +2662 @@
                 char norm = 'I';
                 double Jnorm = 1.0, Jcond = 1.0;
-                dgecon_(&norm, &N, solver->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
+                dgecon_(&norm, &N, block->J_scale->data, &N, &Jnorm, &Jcond, solver->lapack_work, solver->lapack_iwork,&info);       
 
                 if(tol * Jcond < UNIT_ROUNDOFF) {
@@ -2893 +2702 @@
     
     if(block->init) {
-        jmi_setup_f_residual_scaling(block);
         flag = jmi_kinsol_init(block);
         if(flag) return flag;
@@ -2925 +2733 @@
 
     /* update the scaling only once per time step */
-    if(block->init || (block->options->rescale_each_step_flag && (curtime > solver->kin_scale_update_time)) || block->force_rescaling) {
+    if(block->init || (block->options->rescale_each_step_flag && (curtime > block->scale_update_time)) || block->force_rescaling) {
         jmi_update_f_scale(block);
+        jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block);
     }
     
@@ -2971 +2780 @@
         solver->use_steepest_descent_flag = 1;
         
-        flag = KINSol(solver->kin_mem, solver->kin_y, KIN_LINESEARCH, solver->kin_y_scale, solver->kin_f_scale);
+        flag = KINSol(solver->kin_mem, solver->kin_y, KIN_LINESEARCH, solver->kin_y_scale, block->f_scale);
         if(flag == KIN_INITIAL_GUESS_OK) {
             flag = KIN_SUCCESS;
@@ -2998 +2807 @@
         /* Update the scaling  */
         jmi_update_f_scale(block);
+        jmi_regularize_and_do_condition_estimate_on_scaled_jacobian(block);
         
         /* For the second solve set tight tolerance on the step and specified tolerance on the residual */
@@ -3121 +2931 @@
             DenseCopy(solver->saved_state->J_modified, solver->J_LU);
     }
-    DenseCopy(solver->saved_state->J, solver->J);
-    
-    solver->kin_scale_update_time = solver->saved_state->kin_scale_update_time;
+    DenseCopy(solver->saved_state->J, block->J);
+    
+    block->scale_update_time = solver->saved_state->kin_scale_update_time;
     block->force_rescaling        = solver->saved_state->force_rescaling;
 
@@ -3130 +2940 @@
     solver->J_is_singular_flag = solver->saved_state->J_is_singular_flag;
     ((struct KINMemRec *)solver->kin_mem)->kin_nnilset = ((struct KINMemRec *)solver->kin_mem)->kin_nni - solver->saved_state->mbset;
-    memcpy(N_VGetArrayPointer(solver->kin_f_scale), N_VGetArrayPointer(solver->saved_state->kin_f_scale), block->n*sizeof(realtype));
+    memcpy(N_VGetArrayPointer(block->f_scale), N_VGetArrayPointer(solver->saved_state->kin_f_scale), block->n*sizeof(realtype));
     memcpy(N_VGetArrayPointer(solver->kin_y_scale), N_VGetArrayPointer(solver->saved_state->kin_y_scale), block->n*sizeof(realtype));
     memcpy(solver->lapack_ipiv, solver->saved_state->lapack_ipiv, (block->n+2)*sizeof(int)); 
@@ -3139 +2949 @@
         jmi_log_ints(block->log, node, logInfo, "is_singular", &(block->force_rescaling), 1);
         jmi_log_ints(block->log, node, logInfo, "force_rescaling", &(solver->J_is_singular_flag), 1);
-        jmi_log_reals(block->log, node, logInfo, "residual_scaling_update_time", &(solver->kin_scale_update_time), 1);
+        jmi_log_reals(block->log, node, logInfo, "residual_scaling_update_time", &(block->scale_update_time), 1);
         jmi_log_reals(block->log, node, logInfo, "iv_scaling", N_VGetArrayPointer(solver->kin_y_scale), block->n);
-        jmi_log_reals(block->log, node, logInfo, "residual_scaling", N_VGetArrayPointer(solver->kin_f_scale), block->n);
+        jmi_log_reals(block->log, node, logInfo, "residual_scaling", N_VGetArrayPointer(block->f_scale), block->n);
         jmi_log_leave(block->log, node);
     }
@@ -3180 +2990 @@
             DenseCopy(solver->J_LU, solver->saved_state->J_modified);
         }
-        DenseCopy(solver->J, solver->saved_state->J);
+        DenseCopy(block->J, solver->saved_state->J);
         
-        solver->saved_state->kin_scale_update_time = solver->kin_scale_update_time;
+        solver->saved_state->kin_scale_update_time = block->scale_update_time;
         solver->saved_state->force_rescaling       = block->force_rescaling;
         solver->saved_state->force_new_J_flag      = solver->force_new_J_flag;
@@ -3189 +2999 @@
         solver->saved_state->J_is_singular_flag = solver->J_is_singular_flag;
         solver->saved_state->mbset = ((struct KINMemRec *)solver->kin_mem)->kin_nni - ((struct KINMemRec *)solver->kin_mem)->kin_nnilset;
-        memcpy(N_VGetArrayPointer(solver->saved_state->kin_f_scale),N_VGetArrayPointer(solver->kin_f_scale), block->n*sizeof(realtype));
+        memcpy(N_VGetArrayPointer(solver->saved_state->kin_f_scale),N_VGetArrayPointer(block->f_scale), block->n*sizeof(realtype));
         memcpy(N_VGetArrayPointer(solver->saved_state->kin_y_scale),N_VGetArrayPointer(solver->kin_y_scale), block->n*sizeof(realtype));
         memcpy(solver->saved_state->lapack_ipiv, solver->lapack_ipiv, (block->n+2)*sizeof(int)); 
@@ -3198 +3008 @@
             jmi_log_ints(block->log, node, logInfo, "is_singular", &(solver->J_is_singular_flag), 1);
             jmi_log_ints(block->log, node, logInfo, "force_rescaling", &(solver->J_is_singular_flag), 1);
-            jmi_log_reals(block->log, node, logInfo, "residual_scaling_update_time", &(solver->kin_scale_update_time), 1);
+            jmi_log_reals(block->log, node, logInfo, "residual_scaling_update_time", &(block->scale_update_time), 1);
             jmi_log_reals(block->log, node, logInfo, "iv_scaling", N_VGetArrayPointer(solver->kin_y_scale), block->n);
-            jmi_log_reals(block->log, node, logInfo, "residual_scaling", N_VGetArrayPointer(solver->kin_f_scale), block->n);
+            jmi_log_reals(block->log, node, logInfo, "residual_scaling", N_VGetArrayPointer(block->f_scale), block->n);
             jmi_log_leave(block->log, node);
         }
@@ -3206 +3016 @@
     return 0;
 }
-

trunk/RuntimeLibrary/src/jmi/jmi_kinsol_solver.h

@@ -58 +58 @@
 int jmi_kinsol_restore_state(jmi_block_solver_t* block);
 
-/**< \brief Retrieve residual scales used in Kinsol solver */
-double* jmi_kinsol_solver_get_f_scales(jmi_kinsol_solver_t* solver);
-
 /**< \brief Convert Kinsol return flag to readable name */
 const char *jmi_kinsol_flag_to_name(int flag);
@@ -87 +84 @@
 
     N_Vector kin_y_scale;          /**< \brief Work vector for Kinsol scaling of y */
-    N_Vector kin_f_scale;          /**< \brief Work vector for Kinsol scaling of f */
     N_Vector gradient;              /**< \brief Steepest descent direction */
-    realtype* residual_nominal;      /**< \brief Vector for reading in manual scaling factors for f  */
-    realtype* residual_heuristic_nominal; /**< \brief Vector for reading in heuristic scaling factors for f  */
-    realtype kin_scale_update_time; /**< \brief The last time when Kinsol scale was updated */
     realtype kin_jac_update_time; /**< \brief The last time when Jacobian was updated */
     realtype kin_ftol;             /**< \brief Tolerance for F */
@@ -97 +90 @@
     realtype kin_reg_tol;          /**< \brief Regularization tolerance */
     
-    DlsMat J;                       /**< \brief The Jacobian matrix  */    
     DlsMat JTJ;                     /**< \brief The Transpose(J).J used if J is singular */
     int J_is_singular_flag;         /**< \brief A flag indicating that J is singular. Regularized JTJ is setup */
     int use_steepest_descent_flag;  /**< \brief A flag indicating that steepest descent and not Newton direction should be used */
     int force_new_J_flag;           /**< \brief A flag indicating that J needs to be recalculated */
-    int using_max_min_scaling_flag; /**< \brief A flag indicating if either the maximum scaling is used of the minimum */
     int updated_jacobian_flag;      /**< \brief A flag indicating if an updated Jacobian is used to solve the system */
     int handling_of_singular_jacobian_flag; /**< \brief A flag for determining how singular systems should be treated */
     DlsMat J_LU;                    /**< \brief Jacobian matrix/it's LU decomposition */
-    DlsMat J_scale;                 /**< \brief Jacobian matrix scaled with xnorm for used for fnorm calculation */
     DlsMat J_sing;                  /**< \brief Jacobian matrix/it's right singular vectors */
     DlsMat J_SVD_U;                 /**< \brief The left singular vectors */

trunk/RuntimeLibrary/src/jmi/jmi_linear_solver.c

@@ -40 +40 @@
     /* Initialize work vectors.*/
     solver->factorization = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t));
-    solver->jacobian = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t));
     solver->dependent_set = (jmi_real_t*)calloc(n_x*n_x,sizeof(jmi_real_t));
     solver->jacobian_temp = (jmi_real_t*)calloc(2*n_x*n_x,sizeof(jmi_real_t));
@@ -82 +81 @@
     for (i = 0; i < n_x; i++) { solver->dependent_set[(n_x-1)*n_x+i] = 0; }
     
-    memcpy(solver->jacobian_temp, solver->jacobian, n_x*n_x*sizeof(jmi_real_t));
+    memcpy(solver->jacobian_temp, block->J->data, n_x*n_x*sizeof(jmi_real_t));
     /* Compute the null space of A */
     /*
@@ -197 +196 @@
     if (solver->n_extra_rows > 0) {
         for (j = 0; j < n_x; j++) {
-            memcpy(&solver->jacobian_temp[j*(n_x+solver->n_extra_rows)], &solver->jacobian[j*n_x], n_x*sizeof(jmi_real_t));
+            memcpy(&solver->jacobian_temp[j*(n_x+solver->n_extra_rows)], &block->J->data[j*n_x], n_x*sizeof(jmi_real_t));
             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));
         }
@@ -246 +245 @@
              A regularization strategy for simple cases singular jac should be introduced.
           */
-        info = block->F(block->problem_data,NULL,solver->jacobian,JMI_BLOCK_EVALUATE_JACOBIAN);
-        memcpy(solver->factorization, solver->jacobian, n_x*n_x*sizeof(jmi_real_t));
+        info = block->F(block->problem_data,NULL,block->J->data,JMI_BLOCK_EVALUATE_JACOBIAN);
+        memcpy(solver->factorization, block->J->data, n_x*n_x*sizeof(jmi_real_t));
         if(info) {
             if(block->init) {
@@ -275 +274 @@
             for (j = 0; j < n_x; j++) {
                 /* Unrecoverable error*/
-                if ( solver->jacobian[i*n_x+j] - solver->jacobian[i*n_x+j] != 0) {
+                if ( block->J->data[i*n_x+j] - block->J->data[i*n_x+j] != 0) {
                     jmi_log_node(block->log, logError, "NaNOutput", "Not a number in the Jacobian <row: %I> <col: %I> from <block: %s>", 
                             i,j, block->label);
@@ -289 +288 @@
                                      "Linear solver invoked for <block:%s>", block->label);
         jmi_log_reals(block->log, destnode, logInfo, "ivs", block->x, block->n);
-        jmi_log_real_matrix(block->log, destnode, logInfo, "A", solver->jacobian, block->n, block->n);
+        jmi_log_real_matrix(block->log, destnode, logInfo, "A", block->J->data, block->n, block->n);
     }
 
@@ -386 +385 @@
             dgelss_(&n_rows, &n_x, &i, solver->jacobian_temp, &n_rows, solver->rhs, &n_rows ,solver->singular_values, &rcond, &rank, solver->rwork, &iwork, &info);
         } else {
-            memcpy(solver->jacobian_temp, solver->jacobian, n_x*n_x*sizeof(jmi_real_t));
+            memcpy(solver->jacobian_temp, block->J->data, n_x*n_x*sizeof(jmi_real_t));
             dgelss_(&n_x, &n_x, &i, solver->jacobian_temp, &n_x, solver->rhs, &n_x ,solver->singular_values, &rcond, &rank, solver->rwork, &iwork, &info);
         }
@@ -395 +394 @@
         }
         
-        if(block->callbacks->log_options.log_level >= 5){
-            jmi_log_node(block->log, logInfo, "Info", "Successfully calculated the minimum norm solution to the linear system in <block: %s>", block->label);
-        }
     }else{
         /*
@@ -449 +445 @@
     block->F(block->problem_data,block->x, solver->rhs, JMI_BLOCK_EVALUATE);
 
+    /* Check if the calculated solution from minimum norm is a valid solution to the original problem */
+    if(solver->singular_jacobian==1) {
+        double scaledMaxNorm;
+        jmi_update_f_scale(block);
+        scaledMaxNorm = calculate_scaled_residual_norm(solver->rhs, N_VGetArrayPointer(block->f_scale), block->n);
+        if(scaledMaxNorm <= block->options->res_tol) {
+            if(block->callbacks->log_options.log_level >= 5){
+                jmi_log_node(block->log, logInfo, "Info", "Successfully calculated the minimum norm solution to the linear system in <block: %s>", block->label);
+            }
+        } else {
+            info = -1;
+            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);
+            jmi_log_reals(block->log, destnode, logError, "residual", solver->rhs, block->n);
+            jmi_log_reals(block->log, destnode, logError, "scaled_max_norm", &(scaledMaxNorm), 1);
+            jmi_log_leave(block->log, destnode);
+        }
+
+    }
+
     return info==0? 0: -1;
 }
@@ -456 +471 @@
     free(solver->ipiv);
     free(solver->factorization);
-    free(solver->jacobian);
     free(solver->singular_values);
     free(solver->singular_vectors);

trunk/RuntimeLibrary/src/jmi/jmi_linear_solver.h

@@ -58 +58 @@
     int* ipiv;                     /**< \brief Work vector needed for dgesv */
     jmi_real_t* factorization;      /**< \brief Matrix for storing the Jacobian factorization */
-    jmi_real_t* jacobian;         /**< \brief Matrix for storing the Jacobian */
     jmi_real_t* jacobian_temp;         /**< \brief Matrix for storing the Jacobian */
     jmi_real_t* singular_values;  /**< \brief Vector for the singular values of the Jacobian */