优化雅可比矩阵计算；端口转发情况下仿真结果传输方式优化

native/runtime/common.c

@@ -54,6 +54,19 @@ int native_rhs(NativeRun *r, double t, const double *y, double *dy) {
     return model_eval(t,y,dy,w);
 }
 
+/* Derivative probes use deterministic property evaluation: the generated
+ * Jacobian entry point omits cross-storage cache seeds, while the ordinary RHS
+ * retains all existing physical expressions and state-property reuse. */
+int native_jacobian_rhs(NativeRun *r, double t, const double *y, double *dy) {
+#if defined(MODEL_JACOBIAN_CANONICAL_RHS) && MODEL_JACOBIAN_CANONICAL_RHS
+    double w[NOUTPUTS];
+    r->nfev++;
+    return model_eval_jacobian(t,y,dy,w);
+#else
+    return native_rhs(r,t,y,dy);
+#endif
+}
+
 int native_append(NativeRun *r, double t, const double *y) {
     r->final_time=t; memcpy(r->final_state,y,NSTATES*sizeof(double));
     if (!r->options.record_samples) return 1;

native/runtime/cvode_solver.c

@@ -1,5 +1,7 @@
-/* CVODE owns its default numerical Jacobian and dense linear solver.
- * No project Jacobian, sparsity or derivative policy is changed here.
+/* CVODE retains its Jacobian refresh policy and dense matrix/LU solver.
+ * A compiler-proven sparsity pattern groups independent finite differences;
+ * unsupported models and patterns without grouping benefit retain CVODE's
+ * default callback.
  */
 #include "runtime.h"
 #include <cvode/cvode.h>
@@ -9,6 +11,10 @@
 #include <sunlinsol/sunlinsol_dense.h>
 #include <math.h>
 #include <string.h>
+#include <sundials/sundials_math.h>
+#ifndef MODEL_JACOBIAN_COLORED
+#define MODEL_JACOBIAN_COLORED 0
+#endif
 
 typedef struct { void *solver; N_Vector scratch; } CvDense;
 static int cv_dense(void *context, double t, double *out) {
@@ -17,16 +23,176 @@ static int cv_dense(void *context, double t, double *out) {
     memcpy(out,N_VGetArrayPointer(d->scratch),NSTATES*sizeof(double));
     return 1;
 }
+typedef struct {
+    NativeRun *run;
+    void *solver;
+    N_Vector weights;
+    SUNMatrix reference;
+    int colored;
+} CvContext;
 static int cv_rhs(sunrealtype t, N_Vector y, N_Vector dy, void *context) {
-    NativeRun *r=context;
+    NativeRun *r=((CvContext *)context)->run;
     if (!native_poll(r,r->final_time)) return -1;
     return native_rhs(r,t,N_VGetArrayPointer(y),N_VGetArrayPointer(dy)) ? 0 : 1;
 }
+
+#if MODEL_JACOBIAN_COLORED
+/* Check the generated CSC bounds and coloring once, before installing it. */
+static int valid_coloring(void) {
+    if (MODEL_JACOBIAN_COLOR_COUNT<=0 || MODEL_JACOBIAN_COLOR_COUNT>=NSTATES ||
+        model_jacobian_col_ptr[0]!=0 || model_jacobian_col_ptr[NSTATES]!=MODEL_JACOBIAN_NNZ) return 0;
+    for (int j=0;j<NSTATES;j++) {
+        if (model_jacobian_column_color[j]<0 || model_jacobian_column_color[j]>=MODEL_JACOBIAN_COLOR_COUNT ||
+            model_jacobian_col_ptr[j]<0 || model_jacobian_col_ptr[j+1]<model_jacobian_col_ptr[j] ||
+            model_jacobian_col_ptr[j+1]>MODEL_JACOBIAN_NNZ) return 0;
+        int previous=-1;
+        for (int k=model_jacobian_col_ptr[j];k<model_jacobian_col_ptr[j+1];k++) {
+            int row=model_jacobian_row_index[k];
+            if (row<=previous || row>=NSTATES) return 0;
+            previous=row;
+        }
+    }
+    for (int color=0;color<MODEL_JACOBIAN_COLOR_COUNT;color++) {
+        unsigned char used[NSTATES]={0};
+        for (int j=0;j<NSTATES;j++) if (model_jacobian_column_color[j]==color) {
+            for (int k=model_jacobian_col_ptr[j];k<model_jacobian_col_ptr[j+1];k++) {
+                int row=model_jacobian_row_index[k];
+                if (used[row]) return 0;
+                used[row]=1;
+            }
+        }
+    }
+    return 1;
+}
+
+static int jac_rhs(CvContext *context, sunrealtype t, N_Vector y, N_Vector f) {
+    if (!native_poll(context->run,context->run->final_time)) return -1;
+    context->run->jacobian_rhs++;
+#if defined(MODEL_JACOBIAN_CANONICAL_RHS) && MODEL_JACOBIAN_CANONICAL_RHS
+    return native_jacobian_rhs(context->run,t,N_VGetArrayPointer(y),N_VGetArrayPointer(f)) ? 0 : 1;
+#else
+    return native_rhs(context->run,t,N_VGetArrayPointer(y),N_VGetArrayPointer(f)) ? 0 : 1;
+#endif
+}
+
+/* Same perturbation and reciprocal-multiply order as SUNDIALS 7.4 dense DQ.
+ * No constraints are set by this runtime. If that changes, carry their vector
+ * into this context and apply CVODE's sign rule before enabling coloring. */
+static int jac_increments(CvContext *context, N_Vector y, N_Vector fy, double *increments) {
+    sunrealtype step;
+    if (CVodeGetErrWeights(context->solver,context->weights)<0 ||
+        CVodeGetCurrentStep(context->solver,&step)<0) return -1;
+    double norm=N_VWrmsNorm(fy,context->weights);
+    double minimum=norm!=0 ? 1000.0*fabs(step)*SUN_UNIT_ROUNDOFF*NSTATES*norm : 1.0;
+    double *state=N_VGetArrayPointer(y), *weight=N_VGetArrayPointer(context->weights);
+    double square_root=sqrt(SUN_UNIT_ROUNDOFF);
+    for (int j=0;j<NSTATES;j++) {
+        increments[j]=fmax(square_root*fabs(state[j]),minimum/weight[j]);
+        if (!(increments[j]>0) || !isfinite(increments[j])) return 1;
+    }
+    return 0;
+}
+
+static int dense_difference(CvContext *context, sunrealtype t, N_Vector y, N_Vector fy,
+                            SUNMatrix matrix, N_Vector trial, N_Vector ftrial, const double *increments) {
+    double *state=N_VGetArrayPointer(y), *test=N_VGetArrayPointer(trial);
+    double *base=N_VGetArrayPointer(fy), *value=N_VGetArrayPointer(ftrial);
+    memcpy(test,state,NSTATES*sizeof(double));
+    for (int j=0;j<NSTATES;j++) {
+        test[j]=state[j]+increments[j];
+        int flag=jac_rhs(context,t,trial,ftrial);
+        test[j]=state[j];
+        if (flag) return flag;
+        double inverse=1.0/increments[j];
+        double *column=SUNDenseMatrix_Column(matrix,j);
+        for (int i=0;i<NSTATES;i++) column[i]=inverse*(value[i]-base[i]);
+    }
+    return 0;
+}
+
+static int colored_difference(CvContext *context, sunrealtype t, N_Vector y, N_Vector fy,
+                              SUNMatrix matrix, N_Vector trial, N_Vector ftrial, const double *increments) {
+    double *state=N_VGetArrayPointer(y), *test=N_VGetArrayPointer(trial);
+    double *base=N_VGetArrayPointer(fy), *value=N_VGetArrayPointer(ftrial);
+    if (SUNMatZero(matrix)) return -1;
+    for (int color=0;color<MODEL_JACOBIAN_COLOR_COUNT;color++) {
+        memcpy(test,state,NSTATES*sizeof(double));
+        for (int j=0;j<NSTATES;j++) if (model_jacobian_column_color[j]==color) test[j]+=increments[j];
+        int flag=jac_rhs(context,t,trial,ftrial);
+        if (flag) return flag;
+        for (int j=0;j<NSTATES;j++) if (model_jacobian_column_color[j]==color) {
+            double inverse=1.0/increments[j];
+            double *column=SUNDenseMatrix_Column(matrix,j);
+            for (int k=model_jacobian_col_ptr[j];k<model_jacobian_col_ptr[j+1];k++) {
+                int i=model_jacobian_row_index[k];
+                column[i]=inverse*(value[i]-base[i]);
+            }
+        }
+    }
+    context->run->jacobian_colored_evals++;
+    return 0;
+}
+
+static int cv_jacobian(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix matrix,
+                        void *user, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) {
+    CvContext *context=user;
+    NativeRun *r=context->run;
+    double increments[NSTATES];
+    (void)tmp3;
+    int flag=jac_increments(context,y,fy,increments);
+    if (flag) return flag;
+#if defined(MODEL_JACOBIAN_CANONICAL_RHS) && MODEL_JACOBIAN_CANONICAL_RHS
+    /* The CVODE fy belongs to the ordinary RHS cache path. Recompute the
+       baseline using the same deterministic path as every perturbed probe;
+       mixing the two baselines would amplify cache roundoff by 1/increment. */
+    flag=jac_rhs(context,t,y,tmp3);
+    if (flag) return flag;
+    fy=tmp3;
+#endif
+    if (!context->colored) return dense_difference(context,t,y,fy,matrix,tmp2,tmp1,increments);
+    flag=colored_difference(context,t,y,fy,matrix,tmp2,tmp1,increments);
+    if (flag<0) return flag; /* cancellation/timeout must not trigger retries */
+    if (flag>0) {
+        r->jacobian_fallbacks++;
+        return dense_difference(context,t,y,fy,matrix,tmp2,tmp1,increments);
+    }
+    if (context->reference) {
+        /* Diagnostic mode validates every entry, not a sampled submatrix. */
+        flag=dense_difference(context,t,y,fy,context->reference,tmp2,tmp1,increments);
+        if (flag) return flag;
+        r->jacobian_checks++;
+        int mismatch=0;
+        for (int j=0;j<NSTATES;j++) for (int i=0;i<NSTATES;i++)
+            if (SM_ELEMENT_D(matrix,i,j)!=SM_ELEMENT_D(context->reference,i,j)) mismatch=1;
+        if (mismatch) {
+            fprintf(stderr,"{\"event\":\"jacobian-verification-mismatch\",\"time\":%.17g,\"entries\":[",(double)t);
+            int written=0;
+            for (int j=0;j<NSTATES;j++) for (int i=0;i<NSTATES;i++) {
+                double actual=SM_ELEMENT_D(matrix,i,j), expected=SM_ELEMENT_D(context->reference,i,j);
+                if (actual!=expected && written<12) {
+                    fprintf(stderr,"%s{\"row\":%d,\"column\":%d,\"colored\":%.17g,\"dense\":%.17g}",written?",":"",i,j,actual,expected);
+                    written++;
+                }
+            }
+            fprintf(stderr,"],\"state\":[");
+            for (int i=0;i<NSTATES;i++) fprintf(stderr,"%s%.17g",i?",":"",N_VGetArrayPointer(y)[i]);
+            fprintf(stderr,"]}\n");
+            r->jacobian_mismatches++; r->jacobian_fallbacks++;
+            context->colored=0; r->jacobian_colored=0;
+            if (SUNMatCopy(context->reference,matrix)) return -1;
+        }
+    }
+    return 0;
+}
+#endif
+
 static void counters(NativeRun *r, void *solver) {
     long int value=0;
     CVodeGetNumErrTestFails(solver,&value); r->rejected+=(unsigned long)value;
     CVodeGetNumJacEvals(solver,&value); r->njev+=(unsigned long)value;
     CVodeGetNumLinSolvSetups(solver,&value); r->nlu+=(unsigned long)value;
+    CVodeGetNumRhsEvals(solver,&value); r->cvode_rhs+=(unsigned long)value;
+    CVodeGetNumLinRhsEvals(solver,&value); r->linear_rhs+=(unsigned long)value;
 }
 
 int native_bdf(NativeRun *r) {
@@ -35,6 +201,7 @@ int native_bdf(NativeRun *r) {
     N_Vector y=N_VNew_Serial(NSTATES,ctx), atol=N_VNew_Serial(NSTATES,ctx), scratch=N_VNew_Serial(NSTATES,ctx);
     SUNMatrix matrix=NULL; SUNLinearSolver linear=NULL; void *solver=NULL;
     int success=0;
+    CvContext context={.run=r};
     if (!y || !atol || !scratch) goto cleanup;
     memcpy(N_VGetArrayPointer(y),r->final_state,NSTATES*sizeof(double));
     memcpy(N_VGetArrayPointer(atol),model_atol,NSTATES*sizeof(double));
@@ -44,11 +211,24 @@ int native_bdf(NativeRun *r) {
     if (!linear) goto cleanup;
     solver=CVodeCreate(CV_BDF,ctx);
     if (!solver) goto cleanup;
+    context.solver=solver;
     double t=r->options.start;
-    if (CVodeInit(solver,cv_rhs,t,y)<0 || CVodeSetUserData(solver,r)<0 ||
+    if (CVodeInit(solver,cv_rhs,t,y)<0 || CVodeSetUserData(solver,&context)<0 ||
         CVodeSVtolerances(solver,r->options.rtol,atol)<0 ||
         CVodeSetLinearSolver(solver,linear,matrix)<0 ||
         CVodeSetMaxStep(solver,r->options.max_step)<0) goto cleanup;
+#if MODEL_JACOBIAN_COLORED
+    if (valid_coloring()) {
+        context.weights=N_VClone(y);
+        if (!context.weights) goto cleanup;
+        if (r->jacobian_verify) {
+            context.reference=SUNDenseMatrix(NSTATES,NSTATES,ctx);
+            if (!context.reference) goto cleanup;
+        }
+        context.colored=1; r->jacobian_colored=1;
+        if (CVodeSetJacFn(solver,cv_jacobian)<0) goto cleanup;
+    }
+#endif
     r->starts++;
     CvDense dense={solver,scratch};
     while (t<r->options.stop) {
@@ -87,6 +267,8 @@ int native_bdf(NativeRun *r) {
     success=1;
 cleanup:
     if (solver) { counters(r,solver); CVodeFree(&solver); }
+    if (context.reference) SUNMatDestroy(context.reference);
+    if (context.weights) N_VDestroy(context.weights);
     if (linear) SUNLinSolFree(linear);
     if (matrix) SUNMatDestroy(matrix);
     if (y) N_VDestroy(y);

native/runtime/main.c

@@ -59,6 +59,13 @@ static int write_result(NativeRun *r, const char *path, const char *index_path)
     fprintf(f,"{\"success\":%s,\"status\":",!r->status?"true":"false");
     json_string(f,r->status==0?"completed":r->status==1?"cancelled":"failed");
     fprintf(f,",\"message\":"); json_string(f,r->message);
+    fprintf(f,",\"jacobianMode\":\"%s\",\"jacobianRhsCalls\":%lu,"
+              "\"jacobianColoredEvals\":%lu,\"jacobianFallbacks\":%lu,"
+              "\"jacobianChecks\":%lu,\"jacobianMismatches\":%lu,"
+              "\"cvodeRhsCalls\":%lu,\"cvodeLinearRhsCalls\":%lu",
+            !r->options.bdf?"not-used":r->jacobian_colored?"colored-difference":"dense-difference",
+            r->jacobian_rhs,r->jacobian_colored_evals,r->jacobian_fallbacks,
+            r->jacobian_checks,r->jacobian_mismatches,r->cvode_rhs,r->linear_rhs);
     fprintf(f,",\"backend\":\"native-c\",\"method\":\"%s\",\"solver\":\"%s\",\"sundialsVersion\":\"%s\","
               "\"simulatedUntil\":%.17g,\"solveSeconds\":%.17g,\"solveCpuSeconds\":%.17g,"
               "\"nfev\":%lu,\"acceptedSteps\":%lu,\"rejectedSteps\":%lu,\"stateTransitions\":%lu,"
@@ -116,6 +123,7 @@ int main(int argc, char **argv) {
             vector(stdout,y,NSTATES); fputc('\n',stdout); return 0;
         }
         if (!strcmp(arg,"--solve-only")) { r.options.record_samples=0; continue; }
+        if (!strcmp(arg,"--verify-jacobian")) { r.jacobian_verify=1; continue; }
         if (i+1==argc) return 64;
         const char *value=argv[++i];
         if (!strcmp(arg,"--output")) output=value;