003x-1101/x5 $3.00 + .m % 19X5 Pergamon Press Ltd
INVARIANT IMBEDDING IN SEMICONDUCTOR DEVICE SIMULATION P. T. LAI and Y. C. CHENG Department of Electrical Engineering, University of Hong Kong, Hong Kong (Recemed 13 March 1984; in revisedform 18 Muy 1984) Abstract-A general, fail-safe and fool-proof numerical algorithm called invariant imbedding is used to successfully solve both the Poisson equation and the continuity equations for two-dimensional (2-D) problems. This direct method avoids the convergence problem commonly encountered in iterative methods, and is computationally more efficient than the classical Gaussian elimination method especially when the insulator occupies a large portion of the semiconductor device under simulation, and/or more accurate treatments are made for the interfacial region.
1. INTRODUCTION In the simulation of semiconductor devices, one normally has to cope with the self-consistent solution to a system of three coupled nonlinear partial differential equations describing the intrinsic physics of the devices. The discretized and linearized version of such a system results in one to three linear matrix equation(s) of very high rank, whose solution requires excessive computation time. Iterative techniques such as the classical iterative methods and the iterative factorization methods [l] are quite often used to solve a large system of linear equations owing to their small memory space required, and probable faster speed especially for problems of very large size. But they may show a very slow convergence rate in case of non-optimal parameters estimated, or even divergence, especially in solving the continuity equation since the iteration matrix involved is not that well conditioned numerically. One way to get around this problem is to use a direct method for the continuity equation while the iterative method is applied for the Poisson equation [2]. However, the use of two algorithms increases both the programming effort and the required storage area. In fact, the convergence problem becomes more pronounced for the Poisson equation also if the device structure has a nonplanar insulator-semiconductor interface. In that case, the use of singular functions to avoid excessively fine mesh in the vicinity of the corner(s) along the interface deteriorates the conditioning of the iteration matrix[3,4]. Last but not least, iterative methods have the possibility of premature termination of the iteration cycle and convergence to a wrong solution. Therefore, in order to increase its strength and widen its area of application, the simulation program may have to resort to direct methods all the way through. The most simple and easy available direct method of course is the Gaussian elimination method. Nevertheless, it is generally criticized for its inefficiency both in terms of speed and storage, for problems of large size.
435
It is the aim of this work to study the feasibility of another direct method called invariant imbedding[5], with a view to alleviating the problems of the conventional elimination method while preserving its easy implementation. 2. NUMERICAL SOLUTION to find solutions to the Poisson equation and the two continuity equations in a self-consistent way, each of these nonlinear partial differential equations is first Newton-linearized to a quasi-harmonic equation of the form In order
(1) for 2-D cases, with f(x, y) representing either the electrostatic potential u, the exponent of the negative of the electron quasi-Fermi potential e-l‘ and the exponent of the hole quasi-Fermi potential eW. In case of the breakdown of the Einstein relation, u and w have to be muItiplied by the ratio of the mobility to the diffusion coefficient[6]. Functions r, s, z in eqn (1) can be easily related back to the original nonlinear partial differential equations. The region being analysed is then discretized by imposing a tensor-product grid upon it. By applying the standard five-point difference formula to a node of row index i and column index j, the discrete equation for the unknown f,, at this mesh point can be written as
c,,fi,
-
c,~I,fr~l,-c,+I,f’+I,-ci,-lfr,-l
-c,,+~fr,+~
=e,,
for some coefficients c’s and e,,. Doing the same for a whole column knowns f,=(r,,,r,,....,r,,)’
(21
of M un-
(3)
436
P. T. LAI and Y. C. CHHNG
and assembling
the M equations
together
gives
Q,f,-L,f,~,-R,f,,,=e,,
(4)
I
j =
.
Y ..-10
*
I
1 ... I I
x
where capital letters are matrices and bold letters are vectors, all of dimension M. Note that M may vary for each column in case of irregular structures, To clarify the picture, we consider a simple problem in which a rectangular region with constant Dirichlet boundary condition f, and uniform mesh spacings is governed by the Laplace equation. In this case, the matrices and vectors in eqn (4) are given by
%,=
-1 i 0
A ,+lf,
e, =
/=l,
f, i 0
otherwise.
that we seek a solution
j
+b,+lt
Fig. 1. An MOS structure chosen for numerical simulation. Bold boundaries have Dirichlet conditions and the rcmaining parts arc of Neuman type. The arrows below the region show the procedure used in solving for the unknowns f,.
I=M,
=
O-N (left-right) (6)
by introducing some intermediate unknown matrices A, and vectors b, which are to be found as follows. By substitution of eqn (6) eqn (4) can be transformed to f,=(Q,-
R,A,+l)ml(L,f,m,
Viewing this equation by j - 1 gives
+ R,b,+,
+e,>.
as eqn (6) but with j replaced
b,=(Q,-R,A,+,)~l(R,b,+,+e,).
(7)
These matrix and vector recurrence relations combined with eqn (6) are sufficient for evaluating all the unknowns f,,. The procedure involved is demonstrated by the arrows underneath the silicon region in Fig. 1. Supposing the right boundary condition is given by f,< ,r and making use of eqn (6) results in A N + I = 0,
b ,%+I =f,v+1
I
(5)
Invariant imbedding means to eqn (4) in the form of =
:/
otherwise,
f,.
f J+l
for fixed boundary
Then successive applications of eqn (7) provide us with A,., b,v,. . , A,, b,, respectively, and can bc viewed as a form of block Gaussian elimination. If the left boundary condition is f,,, then all the unknowns f, ,f,, ,f, can be evaluated in turn by repeated use of eqn (6). Thus the original boundary value problem can be replaced by an initial value problem through the idea of invariant imbedding with eqn (8) as the initial conditions, or it may be said that the 2-D problem is broken down into a sequence of 1-D problems. For the more general case of irregular region and/or mixed boundary condition(s), a detailed discussion can be found in Ref. [5]. It should be noted that the use of the standard five-point difference scheme together with the transformation of variable given by eqn (6) is efficient only for the case of the Poisson problem with grid points equally spaced in each column. In such a situation, L, and R, are proportional to an identity matrix and the Q, are symmetric, rendering all the A, symmetric, as can be deduced from eqns (7) and (8). This can reduce both the computation time and storage by approximately half as compared to the case of asymmetric A,. For more general cases, in order to retain the symmetry of A, (if possible), it is necessary to discretize eqn (1) using the standard box integration method and to replace eqn (6) with a different transformation
(8)
(Note
that eqn (9a) has to be replaced
A ,v+, = I,
L,f, for free boundary
condition.
( 9a)
R,f,+i=A,,if,+b,,,.
condition
or
bNi 1 = 0
SUBSTRAI'E
Pml=l,
L, = R, = I,
e, = (e,),
:;i
SiO
I=m,
4 Q,=(a,,),
GATE
if the direction
I =A,
of marching
,f, +b,
by
I
is reversed,
(Yb) as in the
Invariant imbedding in semiconductor device simulation oxide region of Fig. 1.) These equations define a new set of auxiliary unknowns A, and b,, which can be calculated by the new recurrence relations
be calculated first. Using a standard box integration for the interfacial column (j = 0), eqn (4) can be written as Q,f,
b, = R,. l(Q, -A,+l)-l(b,+l
+e,>,
431
Substitution
(10)
- L,f_,
- R,f,
= e,.
of eqns (9a) and (9b) with j = 0 gives
and the new initial condition or
(13)
A Nil = 0, b N * 1 = R,vfr++l
for fixed boundary
f,=(Q,-A_,-A,)-‘(e,+b_,
condition
or
With f, at hand, we can start marching back through both regions, collecting all the unknowns f, by making use of eqns (9a) and (9b).
(11) A .v+ L =&, b N+ L = 0
for free boundary
condition.
3. DISCUSSION
For self-adjoint differential operators such as the one in eqn (1) with f(x, y) as the choice of dependent variable, Q, are always symmetric and R,_ 1 and L, are related by R ,_I = LT.
A 2-D analysis program was developed to implement the numerical scheme described in the previous section. It was checked for accuracy against another program using the one-step Newton-SLOR method[2]; the results showed good agreement. No numerical instability was encountered during the matrix inversions in eqn (7) or eqn (10). Detailed discussion on stability can be found in Ref.[5]. Its validity is further enhanced from the calculated surface potential which is obtained by covering the whole left boundary of Fig. 1 with gate metal. (The results in Table 1 show little difference between the numerical values and the analytical values (two times quasi-Fermi voltage).) Such a quasi-l-D version was applied in the study of the influence of nonuniform substrate doping on the C - V curve of MOS structures. A typical graph is illustrated in Fig. 2. This program has also been successfully used in modeling the narrow width effect for MOSFETs[7], thus gaining further credibility. Some simulation results are presented in Fig. 3. For the purpose of comparison among numerical methods, the rectangular device structure shown in
(12)
By inspection of the form of A, ~, in eqn (lo), and making use of eqns (11) and (12) we see that the symmetry of A, can be preserved. The procedure demonstrated above resembles that of marching from one boundary to the opposite boundary and then back. This applies well to the continuity equations as they only concern the semiconductor portion. On the other hand, the solution to the Poisson equation is also needed in the insulator region. As a result, marching has td be started from both the gate and the substrate towards the interface. After satisfying the Gauss law there, we then go back to the two boundaries. The whole process is clearly shown by all the arrows in Fig. 1. Before the unknowns f, in both regions can be evaluated, the unknown at the interface, i.e. f, must
Table 1. Comparison of numerical and analytical surface potentials (U,, and (/,.) under the condition of strong inversion in MOS structures with various combinations of oxide thickness (rox), substrate doping (N,) and back-bias ( Vss)
1 a
+b,).
( 500 1 lFl51
0
I.5854l.58581
.07$
1
P. T. LAI and Y. C. CHENG
43x
CAPACITANCE/(32.4 nF cm-z)
DOPING/m
-3
DISTA?JCS "I:OH INTERFACE
0 ').7.m
1
L’+__
I
-0.92
-1.50
lo.33
I
1.2:
1 1 i,!i
I
0.83
1.4,
GAT?. VO!.':AI:E: (Y) Fig. 2. Simulated CV curves illustrating the effect of oxidation-induced surface dopant depletion in MOS structures with oxide thickness 1000 A and substrate doping 1 X lOI cm ‘. The inset shows the surface-depleted profile which assumes an erfc shape.
1.60
.---.
‘\,+ -3-t:
--2+ - -._
0.00
=
11
v
+++++experimental
Vtrc. = 11 V L1
x---x numerical
vt;:; =
3 V
ooooo experimental
Vjc
? V
=
2.00
++
-:a. %-__ --...__+o--o - - - - - _ _+?+ --_ I
00
VBS
++
\+
CI.
numerical
4.00
1
6.00
I
I
8.00
10 'dlJ
CHANNEL WIDTH ( m) Fig. 3. Comparison of numerical and experimental threshold voltage changes of conventional MOSFETs vs channel width for two back-biases Vss. Oxide thicknesses for gate and field regions are 500 and X500 A. Substrate doping has a uniform value of 1 X lOI5 cm-‘. (Experimental data arc extracted from Ref. [9].)
Invariant
imbedding
in semiconductor
Fig. 1 is adopted. In practice, the structure is simply a metal pad isolated from the semiconductor underneath which can be used to study the effect of fringe capacitance. The gate voltage is set equal to the threshold voltage of a large device with no backbias. Other parameters are as follows: substrate doping (5E15 cme3), oxide thickness (500 A), and gate width (1 pm). Such a choice is arbitrary because it has little effect on the computing effort by the direct methods. The computation time required to solve a system of N X N linear equations Gx = c using two very simple direct methods is shown in Table 2. The values listed are the CPU time in seconds on SPERRY (UNIVAC) 1100. The classical elimination method adopts a standard band solver subroutine while invariant imbedding makes use of another standard subroutine for the inversion of a full matrix (eqn (7) or eqn (10)). As for all the conventional direct methods, these two methods also require 0( N4) arithmetic operations. Note that the simulation time is assumed to be dominated by the matrix manipulations. 3.1. Poisson equation The matrix G obtained by applying the box integration method to the Poisson equation is usually symmetric. But this symmetry is lost if more special treatments are made at the interface, e.g. Ref.[8] or singular functions are used in the vicinity of comers (where the electric field is unbounded) [3,4], in order to prevent prohibitively small mesh size there. These treatments also may approximately double the bandwidth, rendering the band solver routine very inefficient due to about a four-fold increase in computation work. However, for invariant imbedding, the matrices A, in eqn (7) or eqn (10) do not depend on what happens at the interface because they are calculated from the left and the right boundaries towards the interface. As viewed from eqn (7) their symmetry is
device simulation
hooked to that of Q, since both LJ and R, are proportional to an identity matrix for the Poisson equation. It can be easily shown that Q, and therefore A, are symmetric provided that the rows of unknowns are spaced equally apart. This can be the case where the lateral electric field in the device is not very strong, e.g. in the simulation of the width cross section of a conventional MOSFET. For meshes nonuniform in both x- and y-directions, symmetry in A, can always be assured if eqn (10) is used instead of eqn (7). In fact, it can be easily shown that symmetric global matrix G implies symmetric A, but that the reverse is not always true. Consider the most promising case for invariant imbedding in which an asymmetric G is used for the band solver method. As seen from Table 2, the former is approximately four times faster than the latter. This can be further enhanced due to the fact that A, for the insulator region needs to be calculated only once during the Newton iterations. For example, if about half of the nodes are allocated on the insulator side of the interface, e.g. structures with very thick field oxide, the solution time for each iteration may be reduced by a factor of two approximately. Multiplying out, invariant imbedding in such a case can have a speed gain of about eight over classical elimination. Furthermore, the irregular structure of a given region is already taken into account during the evaluations of A,. This is not the case for classical elimination because the resulting variable bandwidth of G is not considered by the standard band solver routine. We now consider the least favorable situation for invariant imbedding, in which it uses symmetric A, while symmetry is also present in G of the other method. Invariant imbedding is a bit faster as shown in Table 2, and this can be further enhanced by its two advantages just mentioned. In addition classical elimination is less accurate here because it cannot use special treatments for both the interface and the comers if G is required to be symmetric in this case.
Table 2. CPU time in seconds on SPERRY 1100 for each solution of a linear matrix equation of N x N unknowns resulting from the discretization and linearization of the semiconductor equations governing the region shown in Fig. 1
N
CLASSICAL ELIMINATION
439
440
P. T. LAI and Y. C.
3.2. Continuity equation As for the Poisson equation, the global matrix G is normally symmetric with em” or e” as variables since the differential operator in eqn (1) is of a self-adjoint form. But the coefficients of G then contain factors like e“ or e ’ which can very easily underflow or overflow numerically. This scaling problem can be mitigated a bit while retaining the symmetry of G by introducing peculiar forms of unknowns such as err/2m” and e” ‘/I, which are less volatile[6]. But to be more safe and general, this numerical instability is best suppressed by taking the electron concentration n and the hole concentration p as the variables. The resulting matrix G then loses its symmetry. Asymmetry in A, also arises for invariant imbedding because L, and R, in eqn (7) arc not proportional to an identity matrix. In fact the symmetry of A, simply follows that of G if we resort to eqn (10). Therefore, if both methods involve asymmetric matrices, invariant imbedding is again two to three times faster than the band solver method. In addition, as in the Poisson case, invariant imbedding may be said to take into consideration part of the sparsity present in G by inverting matrices of different orders. It should be noted that a complicated treatment is usually not necessary at the interface and/or the corners associated because the continuity equations are solved for the semiconductor region only. 3.3. Memory requirement To simplify the comparison, a rectangular region of M rows and N columns of unknowns (M I N) is considered. From the structures of the mat&es, it can be shown easily that matrices A, occupy a total area of NM2/2 and NM’ for the symmetric and asymmetric cases, respectively, while G in the conventional elimination method requires an area of NM2 and 2NM’ correspondingly. In view of the discussions in Sections 3.1 and 3.2, and assuming the program size is dominated by these matrices, invariant imbedding may have a storage space one fourth that of the other method at the best and approximately one half in the less favorable situation. This advantage is further enhanced in case of irregular regions because the matrices A, then are of variable order. This again demonstrates that invariant imbedding exploits the sparse nature of G to a certain extent. 3.4. Other direct methods It should be emphasized that the invariant imbedding method is proposed here mainly for simplicity in understanding, ease of implementation, power, and general application. Of course, there exist many other direct algorithms which work at a higher speed and avoid very large core storage, but they usually lack some of the advantages just mentioned. One such method is the class of ordering techniques which exploits the sparsity of the matrix G
CHENG
during the elimination procedure. GEMINI [2] and FIELDAY [9] are examples of device simulators which use this method. Among such techniques the George nested dissection algorithm[lO] shows the best performance with the asymptotic bounds of 0( N’) for work and O(N’ log N) for storage which have been shown to be optimal for direct methods using symmetric Gaussian elimination. However, the construction of such a routine is by no means trivial (e.g. the involved data structures and the symbolic LU factorization) and often the best dissection has to be done by examining the original structure of a given region, thus making automation difficult, In consequence of this, George [ll], proposed a one-way dissection with greater ease of implementation but the asymptotic bounds are 0( N’5) and 0( N’“) lying between those of band Choleski and nested dissection. The main condition for the success of these dissection methods is that no row interchanges are used to preserve stability in the elimination. Therefore, definiteness or diagonal dominance of matrix G is clearly desirable. For best performance they may also require a regular region with grid size to be a power of two. An even better method which reduces work and storage to 0( N’ log N) and 0( N’), respectively. is the method of marching (also termed shooting or frontal elimination) resurrected by Roache[l2] but notorious for its instability. This method has been generalized by Bank and Rose[13] to provide better stability bounds, but its value for particular and more general elliptic equations is yet uncertain. There are also fast methods based on the fast Fourier transform technique of Hockney[l4] or the reduction method of Buneman[lS]. For example, Reiser [ 161 applied this in his transient simulation of semiconductor devices. Although both show the same asymptotic bounds as the marching method, their success depends on particular forms of the dilferential equation (e.g. separable equations only, thus inapplicable for the continuity equations). and to a large extent, of boundary conditions and especially of boundary shapes, even though they can be cxtended to work on non-rectangular domains by the capacitance matrix method[17]. Even worst is that these two algorithms do not work for the Newtonlinearized Poisson problem. 3.5. Iterutiw meth0d.s Now let us turn our discussion to the iterative methods which require small storage. Due to the well conditioning of the matrix G for the Poisson equation, execution time is usually shorter than direct methods especially for very large problems, and is cut further for nonlinear equations by using the one-step NewtonSLOR [2]. Some computational times based on the structure in Fig. 1 with the same parameters as for Table 2 are listed in Table 3. It should be borne in mind that the comparison is rather rough because the performance of the iterative method depends heavily on the biasing conditions.
Invariant
imbedding
in semiconductor
device simulation
441
Table 3. Comparison of CPU times in seconds on SPERRY 1100 for solving the Poisson equation governing the structure in Fig. 1 with N x N unknowns. Error bound is set at 1 X 10m4
N
18
22
26
30
INVARIANT ItBEDDING
7.0
14.7
27.5
47.0
ONE-STEP NEWTON-SLOR
8.0
17.8
23.9
43.5
the device structure, and other parameters. In fact the example given may represent approximately the best case condition for the iterative method because higher values of biases, irregular geometries or the bad conditioning of the continuity equation coefficient matrix can deteriorate its convergent rate by a large amount. In addition, the advantages of invariant imbedding in cases of structures with a large proportion of insulator or inclusion of singular functions at comers are not yet explored. After all, iterative methods pin their success heavily on the estimation of optimal parameter(s), values of initial guests, and the numerical conditioning of the iteration matrix concerned. Furthermore, their linear convergence rate is very slow and limits considerably the available accuracy. This very last property implies that they may not be suitable for transient simulations in which very accurate solution at each time step is required to prevent any build-up of error, and also the solution may not be improved by a deferred approach to the limit [l] since the discretization error may not dominate. Other disadvantages are the possibility of non-unique solution even if no divergence occurs, and the possibility of inapplicability when singular functions are used for comer regions where the electric field tends to infinity[3,4]. Finally, in order to be fair, it must be stressed that for sufficiently large problem size, the iterative methods, if convergent, are superior to almost all direct methods. This fact can be illustrated in Table 3 by the slower increase of solution time for the one-step Newton-SLOR method with respect to N. 3.6. Further comments Invariant imbedding, as for other direct methods which store the decomposed matrices during the LU factorization, is very efficient for the many-off cases in which a linear partial differential equation with constant coefficients is solved many times for a given structure with varying right-hand side or different boundary conditions. The reason is that all the matrices A, or G do not change during the whole course of analysis and therefore, need to be calculated only once. Some examples are the numerical modeling of the narrow width effect in MOSFET where the nonlinear carrier terms are neglected in the Poisson equation[l8] or the Picard iteration used instead of the Newton iteration, and the transient
analysis of the semiconductor device where the Poisson equation need not be linearized due to the known carrier concentrations obtained from the previous time step[16]. For other problems, invariant imbedding may show a slower speed which can still be tolerated provided the number of mesh points is not very large. The use of sophisticated formulation at the interface and the comers can therefore do a good job here. In addition, the use of the concept of pipelining during the evaluation of the matrices A, on vector machines and the advent of faster computers can further enlarge the size of the problem while keeping the program run within an acceptable time. As regards to the memory requirement, modem large computers of high computing resources can usually accommodate the storage of the matrices A,. For even larger problems, out-of-core storage can be used very efficiently and easily since A, are only required sequentially. Two final comments are as follows. In the case of 3-D device simulation where the number of unknowns is enormous, almost all the direct methods including invariant imbedding have to succumb to the iterative methods (if convergent) in terms of both speed and storage. The other comment is that invariant imedding provides another example where the finite difference method is more efficient than the finite element method if a tensor-product grid is used. The reason is that the finite element method will normally generate tri-diagonal L, and R, in eqn (7), thus rendering all A, asymmetric no matter how the global matrix G behaves. Even though this is not the case if eqn (10) is used, a tri-diagonal matrix equation (9) has to be solved in order to find the final unknowns f, 4.
CONCLUSION
The possibility of invariant imbedding in 2-D device simulation is demonstrated and discussed. This method compares favorably with the classical elimination method in terms of speed and storage, while retaining the good properties of the latter, i.e. the simplicity in understanding, the ease of program construction and the wide range of applications, some of which may be absent in other fast direct methods. It also eliminates any convergence problem accompanying any iterative method even though the latter shows less storage and may be faster.
442
l’ T L,zr and Y. C’. C‘HLNC;
AckrzoM;ledgeme,lr.s~The authors are indebted to Dr L. G. Tham for helpful discussions. This work was supported in part by the Industrial Development Board of Hong Kong. REFERENCES
1. D. Jacobs, The State o/the Art ttt ~Nunwrcul Attu!~w.\. Academic Press, New York (1977). 2. J. Greenfield, C. Price and R. Dutton. NATO Advanced Study Institute on Process and Device Simulation for MOS-VLSI Circuits Sogcsta-Urhino. Italy (19X2). 3. G. BirkholT, J. Appro.vnrcltio!t Theog~ 6, 215 (lY72). 4. L. Fox and R. Sankar, J. Imt. Math Applic,. 5. 340 (1969).
E. Angel and R. Bellman. Drvrunnc Progvomnzing urid PDE. Academic Press, New York (1972). 6. R. E. Bank, D. .I. Rose and W. Fichtner, IEEE Trcm.~ Electron Lkices ED-30. 1031 (1983). 7. P. T. Lai and Y. C. Chcng, Solid-Sr. Electron. (1984). in press.
5.
8. M. Szuhar, SoldSt. Ektron. 25, Y63 (19X2). 9. E. M. Buturla. P. E. Cottrell, B. M. Grossman and K. A. Salsburg, IRMJ. Res. Dell. 25, 21X (1981). 10. J. A. George. SlAMJ. Nunier. Atrul. 10, 345 (1973). 11. J. A George. Full Jotut Camp Cmferemc. AFIPS Press. New Jersey (lY72). 12. P. J. Roachc, I.ecture Note.\ 1,~ P/~~c~.s (Edited by M. Holt). (8th cdn.) Springer, New, York (1971). 13. R. E. Bank and D. J. Rose, SIAMJ. Mutiter. Ad. 12. 529 (1075). 14. R. W. Hocknq. .J. Assoc~. C‘omput. Mac~h. 12, 95 (1965). 15. 0. Buneman, Rep. 294. Stanford University Institute for Plasma Rescarch, Stanford, California (1969). 16. M. Rciser. (‘on~p. Method,\ ApI. Med. Ettpg. 1, 17 (1972). 17. B. L. Buzbee. F. W. Dorr, J. A. George and Ci. H. Golub. SIAM J. Numer. Awl. 8, 722 (1971). 1X. C.-R. Ji and C.-T. Sah. II:‘I:‘E Trcrns. Electron Der~~e.v ED-30. 635 (1983).