pa-111

src/include/ngspice/cktdefs.h

@@ -270,6 +270,11 @@ struct CKTcircuit {
     unsigned int CKTkeepOpInfo:1; /* flag for small signal analyses */
     unsigned int CKTcopyNodesets:1; /* NodesetFIX */
     unsigned int CKTnodeDamping:1; /* flag for node damping fix */
+    unsigned int CKTlinesearch:1;  /* Enhancement-111: adaptive damped-Newton line search */
+    double CKTlsMerit;             /* E-111: this iteration's residual merit ||F|| = ||G*x-b|| */
+    double *CKTlsXk;               /* E-111: line-search scratch (saved x_k) */
+    double *CKTlsD;                /* E-111: line-search scratch (Newton step d) */
+    int CKTlsBufSz;                /* E-111: allocated size of the LS scratch buffers */
     unsigned int CKTcheckpoint:1;  /* Enhancement-131: DCtran should continue from a restored
                                       checkpoint (keep loaded state, build a fresh output plot) */
     double CKTabsDv;            /* abs limit for iter-iter voltage change */

src/include/ngspice/optdefs.h

@@ -136,6 +136,7 @@ enum {
     OPT_LTEABSTOL,
     OPT_LTETRTOL,
     OPT_NEWTRUNC,
+    OPT_LINESEARCH,     /* Enhancement-111: adaptive damped-Newton line search */
 };
 
 #ifdef XSPICE

src/include/ngspice/tskdefs.h

@@ -69,6 +69,7 @@ struct TSKtask {
     unsigned int TSKkeepOpInfo:1; /* flag for small signal analyses */
     unsigned int TSKcopyNodesets:1; /* flag for nodeset copy */
     unsigned int TSKnodeDamping:1;  /* flag for node damping */
+    unsigned int TSKlinesearch:1;   /* Enhancement-111: adaptive damped-Newton line search */
     unsigned int TSKnoopac:1; /* flag for no OP calculation before AC */
     double TSKabsDv;                 /* abs limit for iter-iter voltage change */
     double TSKrelDv;                 /* rel limit for iter-iter voltage change */

src/maths/ni/niiter.c

@@ -110,6 +110,29 @@ NIiter(CKTcircuit *ckt, int maxIter)
                 ckt->CKTniState |= NISHOULDREORDER;
             }
 
+            /* Enhancement-111: residual merit ||F(x_k)|| = ||G*x_k - b|| for the
+             * globalized Newton line search. Computed here -- after the matrix
+             * is loaded and preordered but BEFORE it is LU-factored (spMultiply
+             * requires an unfactored matrix). G is the just-loaded Jacobian,
+             * x_k = CKTrhsOld, b = CKTrhs. F is the KCL residual (a current) --
+             * the merit function ngspice's iterate-based Newton otherwise lacks.
+             * CKTrhsSpare is scratch (free until the solve below). */
+            if (ckt->CKTlinesearch && ckt->CKTrhsSpare && (iterno > 1)) {
+                int sz = SMPmatSize(ckt->CKTmatrix);
+                int k;
+                double m = 0.0;
+                SMPmultiply(ckt->CKTmatrix, ckt->CKTrhsSpare, ckt->CKTrhsOld,
+                            NULL, NULL);
+                for (k = 1; k <= sz; k++) {
+                    double resid = ckt->CKTrhsSpare[k] - ckt->CKTrhs[k];
+                    double w = fabs(resid) /
+                        (ckt->CKTabstol + ckt->CKTreltol * fabs(ckt->CKTrhsSpare[k]));
+                    if (w > m)
+                        m = w;
+                }
+                ckt->CKTlsMerit = m;
+            }
+
             if (ckt->CKTniState & NISHOULDREORDER) {
                 startTime = SPfrontEnd->IFseconds();
 
@@ -323,6 +346,83 @@ NIiter(CKTcircuit *ckt, int maxIter)
             }
         }
 
+        /* Enhancement-111: globalized (damped) Newton via Armijo backtracking
+         * line search (option `linesearch`, OFF by default). Runs only on the
+         * non-convergence path. Using the residual merit ||F|| = ||G*x - b||
+         * (the KCL current mismatch, computed above before factorization), it
+         * damps the full Newton step x_k -> x_full by the largest lambda in
+         * {1, 1/2, 1/4, ...} that gives a sufficient decrease of ||F|| (Armijo).
+         * Each trial RE-LOADS the devices at the trial point x_k + lambda*d and
+         * re-evaluates ||F|| on the (SMPclear-reset, unfactored) matrix. At/near
+         * a solution the full step already reduces ||F||, so lambda = 1 is
+         * accepted on the first trial (result-neutral); backtracking only kicks
+         * in on genuine overshoot. This gives ngspice the principled globalized
+         * Newton it lacks -- the merit is the real residual, not the iterate
+         * change. */
+        if (ckt->CKTlinesearch && (ckt->CKTnoncon != 0) &&
+            ((ckt->CKTmode & MODETRANOP) || (ckt->CKTmode & MODEDCOP)) &&
+            (ckt->CKTmode & MODEINITFLOAT) && (iterno > 1))
+        {
+            int sz = SMPmatSize(ckt->CKTmatrix);
+            int k, saved_noncon = ckt->CKTnoncon;
+            double lambda = 1.0, merit_k = ckt->CKTlsMerit;
+
+            if (ckt->CKTlsBufSz < sz + 1) {
+                FREE(ckt->CKTlsXk);
+                FREE(ckt->CKTlsD);
+                ckt->CKTlsXk = TMALLOC(double, sz + 1);
+                ckt->CKTlsD  = TMALLOC(double, sz + 1);
+                ckt->CKTlsBufSz = sz + 1;
+            }
+            /* save x_k and the full Newton step d = x_full - x_k */
+            for (k = 1; k <= sz; k++) {
+                ckt->CKTlsXk[k] = ckt->CKTrhsOld[k];
+                ckt->CKTlsD[k]  = ckt->CKTrhs[k] - ckt->CKTrhsOld[k];
+            }
+            for (;;) {
+                double trial_merit = 0.0;
+                for (k = 1; k <= sz; k++)
+                    ckt->CKTrhsOld[k] = ckt->CKTlsXk[k] + lambda * ckt->CKTlsD[k];
+                /* Reset the device state (junction-voltage limiting reference)
+                 * to x_k before every trial, so each trial load limits relative
+                 * to the SAME point -- SPICE limiting is stateful and would
+                 * otherwise drift across trials and corrupt the iteration. */
+                if (OldCKTstate0 && ckt->CKTstate0)
+                    memcpy(ckt->CKTstate0, OldCKTstate0,
+                           (size_t) ckt->CKTnumStates * sizeof(double));
+                if (CKTload(ckt))       /* trial load failed -> stop backtracking */
+                    break;
+                SMPmultiply(ckt->CKTmatrix, ckt->CKTrhsSpare, ckt->CKTrhsOld,
+                            NULL, NULL);
+                for (k = 1; k <= sz; k++) {
+                    double resid = ckt->CKTrhsSpare[k] - ckt->CKTrhs[k];
+                    double w = fabs(resid) /
+                        (ckt->CKTabstol + ckt->CKTreltol * fabs(ckt->CKTrhsSpare[k]));
+                    if (w > trial_merit)
+                        trial_merit = w;
+                }
+                /* Armijo sufficient-decrease (c = 1e-4); floor lambda at 1/64 */
+                if (trial_merit <= (1.0 - 1.0e-4 * lambda) * merit_k ||
+                    lambda <= 1.0 / 64.0)
+                    break;
+                lambda *= 0.5;
+            }
+            /* Accept the (damped) step. Put x_trial into CKTrhs and restore
+             * CKTrhsOld = x_k so the SWAP below advances to x_trial. Roll the
+             * device state back to x_k so the trial loads leave NO state trace:
+             * the next iteration then loads at x_trial with the x_k limiting
+             * reference -- exactly the cadence a normal (un-line-searched)
+             * iteration would have. Only the chosen x_trial position persists. */
+            for (k = 1; k <= sz; k++) {
+                ckt->CKTrhs[k]    = ckt->CKTrhsOld[k];
+                ckt->CKTrhsOld[k] = ckt->CKTlsXk[k];
+            }
+            if (OldCKTstate0 && ckt->CKTstate0)
+                memcpy(ckt->CKTstate0, OldCKTstate0,
+                       (size_t) ckt->CKTnumStates * sizeof(double));
+            ckt->CKTnoncon = saved_noncon; /* trial loads dirtied it */
+        }
+
         if (ckt->CKTmode & MODEINITFLOAT) {
             if ((ckt->CKTmode & MODEDC) && ckt->CKThadNodeset) {
                 if (ipass)

src/spicelib/analysis/cktdest.c

@@ -86,6 +86,8 @@ CKTdestroy(CKTcircuit *ckt)
     FREE(ckt->CKTrhs);
     FREE(ckt->CKTrhsOld);
     FREE(ckt->CKTrhsSpare);
+    FREE(ckt->CKTlsXk);        /* Enhancement-111 */
+    FREE(ckt->CKTlsD);         /* Enhancement-111 */
     FREE(ckt->CKTirhs);
     FREE(ckt->CKTirhsOld);
     FREE(ckt->CKTirhsSpare);

src/spicelib/analysis/cktdojob.c

@@ -102,6 +102,7 @@ CKTdoJob(CKTcircuit* ckt, int reset, TSKtask* task)
     ckt->CKTkeepOpInfo = task->TSKkeepOpInfo;
     ckt->CKTcopyNodesets = task->TSKcopyNodesets;
     ckt->CKTnodeDamping = task->TSKnodeDamping;
+    ckt->CKTlinesearch = task->TSKlinesearch; /* Enhancement-111 */
     ckt->CKTabsDv = task->TSKabsDv;
     ckt->CKTrelDv = task->TSKrelDv;
     ckt->CKTtroubleNode = 0;

src/spicelib/analysis/cktntask.c

@@ -72,6 +72,7 @@ CKTnewTask(CKTcircuit *ckt, TSKtask **taskPtr, IFuid taskName, TSKtask **defPtr)
         tsk->TSKkeepOpInfo      = def->TSKkeepOpInfo;
         tsk->TSKcopyNodesets    = def->TSKcopyNodesets;
         tsk->TSKnodeDamping     = def->TSKnodeDamping;
+        tsk->TSKlinesearch      = def->TSKlinesearch; /* Enhancement-111 */
         tsk->TSKabsDv           = def->TSKabsDv;
         tsk->TSKrelDv           = def->TSKrelDv;
         tsk->TSKnoopac          = def->TSKnoopac;
@@ -135,6 +136,7 @@ CKTnewTask(CKTcircuit *ckt, TSKtask **taskPtr, IFuid taskName, TSKtask **defPtr)
         tsk->TSKkeepOpInfo      = 0;
         tsk->TSKcopyNodesets    = 0;
         tsk->TSKnodeDamping     = 0;
+        tsk->TSKlinesearch      = 0; /* Enhancement-111: off by default */
         tsk->TSKabsDv           = 0.5;
         tsk->TSKrelDv           = 2.0;
         tsk->TSKepsmin          = 1e-28;

src/spicelib/analysis/cktsopt.c

@@ -160,6 +160,9 @@ CKTsetOpt(CKTcircuit *ckt, JOB *anal, int opt, IFvalue *val)
     case OPT_NODEDAMPING:
         task->TSKnodeDamping = (val->iValue != 0);
         break;
+    case OPT_LINESEARCH: /* Enhancement-111 */
+        task->TSKlinesearch = (val->iValue != 0);
+        break;
     case OPT_ABSDV:
         task->TSKabsDv = val->rValue;
         break;
@@ -359,6 +362,8 @@ static IFparm OPTtbl[] = {
         "Copy nodesets from device terminals to internal nodes" },
  { "nodedamping", OPT_NODEDAMPING, IF_SET|IF_FLAG,
         "Limit iteration to iteration node voltage change" },
+ { "linesearch", OPT_LINESEARCH, IF_SET|IF_FLAG,
+        "Adaptive damped-Newton line search on the weighted step norm" },
  { "absdv", OPT_ABSDV, IF_SET|IF_REAL,
         "Maximum absolute iter-iter node voltage change" },
  { "reldv", OPT_RELDV, IF_SET|IF_REAL,