A new numerical method, the probabilistic method of discrete ordinates (PMDO) for solving multigroup transport equations in three-dimensional complex geometry, is presented. The method can be used for reactor core and shielding calculations. Integral equations are adopted for the angular flux in cells of arbitrary form. They are coupled by means of net currents defined at interfaces. The sphere of directions is arbitrarily subdivided into a number of angular diapasons. These diapasons, along with cell volume and pieces of cell surface, produce elementary phase domains, so the basic PMDO equations are the algebraic analogues of piece-wise coupled integral transport equations. They are written for neutron flux and currents integrated over corresponding phase domains. The coefficients of the equations discretely depend on the angular variable and have the meaning of probabilities of uncollided neutrons being transmitted between different phase domains. On the basis of algebraic equations separately obtained for coarse and fine domains, the global-local iterative PMDO scheme has also been developed specifically for calculations in extensive heterogeneous media. Together with the direct PMDO equations, the system of conjugate equations has been constructed for the calculation of neutron importance function related to various nonlinear functionals. Codes based on the method and some numerical applications, including examples related to criti-cality calculations and deep penetration problems, have been briefly discribed.