A finite element method for the solution of a general nonlinear, nonhomogeneous and anisotropic magnetostatic problem without hysteresis is considered. A weak formulation involving an integral operator on the boundary only is derived, for which the resulting finite element matrix consists of a sparse part corresponding to the interior of the domain and a full part corresponding to the boundary. Existence and uniqueness of the solution of the weak formulation is proved, and error estimates are derived for the finite element space formed by solenoidal piecewise polynomial vector functions. The solution of the problem by successive iteration is analyzed.