We present a theoretical approach to the ab initio calculation of single or multiphoton double electron ionization cross sections, sigma(++)(E), of polyelectronic atoms, near (Wannier region) and far from threshold. The overall computational method is variational, uses functions of real as well as of complex coordinates, and follows the many-electron, many-photon theory proposed by Mercouris and Nicolaides [J. Phys. B 21, L285 (1988); 23, 2037 (1990)]. It incorporates the electronic structure and the pair correlations in the continuum via configuration-interaction techniques. sigma(++)(E) is obtained as the imaginary part of a complex eigenvalue that is computed by diagonalizing a state-specific non-Hermitian matrix constructed from separately optimized function spaces Q and P representing the field-induced resonance state. Q contains correlated wave functions of bound or quasibound states expanded over numerical and analytic orbitals of real coordinates. P is composed, in principle, of subspaces P-1 and P-2, representing the one- and the two-electron channels, respectively, which are optimized separately and then are allowed to mix via the construction of the total non-Hermitian matrix. Both are spanned by basis sets of real coordinates for the ionized core and of complex coordinates for the outgoing part of the one- and the two-electron resonance state. The two-electron square integrable ''continuum'' function space is made orthogonal to the available single electron channels in order for sigma(++)(E) not to include portions of the single electron ionization cross section sigma(+)(E). Application is made to the single photon sigma(++)(E) of the prototypical systems H- and He, but without the mixing of P-2 and P-1, due to numerical instabilities. The two-electron ionization channels were composed of Slater-type orbitals, symmetry-coupled according to (sp), (pd), and (df). Higher symmetries would also be needed at higher energies, with corresponding increase of angular correlation terms in the initial-state wave function. The continuous energy ranged from E=0 to E=250 eV. In the threshold region E=0-2 eV, the length and velocity results are in good agreement with experiment for H- and in reasonable agreement with experiment for He. Far from threshold, there is discrepancy between length and velocity forms in this as well as in previous works by other methods. Apart from whatever inadequacies of the basis functions, this is possibly due to the exclusion of mixing of the single electron open channels into the correlated wave function of the two free electrons. By comparing the results from the use of correlated wave functions with those obtained when the calculation of the transition matrix element is done with wave functions of real coordinates, where the initial state is correlated but the final one is only a product of Coulomb wave functions, the effect of correlation of the two free electrons is deduced for the case of He, without considering the mixing of one- and two-electron channels. Finally, a by-product of the present development was the calculation of the He sigma(+)(E) to the n=1 single ionization threshold. Comparison with previous accurate experimental results reveals very good agreement.