In this paper, the element incidence matrix has been extended to develop a comprehensive three-phase distribution system power flow program for radial topology. Three-phase overhead or underground primary feeders and double-phase or single-phase line sections near the end of the feeder laterals have been considered. Unbalanced loads with different types including constant power, constant current and constant impedance are modeled at the system buses. Substation voltage regulator (SVR) consisting of three single phase units connected in wye or two single-phase units connected in open delta are modeled to satisfy the desired voltage level along the feeder. The mathematical model of distributed generation (DG) connected as PQ and PV buses are integrated into the power flow program to simulate the penetration of DGs in the distribution systems. The proposed method has been tested and compared with different IEEE test feeders result. The developed algorithm has been used to study the impact of both SVR and high penetration of DG on voltage profile and system power losses.