A computer-implemented method for simulating operation of an integrated circuit, the method comprising: using a parallel processing system, generating a system of differential or algebraic equations modeling operation of the integrated circuit; discretizing the system of DAE to form a discretized DAE; partitioning the integrated circuit into load balanced partitions, each partition representing a subsystem of the integrated circuit; for each partition representing a subsystem of the integrated circuit: when the DAE are nonlinear, linearizing the discretized DAE using a nonlinear iterative method to obtain a linear system having a circuit Jacobian matrix; solving the linear system using a linearized circuit Jacobian matrix solver, wherein the solving comprises: splitting the circuit Jacobian matrix into two matrices M and N, wherein M is a matrix suitable for parallel processing and N is a coupling matrix; preconditioning the two matrices M and N to form a preconditioned equation having a matrix of form I + M -1 N, where I is an identity matrix having ones on a diagonal and zeros on off-diagonal positions of the identity matrix; and obtaining a solution to I + M -1 N in the preconditioned equation using a combined direct and iterative solving method; and recursively performing the steps from partitioning the circuit till obtaining a solution to I+M -1 N to reduce sizes of subsystems of the integrated circuit; wherein the generating the system of DAE comprises mapping a design of the integrated circuit to the system of DAE, the system of DAE defining a change in voltage at each of a plurality of nodes in the integrated circuit, the change in voltage occurring over a definable period of time. using a parallel processing system, generating a system of differential or algebraic equations modeling operation of the integrated circuit; discretizing the system of DAE to form a discretized DAE; partitioning the integrated circuit into load balanced partitions, each partition representing a subsystem of the integrated circuit; for each partition representing a subsystem of the integrated circuit: when the DAE are nonlinear, linearizing the discretized DAE using a nonlinear iterative method to obtain a linear system having a circuit Jacobian matrix; solving the linear system using a linearized circuit Jacobian matrix solver, wherein the solving comprises: splitting the circuit Jacobian matrix into two matrices M and N, wherein M is a matrix suitable for parallel processing and N is a coupling matrix; preconditioning the two matrices M and N to form a preconditioned equation having a matrix of form I + M -1 N, where I is an identity matrix having ones on a diagonal and zeros on off-diagonal positions of the identity matrix; and obtaining a solution to I + M -1 N in the preconditioned equation using a combined direct and iterative solving method; and recursively performing the steps from partitioning the circuit till obtaining a solution to I+M -1 N to reduce sizes of subsystems of the integrated circuit; wherein the generating the system of DAE comprises mapping a design of the integrated circuit to the system of DAE, the system of DAE defining a change in voltage at each of a plurality of nodes in the integrated circuit, the change in voltage occurring over a definable period of time. when the DAE are nonlinear, linearizing the discretized DAE using a nonlinear iterative method to obtain a linear system having a circuit Jacobian matrix; solving the linear system using a linearized circuit Jacobian matrix solver, wherein the solving comprises: splitting the circuit Jacobian matrix into two matrices M and N, wherein M is a matrix suitable for parallel processing and N is a coupling matrix; preconditioning the two matrices M and N to form a preconditioned equation having a matrix of form I + M -1 N, where I is an identity matrix having ones on a diagonal and zeros on off-diagonal positions of the identity matrix; and obtaining a solution to I + M -1 N in the preconditioned equation using a combined direct and iterative solving method; and splitting the circuit Jacobian matrix into two matrices M and N, wherein M is a matrix suitable for parallel processing and N is a coupling matrix; preconditioning the two matrices M and N to form a preconditioned equation having a matrix of form I + M -1 N, where I is an identity matrix having ones on a diagonal and zeros on off-diagonal positions of the identity matrix; and obtaining a solution to I + M -1 N in the preconditioned equation using a combined direct and iterative solving method; and recursively performing the steps from partitioning the circuit till obtaining a solution to I+M -1 N to reduce sizes of subsystems of the integrated circuit; wherein the generating the system of DAE comprises mapping a design of the integrated circuit to the system of DAE, the system of DAE defining a change in voltage at each of a plurality of nodes in the integrated circuit, the change in voltage occurring over a definable period of time.