SPICEDiego Circuit Simulation for Post Layout Analysis ChungKuan

SPICEDiego: Circuit Simulation for Post Layout Analysis Chung-Kuan Cheng Department of Computer Science and Engineering University of California, San Diego, CA, USA Contact: ckcheng@ucsd. edu

Outline • Background of Transient Circuit Simulation – Formulations – Problems • Matrix Exponential Circuit Simulation (Mexp) – Stiffness Problem – Skip of Regularization – Flexible Time Stepping – Errors • Experiments • Conclusions 2

Previous Works • SPICE – R. Rohrer, 1960 s, D. Pederson, 1973 – Xyce (parallel processing) – Ngspice (open source) • Matrix Exponential – C. Moler, C. Van Loan, 19 dubious ways to compute the exponential of a matrix, 1978; twenty-five years later, 2003 – A. Nauts, R. E. Wyatt, 1983; T. J. Park, J. C. Light, 1986; Y. Saad, 1992 3

VLSI Circuit: Devices + Interconnect • Multiple layers of metal with tightly coupling C and L (post layout) • Millions to billions of devices • Low margin designs with severe noises • • • C: cap. /ind. matrix G: cond. /imp. matrix x(t): v/i vector, and u(t): input sources n: millions to billions of nodes 4

Simulation Flow • Device Evaluation Circuit netlist – effective caps and current – linearization • Numerical Integration Device Evaluation – solve a system of ordinary differential equations Numerical Integration • Convergence & Error Check – stability and errors Convergence & Error Check re-evaluate stepping Step Control finish • Step Control – adjust step size 5

Transient Circuit Simulation Transient simulation: Numerical integration • Low order approx. using linear multistep methods: Backward Euler, Trapezoidal, Gear’s – Dahlquist stability barrier – Matrix decomposition for each time step h – Implicit process with Newton Raphson iteration – Local truncation error limits the time step • High order approx. using matrix exponential operators: – Matrix decomposition for flexible time step h – Explicit and stable process – Error reduces as step size increases for certain operators 6

Numerical Integration • Forward Euler • Backward Euler • Trapezoidal Method 7

Numerical Integration • Trade off between stability and performance stability high our method Backward Euler SILCA [Li & Shi, ‘ 03] ACES [Devgan & Rohrer, ‘ 97] Telescopic [Dong & Li, ‘ 10] Forward Euler low high computational effort 8

Matrix Exponential Methods • • Analytical Formulation Standard Krylov Method Rational Krylov Method Invert Krylov Method 9

Matrix Exponential Method • Analytical solution perspective – Let A=-C-1 G, b=C-1 u (C can be regularized [TCAD 2012]) • Let input be piecewise linear 10

Matrix Exponential Computation • 11

Matrix Exponential Method • Accuracy: Analytical solution – Approximate e. Ah as (I+Ah) Forward Euler – Approximate e. Ah as (I-Ah)-1 Backward Euler • Scalable: Sparse matrix-vector multiplication (Sp. MV) • Stability: Stable for passive circuits reference solution 12

Property of Standard Krylov Subspace Method • Theorem [Saad ‘ 92]: Standard Krylov method fits the first m terms in Taylor’s expansion. • Standard Krylov subspace tends to capture the eigenvalues of large magnitude • For transient analysis, the eigenvalues of small real magnitude are wanted to describe the dynamic behavior. 14

Rational Krylov Subspace • Spectral Transformation: – Shift-and-invert matrix A – Rational Krylov subspace captures slow-decay components – Use rational Krylov subspace for matrix exponential Important eigenvalue: Component that decays slowly. Not so important eigenvalue: Component that decays fast. 15

Rational Krylov Subspace: e. Ahv Rational Krylov subspace • Arnoldi orthogonalization Vm=[v 1 v 2 … vm] • Matrix reduction (Assume that v 1=v) • Matrix exponential operator 16

Invert Krylov Subspace e. Ahv Invert Krylov space • Arnoldi orthogonalization Vm=[v 1 v 2 … vm] • Matrix reduction (v 1=v) • Matrix exponential operator Comment: The calculation of A-1 vi=G-1 Cvi is much simpler when G is sparse but C is not for post layout simulation 17

Absolute Error: A Sample Test Case • 18

Matrix exponential operators vs. low order scheme • • • FE: Forward Euler Std. M: standard Mexp, m=2 Inv. M: Invert Mexp, m=2 Rat. M: Rational Mexp, m=2 BE: Backward Euler TR: Trapezoidal 19

Matrix exponential operators vs. low order scheme • 20

Matrix exponential operators vs. low order scheme Method FE 2 nd order divergent BE 2 nd order flat TR 3 rd order flat* Std Kry 2 nd order flat drop Inv Kry 1 st order flat drop Rat Kry 1 st order flat drop • 21

Rational Krylov Method vs Trapezoidal Method

Rational Krylov Space Method

Standard and Invert Krylov Methods

Standard and Invert Krylov Methods As m increases, • the curve of Krylov method shifts to the right and converges at the right end • the curve of invert Krylov method shifts to the left and converges at the left end. The desired time step h is around 10 -12 second range which is near the right end.

Three Krylov Methods •

Experiment (Rational Krylov) • Linux workstation – Intel Core i 7 -920 2. 67 GHz CPU – 12 GB memory. • Test Cases – Stiff RC mesh network (2500 Nodes) • Mexp vs. Rational Mexp – Power Distribution Network (45. 7 K~7. 4 M Nodes) • Rational Mexp vs. Trapezoidal (TR) with fixed time step (avoid LU during the simulation) 31

Experiment (I) • RC mesh network with 2500 nodes. (Time span [0, 1 ns] with a fixed step size 10 ps) stiffness definition: • Comparisons between average (mavg) and peak dimensions (mpeak) of Krylov subspace using – Standard Basis: • mavg = 115 and mpeak=264 – Rational Basis: • mavg = 3. 11, and mpeak=10 • Rational Basis-MEXP achieves 224 X speedup for the whole simulation (vs. Standard Basis-MEXP). 32

• PDN Cases Experiment (II) – On-chip and off-chip components – Low-, middle-, and high-frequency responses – The time span of whole simulation [0, 1 ps] 33

Experiment (II) • • • Mixture of low, mid, and high frequency components. 16 X speedups over TR. Difference of MEXP and HSPICE: 7. 33× 10 -4; TR and HSPICE: 7. 47× 10 -4 34

Experiment: CPU time 35

Conclusions on Rational Krylov Approach • Rational Krylov Subspace solves the stiffness problem. – No need of regularization – Small dimensions of basis. – Flexible time steps. • Adaptive time stepping is efficient to explore the different frequency responses of power grid transient simulation (considering both on-chip and off-chip components) – 15 X speedup over trapezoidal method. 36

Experiment (Invert Krylov) • Implemented in MATLAB and C/C++ – The interactions between C/C++ and MATLAB 2013 a are through MATLAB Executable (MEX) external interface with GCC 4. 7. 3. – Devices are evaluated using BSIM 3 device model and matrix stamping using C code. • In a Linux workstation – Intel CPU i 7 3. 4 GHZ – 32 GB memory 37

Numerical Results [1] Results of IEEE EMC&SI [1]. Simulation of PDN with Nonlinear Components [1] H. Zhuang, I. Kang, X. Wang, J. -H. Lin C. K. Cheng. “Dynamic Analysis of Power Delivery Network with Nonlinear Components Using Matrix Exponential Method”, IEEE EMC & SI 38 2015, (to appear)

Numerical Results (II) Results from our DAC 2015 [2] Simulation of large general circuits Due to the C matrix, the factorization complexity of (C+h. G) increased [2] H. Zhuang, W. Yu, I. Kang, X. Wang, C. K. Cheng, “An Algorithmic Framework of Large-Scale Circuit Simulation Using Exponential Integrators, ” ACM/IEEE DAC 2015 39

Conclusions • Rational/Invert Krylov methods solve the stiffness problem. – No need of regularization – Small dimensions of basis. – Explicit integration (no need of NR iteration) – Flexible time steps • Invert Krylov method – Post layout analysis 40

Conclusions: Future Works • C++ code of SPICEDiego • Mutliscale multiphysical analysis – SPICE , full wave, thermal cycling, stress simulation • Theory of matrix exponential • Distributed computation 41

Thanks and Q&A 42

Rational Krylov Approach: gamma = 1 e-14

Rational Krylov Approach: gamma= 5 e-15

Rational Krylov Approach: gamma=5 e-19