21
PLMAP - A Piecewise Linear MOS Circuit Analysis Program Qingjian Yu * and 0. Wing Department
of Electrical
Received 1 November
Engineering,
Columbia
University,
New
York,
NY
10027,
U.S.A.
1983
Abstract. In this paper we describe a MOS circuit simulation program PLMAP (Piecewise Linear MOS Circuit Analysis Program) which performs efficient time domain simulation of LSI MOS circuits. PLMAP uses a piecewise linear model for MOS transistors. It uses a relaxation method to analyze the piecewise linear circuit at one time point after another. The relaxation method consists of a 2-level decomposition which partitions the whole circuit into 2-level subcircuits and a scheduling method which determines the sequence of analysis of subcircuits. At each time point, the subcircuits are analyzed sequentially and iteratively by using Gauss-Seidel iteration together with a prediction formula to determine the initial guess of solution. Experiences have shown that PLMAP is more efficient in MOS circuit analysis than SPICE.
Keywords. scheduling
Piecewise-linear of subcircuits.
MOS
model,
2-level
circuit
decomposition,
relaxation
method,
1. Introduction
Circuit simulation techniques can be’ divided into two categories. One is exemplified by SPICB [l] and ASTAP [2], in which each circuit element is modeled in its finest detail and the circuit is representedby a set of nonlinear algebraic-differential equations. The equations are solved numerically, aided by sparsematrix techniques. The other is representedby MOTIS [3] and its derivatives [4-61, in which the circuit elements are crudely modeled and the circuit equations are solved only to first-order accuracy. While the programs of the first category give accurate and usually reliable results, they are too slow for use in the simulation of * On leave from East China Engineering North-Holland INTEGRATION,
the VLSI journal
0167-9260/84/83.00
0 1984, Elsevier
Institute,
Nanjing,
China.
2 (1984) 27-48 Science Publishers B.V. (North-Holland)
28
Qingjian Yu, 0. Wing / PLMAP - A MOS analysisprogram
large-scale integrated circuits. Programs in the second category are fast, but accuracy is sacrificed and sometimes large error or numerical instability may occur. Recently, there have been new attempts to develop circuit analysis programs which are fast and accurate. An example is the program based on the waveform relaxation method. This method is efficient but it requires a large amount of memory to store the waveforms, and it suffers from convergenceproblem when strong feedback exists [7,8]. In this paper we describea circuit simulation program PLMAP suitable for LSI MOS circuits which has both fast speedand good accuracy.Fast speedis achieved by using a piecewise linear model for MOS transistors and by using a 2-level relaxation method to solve the circuit equations. Accuracy is obtained by controlling the absolute error in the piecewise linear approximation and the maximum relative error in the relaxation process. We shall first describe how a piecewise linear model of a transistor is constructed.The model is derived from a quasi-static representation in which charge is conserved.The relaxation method is presented in Section 3, which describeshow a circuit is partitioned into subcircuits, each of which may be further partitioned. The scheduling algorithm used to decide the sequence of analysis of the subcircuits is described in Section 4. The overall organization of PLMAP is given in Section 5, and test examples are given in Section 6.
m P
Yu was born in Shanghai, China on March 27, 1939. He graduated at East China Engineering Institute in Nanjing, China in 1961. Since then, he has been working there with the Department of Electrical Engineering. At present, he is an Associate Professor of Electrical Engineering. During the period July 1982-June 1984 he was a visiting scholar at Columbia University in New York. His professional interests include circuit theory and computer-aided design of VLSI circuits. Dr. Qingjian Yu is a member of the Electronic Society of China and also of the Automatic Control Society of China. Qingjian
Wing received the B.S. degree from the University of Tennessee, Knoxville. in 1950. the M.S. degree from the Massachusetts Institute of Technology, Cambridge, in 195% and the Engrg.Sci.D. degree from Columbia Uruversity, New York, NY, in 1959. From 1952 to 1956 he was a member of the technical staff at Bell Telephone Laboratories. Since 1959 he has been a member of the faculty of the Department of Electrical Engineering, Columbia University, where he is now Chairman and Professor of Electrical Engineering. He also served as Chairman from 1974 to 1978. He was a Fulbright Visiting Lecturer at the Institute of Electronics, Chiao-Tung University, Taiwan, in 1961, a Ford Foundation Engineering Resident at IBM Thomas J. Watson Research Center in 1966-1967, a Visiting Professor at the Institute of Circuit Theory, Technical University of Denmark in 1973, a Visiting Professor at the Indian Institute of Technology at Kanpur, India in 1978, a Visiting Professor at the South China Institute of Technology, Canton, China in 1979, Fulbright Senior Lecturer at the Technical University of Eindhoven, The Netherlands in 1979-1980, and a Visiting Professor at Jiao-Tong University, Shanghai., China in 1982. Dr. Wmg served as Editor of the IEEE Transactions on Circuits and Systemsin 1974-1976, and President of the IEEE Circuits and Systems Society in 1978. He is a Fellow of the IEEE. Omar
Qingjian
Yu, 0. Wing / PLMAP
2. Piecewise linear model for MOS
- A MOS analysis program
29
transistors
The use of piecewise linear model for MOS transistors has two main advantages: (1) The model is simpler than the general nonlinear model but can still be sufficiently accurate. (2) During transient analysis, from one time point to the next, the changeof the node voltagesin a subcircuit is relatively small, so that in most casesthe operating regions of the transistors will remain unchanged,and the time spent on updating the matrix of the subcircuit equations and implementing the LU decomposition for the matrix can be saved.When some of the characteristics of the transistors do change their regions, the update of the circuit equations is much simpler than in the casewhen a general nonlinear model is used. In order to get a relatively precise MOS model, both the accuracy of the drain current model and that of the transistor charge model are important. Furthermore, the charge model must have the property of chargeconservation; otherwise, the model will causenonconvergenceproblems [9]. In PLMAP, a charge conserving model is used. This model is derived from [9] with some modifications, and it consists of four parts: (i) the drain current equations, (ii) the intrinsic charge equations, (iii) the p-n junction charge equations, and (iv) the p-n junction current equations. The drain current equations are as follows: In the accumulation (F/g,,< V,) and subthreshold region (V,,, + V, < l$,s< V,) we have &=O.
In the linear region ( Vslgs > V,, Valgd > V,) we have
Id=k[(vg,- V;)2-(l/gd-u’]In the saturation region ( Vss> V,, V,sd< V,) we have
&=k[(T/,,- v)2+@?3s-v,s,,at)(v,sv’)] =k[(l-d(V,,- v,)‘+ar/ds(vgs01.
(2)
In the above expressions,d, s, g and b are the drain, source, gate and substrate of a MOS transistor, V, is the threshold voltage and
where V, is the flat band voltage, +r the Fermi level and y the body factor; and k = OS&,
- w/l.
The term V,,,,, is the saturation voltage between the drain and source; and 4 v;gs- V,) is the slope of the 1, - V,, characteristic in the saturation region.
30
Qingiian
Yu, 0. Wing / PLMAP
A MOS analysis
-
program
With reference to (1) and (2), there are 3 types of nonlinear where x = I$ -
f (4 = x2, g(x)
functions:
V,,-y;
where x = Vsb+ 2GF,
= 6
and +,y)
where x=
=‘xy,
V,, and y=
VP-
V,.
The piecewise linear approximation of the one-variable functions, e.g., f(x), can be made by controlling the maximum error of the function value in the piecewise linear regions from one to another. For example, suppose in the k th region [xk, xk+i] that f( x) is approximated by fdx)
= %(x
-x/A
with maximum permissible = 0, we use the condition
(4
error 8. Starting from k = 1, x, = 0 and f,(x,)
-fWl = 6
m4fk(xj
to determine
+f/d(x/J
= f(xl)
(5)
uk, and
f (xk+l > -fdxk+l)
(6)
= 6
to determine xk+i until xk+, 2 x,.,.,~,., where x,,, is the maximum value of x within the operating range of the transistor. The 2-variable function T(X, u) = xy can be approximated by using the first order Taylor series expansion in each region D: T~(x,Y)=Yjx+xjY-xiYj, D:
Ix-y-(xi-yj)Ia2,
xi=
(i-
(7) [x+Y-(xi+Yj)I~a~
(8)
where l)a,
yj=
a=2& and 6 is the maximum the approximation are shown in Appendix A. The intrinsic charge following notations are
(j-
l)a,
i,j=1,2
equations used:
of a MOS transistor
Qgd= -Qa Qgb= -Qb,
and co = Co,wl.
Qgs= Q,,, =
region we have 0,
(9)
00) permissible error for the approximation. The regions for shown in Fig. 1. The derivation of this approximation is
Q,=Q,s+Q,+Q,,, Qgs= -Q,,
In the accumulation
,-**,
Qgb=
Cdv,,
- G,>.
are as follows,
where the
Qingjian
Yu, 0. Wing / PLMAP
-A
MOS analysis program
31
Y Sa
3a
a b3bL 0
a
3a
X
Sa
Fig. 1.
In the subthreshold region we have Qgs = Qgd =
Qgb= 0.5Coy2[ - 1 + ,/l -I-4( Vgb- v&u”].
0,
(11)
In the saturation region we have Q,d = 0,
Q,,=+c,(v,s-
v,),
Qg,, = vCo(K,
+ 2dd’2-
04
In the linear region we have Q,b
= vCo( K, + %)1’2,
(vps - mv, - vo v,)-(v,- v,>+(yy v,)+(v,,-K>I'
2&s-
v,)-
("-
(v,-
y)(l/gd-
v)+&d-
')
vt>
03)
1 *
The piecewise linear approximation of the nonlinear functions in the charge model is made by controlling the maximum error of both the function values and their derivatives. For example, suppose in region [xk, xk+i] that the nonlinear function h(x) = (1 + x)li2 in (ll), where x = 4( I& - V,)/y2, is approximated by hkb) = Ukb - Xk) +hk-lbk). Starting from k = 1, xi = 0 and h,(x,) = h(x,) = 1, we use the equations hk(x)=h(x) and u,=h(x) to determine ak, and use the equations h,(x:+,)
-+:+,)
= a,, ak - +;+I)
= 6,
(14
and
05) to determine xk+i, where 6, and 8, are the maximum permissible errors for the function value and its derivative. The 2-variable nonlinear function xy/(x + y) in (13), where x = Vss- V, and y = Vsd- V,, can be approximated in the following way. xk+l
=
Inin{
x:+1,
xi+1
1
32
Qingiian
Yu, 0. Wing / PLMAP
- A MOS analysis program
Let xy/(x+y)=x*z/(l+z) where z=y/x, O
of circuit equations
and decomposition
of circuit
The use of relaxation method to solve the set of piecewise linear circuit equations is another reasonwhy PLMAP is fast. In this program, the Gauss-Seidel iteration is used at each time point. At each iteration, the initial guessis obtained from a first- or second-orderpredictor. At some time point tk+l, let the given circuit be replaced by its companion network based on either the Backward-Euler or the trapezoidal formula, and let the set of nodal equations of the companion network be as follows:
j&,x2,
. ..9&J=o>
wheref,andxi(i=1,2, . . . . n) are vectors. The components of vector xi are node voltages, and fi( x) = 0 is the set of node equations corresponding to xi. Then the (m + 1)st Gauss-Seidel iteration can be expressedas follows:
fl(Xl,-q, . ...x.“)=o, f2(x;1+1,x2, . ... x:)=0,
(17)
fn( x;“+1, xy+l, . . ., x.) = 0, where the first guessx0 is determined by the prediction formula x”=x(tk)+h*i(tk)+P*0.5*h2*2(tk),
where /3 = 6 or 1 depending on whether the Backward-Euler or the trapezoidal companion network is used.
Qingiian
Yu, 0. Wing / PLMAP
-
A MOS
analysicprogram
33
This relaxation method can be regarded as a special case of the waveform relaxation method. It is known that under mild conditions [7], which are satisfied by most practical MOS circuits, the waveform relaxation method will converge;so will this method. Note that at time tk+,, the first guessof this method is based on the true solution at the previous time and it will be near the solution x( t,,,) of (16). Therefore, this method needs fewer iterations than the waveform relaxation method. The relaxation from (16) to (17) is equivalent to the decomposition of the whole circuit into subcircuits [7], where h(x) corresponds to subcircuit i. The nodes with unknown voltages in h(x) = 0 are the internal nodes of the subcircuit, and those with known node voltages are its input nodes. It is desired that the size (the number of the internal nodes) of each subcircuit be as small as possible to reduce the time for solving its node equations. On the other hand, it is also desired that the basic logic cells in a circuit, e.g., NAND gates,NOR gates, etc., not be decomposed,becausetheselogic cells are individually strongly coupled. In this program, 24evel decomposition is used. A first level subcircuit is called a compound subcircuit which may have one or more internal nodes and some small subcircuits in it. A second level subcircuit has only one internal node and is part of a first-level subcircuit. There are two cases on how a MOS transistor is included in a subcircuit. In one case, its drain current passes through some internal node of the subscircuit, e.g., transistor Ml in subcircuit 1 in Fig. 2. Such a transistor is called a basic transistor of the subcircuit. The other caseis that only the gate of the transistor is an internal node of the subcircuit, e.g., transistor M2 in subcircuit 1 in Fig. 2. In this case only its capacitancesconnected to the gate are included in the subcircuit. Such a transistor is called a coupling transistor of the subcircuit. Accordingly, the input nodes of a subcircuit are divided into two groups. One group consists of the drains, sourcesand substrates of the coupling transistors, provided that they are not the terminals of any ideal voltage sources. These are called feedback input nodes. All the other input nodes are called basic input nodes. The first-level decomposition is implemented by forming the subcircuits one after another. Its process is as follows. Choose any node i whose voltage is unknown and which has not been included in any subcircuit as an internal node. Declare a new subcircuit with i as its first internal node. Search all the elements with i as one of their terminals (for a MOS transistor, it is either the drain or the source) and include them as basic elements in the subcircuit. For each element included, check its other terminalj to seeif it has been designated as an internal node of the subcircuit. If it has, then abandon it; otherwise, do further check. Ifj is a terminal of an ideal voltage sourceor an internal node of another subcircuit, then its voltage will be known when the subcircuit is being analyzed. In this case, include it as a basic input node of the subcircuit. If the element is a 2-terminal element or a pass transistor, then includej as a basic input node of the subcircuit, too. (In PLMAP, all pass transistors are identified as such by the user as part of the input data.) This means that j will be included in another subcircuit as an internal node and this subcircuit will stop being expanded through the element. In
34
Qingjian
Yu, 0.
Wing
/ PLMAP
-
A MOS
program
M4
n3
2
1
“c
analysis
n2
Ml
==
c
gs
0
I4
a
4 n3
C
Fig. 2. (a) For simplicity, only two of the capacitances of M2, Css and Csd, are explicitly drawn. (b) Subcircuit 1: internal node - 1; basic input nodes - 3,4; feedback input node - 2; basic transistors - Ml, M3; couplingfransistor - M2. (c) Subcircuit 2: internal node - 2; basic input nodes -1, 4; basic transistors - M2, M4.
other cases,include j in the subcircuit as an internal node and search the basic elements connected to it. This process continues until no new internal node can be found for the subcircuit. Then, check each gate node of the transistors in the subcircuit. If it is not an internal node of the subcircuit, then declare it as a basic input node of the subcircuit. Finally, for each internal node of the subcircuit, search the transistors with their gate nodes connected to it but not included in the subcircuit as basic transistors, and include them in the subcircuit as its coupling transistors together with their drains, sourcesand substratesas the feedback input nodes of the subcircuit if they are not the terminals of any ideal voltage sources. Repeat the above process for all the unknown voltage nodes. The whole circuit will be decomposed into several first level subcircuits, and each node with unknown voltage will be included in one and only one subcircuit as an internal node. In the above process, a subcircuit cannot be expanded through a 2-terminal element or a pass transistor. This is based on the consideration that usually the
Qingiian
Yu, 0. Winn / PLMAP
- A MO.!? analysis program
35
C2
0
Fig. 3.
two parts of a circuit which are separated by such elements are not in the same basic logic .cell, and sometimes a pass transistor may be connected to too many transistors as in the case of a memory circuit. If a subcircuit is expanded through a pass transistor, then the subcircuit may be too large. An example of the implementation of the algorithm is shown in Fig. 3 and Table 1. It can be seen that the combinational logic cell will be included in one subcircuit. If a compound subcircuit has too many internal nodes, then a second level decomposition will be done for it. In this case, each internal node in the compound subcircuit will correspond to a second level subcircuit. The elements and the input nodes of the second level subcircuits can be determined by examining the elements in the compound subcircuit one after another. An example of the second level decomposition is shown in Fig. 4 and Table 2, where
Table 1 Subcircuit
Internal nodes
Basic elements
1
Ml, M2, M4, Cl
2
M3
3
M5, M6, C2
1
2
Basic input nodes
Coupling transistors
Feedback input nodes
4,5,6,7
M5
3
4 1
36
Qingjian
Yu. 0. Wing / PLMAP
- A MOS analysis program
Ml 1
Ml5 9 Ml6
Fig. 4.
some elements appear in 2 subcircuits as basic elements. After the whole circuit is decomposed into first and second level subcircuits, the Gauss-Seidel iteration will be applied in two nested loops. The outer loop is for the compound subcircuits and the inner loop for the second level subcircuits in each compound subcircuit. When the circuit consists of some comphcated logic cells, this 2-level decomposition and iteration algorithm will be more efficient than when only one level decomposition is used.
Table 2 Subcircuit
Internal node
Basic transistors
Input nodes
1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10
Ml, M2, M5, M7, Ml1 M2, M3 M3, M4 MS, M6 M7, M8, M9 M9, Ml0 M12, M13, Ml5 M13, Ml4 M15, Ml6 Mll, Ml2
11, 12,15, 18, 23, 2, 4, 5, 10 12,13, 1, 3 13, 14, 2 15, 16, 1 17, 18, 19, 1, 6 19, 20, 5 24, 21,25, 8, 9, 10 21,22, 7 25,26, 7 23, 24, 1, ‘!
Qingjian
4. Scheduling
Yu, 0. Wing / PLMAP
- A MOS analysis program
37
of subcircuits
For a given circuit, the speed of convergenceof the Gauss-Seidel iteration mainly depends on the sequenceof analysis of the subcircuits. It is desired that the sequencecoincides with the way in which the signal flows. In this program, scheduling is applied to the first level subcircuits only, since there are no explicit signal flow paths among the secondlevel subcircuits and the sequenceof analysis for them is arbitrary. There exist known scheduling algorithms. The one given in [lOI is based on the one-way model and is not suitable to circuits with feedbacks. The algorithms given in [11,13] remove this restriction, but sometimes an incorrect scheduling may be given as reported in [12]. The algorithm given in [12] can only treat the subcircuits with no more than one feedback path. The algorithm used in PLMAP has no such restrictions. The scheduling processis implemented in a simple directed graph, G = (V, E), which is used to describe the connections among the subcircuits. Each vertex in v corresponds to a compound subcircuit, and if a node is common to both subcircuits i andj and it is an internal node of subcircuit i and a basic input node of subcircuitj, then there is an edge ( ui, vi) in G. In this case,uj is called a fanout of ui and ui is called a fanin of ui. If there are ideal signal sourcesin the circuit, then a special vertex u,, is added, which corresponds to a virtual subcircuit containing all signal sources,and for each subcircuit k with a signal source input there is an edge (u,,, uk) in G. The introduction of u,, and the associatededges will make the subcircuits with signal sourceinputs be ordered as early as possible. Note that if all the nodes common to subcircuit i and j are internal nodes of subcircuit i but feedback input nodes of subcircuit j, edge ( ui, uj) does not exist in G, i.e., in the scheduling process,the feedback due to the coupling transistors is ignored. The scheduling algorithm consists of two processes.Processone is an ordering process,which is used to assignorders to as many vertices as possible. Processtwo is a feedback treatment process, which is used to find feedback paths and cut them. The ordering processis basedon the principle that if a vertex in G has no input edgesor all its fanins have been ordered, then the vertex can be ordered next. In order to check if this condition is satisfied, for each v E V, a function in(u) is defined which is the number of unordered fanins of u and is initially set to be its indegree. Let Sl be the set of unordered vertices which is initially the same as V. Let SO be a subset of Sl with in(u) = 0 for each of its vertices u. Let S2 be a subset of Sl each of whose elementsu is the fanout of at least one ordered vertex and has nonzero in(u). S2 is initially an empty set. If there are no feedback paths in the circuit or there are ideal signal sourcesin it, then initially SO is not empty and the ordering processcan be started. Select any vertex ui from SO and order it. After a vertex is ordered, searchits fanouts. For each of its fanouts ui, let in( uj) be replaced by in( uj) - 1. If in( uj) becomes zero, add uj to SO to order it later; otherwise, add%itto S2 for the possible feedback treatment process (see below).
38
Qingjian
Yu. 0. Wing / PLMA P - A MOS analysis program
This processwill be repeateduntil SO is empty. Let fanout( u) be the set of fanouts of u. The ordering processcan be described as follows: for each y E SO do begin order it, and remove it from for each y E fanout( 4) do begin replace in(v) by in(y) if in(y) = 0 then begin add it to SO; if y E S2, remove it that the indegree prior to this step*) end else if y $ S2 add it to S2; end; end
SO and Si;
1;
from S2; (*It is possible of 5 is larger than 1 and v E S2
If there are no feedback paths in G, all the vertices will have been ordered in this way. When there are feedback paths in the circuit, some of the vertices will remain unordered when SO is empty. In this case,the feedback treatment processwill be started to find the loop(s) in G and cut the feedback path(s). This process starts with selecting a starting point u, as follows. Let us first consider the case that S2 is not empty. In this case, any path going from an ordered vertex to unordered vertices must pass through a vertex in S2 first. Therefore, any vertex in S2 can be selected as the starting point. After u, is selected,a depth-first searchis implemented to examine a path emanating from it. When the search goes forward, the path is extended. When the search terminates at a vertex ui, i.e., if ui has no out edgesor all of them have been traversed,it goes back to the predecessorof ui on the path and continues to find a new path from
e4
a Fig. 5.
Qingjiarr
Yu, 0. Wing / PLMAP
- A MOS analysis program
39
this vertex. If the search goes forward from ui to vi and vi has already been on the path being constructed, then a loop is found. In this case, ui is a descendant of both U, and uj. Therefore, edge ( ui, vi) can be regarded as a feedback path to uj and will be cut. After that, the search goes back to ui and continues. This process will be repeated until it terminates at the starting point. For example, for the subgraph of G shown in Fig. 5(a), where all the ordered vertices together with their out edges have been removed and vertices 1 and 4 are in S2, if vertex 1 is selected as the starting point, then the search process may be as follows: vertex
el
e2
1+ 2 + l-2
path
cut e3
+
3 2
2
l-2-3
l-2-3-2
3+
back
2+123+1
back
l-2-3
back
back
l-2
l-3
In this example, when the search goes forward from vertex 1 through e4 to vertex 3, although this vertex has been examined before, it was not on the path being constructed and no loop is found in this case. When the search terminates at the starting point, if there are other vertices in S2 which have not been examined, then select a new starting point from them and begin a new search. In the example shown in Fig. 5(a), vertex 4 will be selected, and a new search process may be as follows. vertex
42 5 2 6 2 4-5
path vertex
4-5-6
e6 + 5 + 4 back
pa&
After changed Note feedback ordering chosen incorrect
4-5
4-5-4
back
4-5-6-l cut e6
4-5-6
4+
cut e8
back
4-5-6-4
6+ 4-5-6
e10
+5
-,4*2+4
back
back
4-5
e13
Fig. 6.
+
back
4-2
the feedback treatment process, the subgraph shown in Fig. 5(a) will be into that shown in Fig. 5(b). that, in this example, for the two’vertices 1 and 4 in S2, vertex 1 has no path to itself, and it is a descendant of vertex 4. Therefore, when the process cannot continue, if any of the vertices in S2, e.g., vertex 1, is to be ordered next as will be done by some algorithm (see [13]), an scheduling may be obtained.
e4
3
e8
-+6
1
6
40
Qingjian
Yu, 0, Wing / PLMAP
- A MOS analysis program
There may be some other cases. One is that S2 is empty, e.g., in the graph of an oscillator circuit. The other is that after all vertices in S2 have been examined, there are still some unexamined vertices in Sl, e.g., in the subgraph shown in Fig. 6, where after vertices 1 and 4 are selected as starting points and vertices l-6 have been examined, vertices 7-9 remain unexamined. Let S4 be the set of unexamined vertices during the feedback treatment process which is initially the same as Sl. Then, in the former case, there must be some loops within the vertices in S4, and in the latter case, such loops may exist. In searching the loops within these vertices, any of them can be selected as a starting point. This is because in these cases none of the vertices in S4 is a descendent of any ordered vertices. Now let us describe the process formally by using the following definitions. For each vertex u in Sl, let SW(U) be the set of unexamined fanouts of u, which is initially the same as fanout( u). Let P(u) be a path currently being constructed from the starting point u, to vertex u and S3 be the set of vertices on P(u). For each vertex uj in S3 other than u,, let pre(u;) be its predecessor on P(u), and let pre(u,) = 0. Then the feedback treatment process consists of the following steps. Algorithm Step 1. If S4 is empty, i.e., if all the vertices in Sl have been examined and all the feedback paths have been cut, go to implement the ordering process; otherwise, go to Step 2. Step 2. Select a starting point u,. Step 2a. If S2 is empty, go to Step 2b; otherwise, select a IJ, from S2 and remove it from the set. Check if it is in S4. If it is, u, is selected as a starting point and the process goes to Step 2c; otherwise, repeat Step 2a. Step 2b. Select a u, from S4 as a starting point and go to Step 2c. Step 2c. Remove u, from S4 and set pm(q) to zero. Let u, be the first vertex in S3 and ui be u,. Step 3. Decide to go forward or backward. If suc(ui) is not empty, select a ui and remove it from suc(uj), and go to step 4; otherwise, go to Step 6. Step 4. Check for feedback path. If ui $ S3, then extend P( ui) to P( ui) by adding ui to S3, set pre( uj) = ui and remove it from S4 (if it is in S4), * replace ui by uj and go to Step 3; otherwise, go to Step 5. Step 5. Cut feedback path. Cut ( ui, uj) by removing uj from fanout (vi) and let in( uj) be replaced by in( uj) - 1. If in( uj) = 0, add uj to SO. Go to Step 3. Step 6. Backtrack. Remove ui from S3, and let ui be replaced by pre(ui). If ui becomes zero, then go to Step 1; otherwise, go to Step 3. * Sometimes, an examined vertex, which has been removed from S4, may be met again, e.g., vertex 3 in Fig. 5(a). As mentioned in the above steps, after all the vertices in Sl have been examined and ail the feedback paths have been cut, the ordering process will be implemented again (Step 1) and all the remaining vertices wilI be ordered. Using the algorithm just described, if there are no feedback paths in G, each edge wilI be traversed forward only once; otherwise, it will be’ traversed at most
Qingjian
Yu, 0. Wing / PLMAP
- A MOS anabsis program
twice: one in the ordering process, and the other in process. Therefore, the traversal and its corresponding time. The other operations, e.g., the construction of S2 starting points, etc., take 0( IV/) time. Because IV1 < IEl + the time complexity of this algorithm is O(lEl).
5. Organization
41
the feedback treatment operations take 0( IED and S4, the selection of 1 for a connected graph,
of PLMAP
PLMAP is written in FORTRAN. It consists of the following processors each of which consists of a number of subroutines: (1) input circuit processor, (3) scheduler, (2) partitioner, (4) analyzer. The input circuit processor reads the description of the circuit from a specified file and stores it in a compact form in the internal array. The input is free format. Both the element level and subcircuit level descriptions are allowed. The program only requires the user to identify the pass transistors in the circuit in order to make an efficient partition. If there are subcircuit descriptions in the input data, the processor will expand the subcircuit descriptions into element descriptions. The partitioner consists of two parts: PART 1 and PART 2. PART 1 decomposes the whole circuit into first level subcircuits by using the algorithm of Section 3. After all the first level subcircuits are ordered by the scheduler, PART 2 will check for each first level subcircuit and do a second level decomposition when necesSWThe scheduler assigns the analysis order to each first level subcircuit and checks for possible strong feedback pairs. Two vertices u; and uj in G are called a strong feedback pair if arcs ( ui, uj) and ( uj, vi) exist, and if the connections between the corresponding subcircuits i andj are through the gate nodes of some transistors in them. The scheduler will combine each strong feedback pair into one subcircuit in order to speed up the convergence of the Gauss-Seidel iteration for the whole circuit. The analyzer consists of modules which implement the analysis of subcircuits and the Gauss-Seidel iteration. It sets up the subcircuit equations at the first time point and updates them later when necessary. It solves these equations and does Katzenelson iteration. After the Gauss-Seidel iteration converges, it checks the local truncation errors of node voltages and updates the time step size and the order of integration when necessary. The circuit equations used in this program are nodal equations. Backward-Euler and trapezoidal methods are used to do integration. The step size and order of integration are adjusted automatically. Because the number of internal nodes in a subcircuit is small, sparse matrix techniques are unnecessary and not used. For each subcircuit with a total of n, nodes (including n internal nodes and ni input nodes), an n,-dimensional node voltage vector and three n-dimensional vectors are constructed. These are: a nodal capacitance constant charge vector whose elements are made of Q,, &, etc., given in Appendix C; a nodal capacitance
42
Qingjian
Yu, 0. Wing / PLiUAP
-
A MOS
analysis
program
charge vector, and a nodal capacitance current vector. Also, for each subcircuit, two n * n, matrices are constructed: a nodal admittance matrix and a nodal capacitance matrix. These matrices and vectors make the update of the nodal equations easy.
6. Examples Four ring oscillators with different number of inverters and simulation time have been tested with PLMAP. The CPU time is shown in Table 3. It can be seen that for the cases tested the CPU time is linearly proportional to the product of the number of inverters and simulation time. Compared to SPICE, for the 7-inverter-oscillator, this program is eight times faster, and when the E-inverteroscillator was tested, SPICE suffered from convergence problem.
7. Conclusion We have described a circuit analysis program for simulation of LSI MOS circuits. It can be seen that the piecewise linear model and the 2-level relaxation method are useful techniques for fast and accurate simulation of LSI MOS circuits. The piecewise linear model used in this program does not include the short channel effect on MOS transistors. When this effect must be considered, multi-variable functions will appear in both the drain current and charge equations [9]. In this case, some reasonable simplifications should be made to avoid making the model too complicated. More research should be done in this respect.
Appendix
A. The piecewise linear approximation
of function
r( x, y) = xy
Consider a square region centering at (xi, uj) with 4 vertices (xi + a, yj), (xi - a, rj), (xi, rj + a) and (xi, uj - a) as shown in Fig. A.l. The first order Taylor series expansion of T(X, y) at point (xi, yj) is r&y)
=yjx
+ xiy - xiyj.
Table 3 Number 3 7 15 31
of inverters
Simulation 60 120 250 500
time (nS)
CPU time (second) 10 42 185 731
Qingiian
Yu, 0.
Wing
/
PLMAP
-
A MOS
analysis
program
43
Let e(x,y)
= [rD(xyy)
-~(x,Y)]~;
then and
e(Xi,Yj)=e(Xi,Y)=e(X,Yj)=O
ih/i3x
z 0,
lk/ay
+ 0
when e # 0,
(&Y)ED.
(A4 (A.2)
From (A.l) and (A.2) it is known that the maximum value of e(x, y) can only be obtained on the boundary of D. It is easy to show that e max= ($2)’ and (10) is obtained. From (A.l) it is known that, at the 4 vertices of the square, rD(x, y) = r(x, y). It is well known that the values of a 2-variable linear function on a straight line are fully determined by the values at any two points on the line. If region D’, which is the neighborhood of D, is also a square centering at (xi + a, yj f a) and has 2 vertices common to D, then rD(x, y) and rD( x, y) defined in D' and D respectively will be continuous on the boundary common to the two regions. This proves the continuity of the piecewise linear approximation. Choose xi = y, = 0; then region D can be expressed by (8) and (9). Given x and y, xi and yj can be determined in the following way. Let u=x+y, urn = xi +yj,
u=x-y
(A-3)
(A.41
un = xi - y,;
Y
(‘i ,Yj+a)
(Xi-a,Yj)
(Xi+a,YJ)
(Xi ,Yj-a) Fig. A.l.
44
Qingjian
Yu. 0.
Wing
/
PLMA
P -
A MOS
analysis
program
then D can be expressed by IU--UU,I
Iu - u,l< a
(A.5)
and v,=2na, m,n=O, 1,2, . . . . urn = 2ma, (A.6) Given x and y, from (A.3), (A.5) and (A.6), u,,, and u,, can be determined, and xi = 0.5( u, + !I,),
Appendix
y, = 0.5( u,, - u,).
B. Charge and current equations
(A-7)
of a p-n junction
There are two p-n junctions in a MOS transistor. Let I, V, C and Q be the current, reversed voltage, capacitance and charge of a p-n junction. We have I=
Is(e-4v/k’
- 1)
where Z, is the reversed saturation
(B.1) current,
and
where C, is the value of C when V = 0, I& is the built-in potential of the junction, and n is a constant (n = f or f). The expression of Q can be derived from (B.2) by integration of C with respect to V: V
Q=/W+Q,.
(B-3)
6, where V0 and Q,, can be chosen arbitrarily since we are now only interested capacitance current dQ/df. If we choose V, = - VP” and Q0 = 0, then
Q=
Co VP”L, ( v/ VP”) 9
in the (B.4)
where f&x)=
+1 + X)‘-n l-n i (1 + x)/0.1”,
- &O.P”,
x > -0.9, (B.5) x<
-0.9.
Qingiian
Appendix
Yu, 0. Wing / PLMAP-
A MOS analysis program
45
C. Piecewise linear model for MOS transistors
There are 11 real variables which determine linearly modelled MOS transistor. They are
x2=
v*-
x:,=v
v,,
ivJ
x4=
-v,,
X8
=
K&n,
x9
=
Vdb/
&+x2’
x5=
V&-X2r
x~ =
Vsb,
i X 11 =
7
1
4
region of a piecewise
VP” 1
x3/x* X 10 =
the operating
x3
’
W)
0,
x* =x3 = 0,
Vgb -
4b)/Y2~
where x4 and x5 correspond to u and v in (A.3), and x6 and x7 are used to determine the operating regions of the currents of the two p-n junctions. The following nonlinear functions appear in the transistor model used in PLMAP and each is approximated by a piecewise linear function defined below, where the coefficients depend on the operating region of the transistor. JV,,+=la,(K,+2+F)+b,9
W-2)
(v,-
~)‘~a&$-
19+b2,
(C.3)
(Vgd-
<)2Aa,(v,,-
v,)+b3,
(C.4
v,,(v,-
v,)*a,(v,-
<)+b‘&-=,b,,
(C.5)
(C.6)
r,(v,,/v,“)~=~(v,,/v,“)+b,, I,(v,,/v,n) &s
A -
v,)&d
a6(v,t&n) -
+
v,)/&,-
(C.7)
b6, v,+
sd-
v,+7(l/gd-
<)+b7(v&-
y), (C.8)
1+ 4(v,, The equations
- Vfb)/Y 2 i 4a,( V+ - I&)/y2 of the piecewise-linear
+ b,.
W.9)
model are as follows.
C. 1. Drain current equations C. I. 1. Linear region (x2 > 0, X3 2 0) I,=G,(V,-Y~,&,)+~~-~&~~-Y~&,)-~~~
(C.10)
where G, = ka,,
G, = ka3,
J2=k{b2-=2[V~+2~F+Y(=,2~~+b,)]}, J3=k{b3-a3[V~+29~+~(a,2~,+b,)1j.
(C.11)
46
Qingiian
C.I.2. Saturation
Yu, 0. Wing / PLMAP
- A MOS analysis program
region (x2 > 0, x3 < 0) (C.12)
Id=G2(Vgs-ya,Vsb)+J2+G~Vds-J4 where G, = k[(l
- a)a2 + aa,], (C.13)
J2=k(l-a){b2-a2[~,+29,+~(a,2~,+b,)]}, G4 = akb4, C.2. Intrinsic
J4 = aka,b,.
charge equations
The general form of the intrinsic
pii=\;
is
charge equations
;
is as follows:
(C.14)
;-[;+[;;-y
where c,,
= 0,
Cgd = -(cs,
+ Cd, + GA + Cd, + CtJtJY
c&g= -(csg+
cd,+
c,,L
Cgb = -(Cs,
cgs= -(cs,+
Gs+
GA
Q,,
C.2.1. Accumulation cpp=
(Qso+ Qdo+ Qd
region (x-, < 0, x,, < 0)
-Cbg=
Other parameters
= -
(C.15)
-C.& = c,, = c,,
Q,, =
-
Qm =
- GG.
(C.16)
in (C.14) are zero.
C.2.2. Subthreshold
region (x2 < 0, x,, > 0)
egg = - C& = - C,, = C,, = 2a,C,, Q,, = - Qbo = [0.5y2(b, Other parameters C. 2.3. Saturation
- 1) - 2a,&]
(C.17)
C,.
in (C.14) are zero. region
Cdg = Cd, = Cd, = Cd, = C,, = C,, = 0,
Qbo=
Qdo = 0,
C,, = - C,, = ya,C,, csg= -$c,,
-Gb~2~F
+ b,),
C,, =
C,, = - h&,
Q,,=,C~[v,+2~,+y(a,2~q+bl)l.
fC,(l + yq),
(C.18)
Qingiian
Yu, 0. Wing / PLMAP
- A MOS analysis program
47
C.2.4. Linear region Csg= -fC,(l+a,+b,), Css=~[(2+b7)+ya,(l+a7+b7)]C,, C,, = +(a, - l>C,, Csb = v,G,~ Q,, = iCo(l
+ a7 + b,)[G
+ 29, + yW%+
+ bA19
(C.19)
Cdg = (a7 + b, - l>C,, Cds= [v,(l -a7-b7)--b71G7 Cd, = (1 - a,)C,, Cd, = ya&,,
~~~=(l-a,-b,)[V,+2~,+y(a,2~,+b,)lC,. Other parameters
are the same as in the saturation
region.
C.3. p-n junction chargeequations (C.20) where C sbsb The equations
=
a,CsbO,
Qsbo
=
b5C,bCIV,n-
for Qdb are similar.
C.4. p-n junction current equations
0, Isb
Kb
2
-
Vdiode~
Kb
<
-
GocIe~
(C.21)
= G(J”h
+
Vdiocte)~
where Vdiode= 0.7~ and G = 0.1 mho. The equation
for I& is similar.
References [l] Nagel, L.W., SPICEZ: A computer program to simulate semiconductor circuits, ERL Memo ERL-M520, Univ. of California, Berkeley, 1975. [2] Weeks, W.T., A.J. Jiminez, Q.W. Mahoney, D. Matha, H. Qassemzadeh and T.R. Scott, Algorithms for ASTAP - a network analysis program, IEEE Trans. Circuit Theory CT-20 (6) (1973) 628-634. [3] Chawla, B.R., H.K. Gummel and P. Kozak, MOTIFS - A MOS timing simulator, IEEE Trans. CAS CAS-22 (1975) 901-910.
48
Qingjian
Yu, 0. Wing / PLMA P -A
MOS analysis program
[4) Fan, S.P., M.Y. Hsueh, A.R. Newton and D.O. Pederson, MOTIS-C: A new circuit simulator for MOS LSI circuits, Proc. 1977 IEEE Internal. Symp. on Circuits and Systems (1977) pp. 700-703. [5] Newton, A.R., Techniques for the simulation of large-scale integrated circuits, IEEE Trans. CAS CAS-26 (9) (1979) 741-749. [6] Jouppi, N., Timing analysis for NMOS VLSI, Proc. 20th Design Automation Cant (1983) pp. 411-418. [7] Lelarasmee, E., A.E. Ruehli and A.L. Sangiovanni-Vincentelli, The waveform relaxation method for time domain analysis of large scale integrated circuits, IEEE Trans. Cohputer Aided Design CAD-l (3) (1982) 131-145. [8] White, J. and A.L. Sangiovanni-Vincentelli, RELAXZ: A modified waveform relaxation approach to the simulation of MOS digital circuits, Proc. 1983 IEEE Internat. Symp. on Circuits and Systems (1983) pp. 756-759. [9] Yang, P. and P.K. Chatterjee, SPICE modeling for small geometry MOSFET circuits, IEEE Trans. CAD CAD-1 (4) (1982) 169-182. [lo] Ruehli, A.E., A.L. Sangiovanni-Vincentelli and G. Rabbat, Time analysis of large scale circuits containing one-way macromodels, IEEE Trans. CAS CAS-29 (3) (1982) 185-189. [ll] Tanabe, N., H. Nakamura and K. Kawakita, MOSTAP: A MOS circuit simulator for LSI circuits, Proc. 1980 IEEE Internat. Symp. on Circuits and Systems (1980) pp. 1035-1038. [12] Wei, Y.P., Large-scale circuit simulation, Ph.D. Thesis, R-977, UILU-ENG, 82-2243, Urbana Illinois, 1982. [13] Lelarasmee, E. and A. Sangiovanni-Vincentelli, RELAX: A new circuit simulator for large scale MOS integrated circuits, Proc. 19th Design Automation Conf: (1982) pp. 682-690.