We present a new procedure to invert for kinematic source parameters on a finite fault. On the basis of the reciprocity relation of the Green's functions, we use a newly developed fourth-order viscoelastic finite-difference algorithm to calculate three-dimensional (3-D) Green's functions (actually the tractions) on the fault. We invert the data for the unknown source parameters at the nodes (or corners) of the subfaults. The source parameters within a subfault area are allowed to vary; this variation is calculated through bilinear interpolation of the four nodal quantities. We have developed a global nonlinear inversion algorithm that is based on simulated annealing methods to solve efficiently for the nodal parameters. We apply this method to the 1989 Loma Prieta, California, M 6.9 earthquake for both a 1-D and 3-D velocity structure. We show (1) the bilinear interpolation technique reduces the dependence of inversion results on the subfault size by naturally including the effects of nearby subfaults. (2) While the number of synthetic seismograms that must be computed is greatly increased by the bilinear interpolation, the structure of the inversion method minimizes the actual numbers of computations. (3) As expected, complexity in the velocity structure is mapped into the source parameters that describe the rupture process; there are significant differences between faulting models derived from 1-D and 3-D structural models.