Changeset 9082

Timestamp:

Jul 18, 2016 10:21:49 AM (3 years ago)

Author:

Christian Andersson

Message:

Cleaned up computation of the scaled vector norm. Note though that this should be moved to jmi_maths in the future. Related to ticket:4815

Location:

trunk/RuntimeLibrary/src/jmi

Files:

• jmi_block_solver.c (modified) (1 diff)
• jmi_block_solver.h (modified) (2 diffs)
• jmi_linear_solver.c (modified) (2 diffs)

trunk/RuntimeLibrary/src/jmi/jmi_block_solver.c

@@ -977 +977 @@
 }
 
-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;
-}
-
+int jmi_scaled_vector_norm(jmi_real_t *x, jmi_real_t *scale, jmi_int_t n, jmi_int_t NORM, jmi_real_t* out) {
+    int ef = 0;
+    
+    if (NORM == JMI_NORM_MAX) {
+        int i;
+        jmi_real_t tmp;
+        *out = 0.0;
+        for (i = 0; i < n; i++) {
+            tmp = JMI_ABS(x[i]*scale[i]);
+            *out = *out >= tmp ? *out : tmp;
+        }
+    } else {
+        ef = -1;
+    }
+    
+    return ef;
+}
+

trunk/RuntimeLibrary/src/jmi/jmi_block_solver.h

@@ -65 +65 @@
 #define JMI_VAR_NOT_USED(x) ((void)x)
 
+#define JMI_NORM_MAX 0
+
 /** \brief Jacobian variability for the linear solver */
 typedef enum jmi_block_solver_jac_variability_t {
@@ -395 +397 @@
 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 Computes the scaled vector norm */
+int jmi_scaled_vector_norm(jmi_real_t *x, jmi_real_t *scale, jmi_int_t n, jmi_int_t NORM, jmi_real_t* out);
 
 /** \brief Check and log illegal iv inputs */

trunk/RuntimeLibrary/src/jmi/jmi_linear_solver.c

@@ -483 +483 @@
     /* Check if the calculated solution from minimum norm is a valid solution to the original problem */
     if(solver->singular_jacobian==1) {
-        double scaledMaxNorm;
+        jmi_real_t scaled_max_norm;
+        int ef = 0;
         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) {
+        
+        ef = jmi_scaled_vector_norm(solver->rhs, N_VGetArrayPointer(block->f_scale), block->n, JMI_NORM_MAX, &scaled_max_norm);
+        if (ef == -1) {
+            jmi_log_node(block->log, logError, "NormFailure", "Failed to compute the scaled residual norm to the linear system in <block: %s>", block->label);
+        }
+        
+        if(scaled_max_norm <= 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);
@@ -494 +500 @@
             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, "residuals", solver->rhs, block->n);
-            jmi_log_reals(block->log, destnode, logError, "scaled_max_norm", &(scaledMaxNorm), 1);
-            jmi_log_reals(block->log, destnode, logError, "tolerance", &(block->options->res_tol), 1);
+            jmi_log_reals(block->log, destnode, logError, "scaled_max_norm", &(scaled_max_norm), 1);
+            jmi_log_reals(block->log, destnode, logError, "tolerance", &(block->options->res_tol), 1);
             jmi_log_leave(block->log, destnode);
         }