From 13210bbf42495ade67e66df99c2ca8be04844bfd Mon Sep 17 00:00:00 2001 From: Meisam Bahadori Date: Sun, 26 Jul 2026 18:39:57 +0200 Subject: [PATCH] pa-111 --- src/include/ngspice/cktdefs.h | 5 ++ src/include/ngspice/optdefs.h | 1 + src/include/ngspice/tskdefs.h | 1 + src/maths/ni/niiter.c | 100 +++++++++++++++++++++++++++++++ src/spicelib/analysis/cktdest.c | 2 + src/spicelib/analysis/cktdojob.c | 1 + src/spicelib/analysis/cktntask.c | 2 + src/spicelib/analysis/cktsopt.c | 5 ++ 8 files changed, 117 insertions(+) diff --git a/src/include/ngspice/cktdefs.h b/src/include/ngspice/cktdefs.h index d9d22327b..0a65c633b 100644 --- a/src/include/ngspice/cktdefs.h +++ b/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 */ diff --git a/src/include/ngspice/optdefs.h b/src/include/ngspice/optdefs.h index 25d61eebd..8890f2672 100644 --- a/src/include/ngspice/optdefs.h +++ b/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 diff --git a/src/include/ngspice/tskdefs.h b/src/include/ngspice/tskdefs.h index f59ccec9b..35c77f863 100644 --- a/src/include/ngspice/tskdefs.h +++ b/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 */ diff --git a/src/maths/ni/niiter.c b/src/maths/ni/niiter.c index 9f06a62df..ed61e9f3c 100644 --- a/src/maths/ni/niiter.c +++ b/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) diff --git a/src/spicelib/analysis/cktdest.c b/src/spicelib/analysis/cktdest.c index 8a868a5f4..80390aafc 100644 --- a/src/spicelib/analysis/cktdest.c +++ b/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); diff --git a/src/spicelib/analysis/cktdojob.c b/src/spicelib/analysis/cktdojob.c index 3a2922e84..5f9a9c54d 100644 --- a/src/spicelib/analysis/cktdojob.c +++ b/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; diff --git a/src/spicelib/analysis/cktntask.c b/src/spicelib/analysis/cktntask.c index e143798cf..083be6672 100644 --- a/src/spicelib/analysis/cktntask.c +++ b/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; diff --git a/src/spicelib/analysis/cktsopt.c b/src/spicelib/analysis/cktsopt.c index 816f603bc..ad7c1023b 100644 --- a/src/spicelib/analysis/cktsopt.c +++ b/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,