This paper presents a new parallel algorithm and performance results of iterative solution methods for three-dimensional MOSFET simulation with Gummel's method. A splitting-up operator method is proposed for incomplete factorization of sparse matrices arising from semiconductor device equations, suitable for parallel computations. This method is combined with the conjugate gradients and BiCGSTAB procedure to obtain a new parallel version of the iterative solution methods. Natural parallelism is realized by developing the solution method according to the natural ordering. In large-scale simulations of greater than 100000 grid nodes, the high parallel efficiency level over 90% can be achieved using a new-type massively parallel computer: ADENART with up to 256 processors. The real performance of the solution methods is superior to those calculated by the vectorized version of the ICCG and ILUBiCGSTAB methods using a vector-type supercomputer.>