In this work, a relativistic version of the state-averaged complete active space self-consistent field method is developed (spin-orbit coupled state-averaged complete active space self-consistent field; CAS-SOC). The program follows a "one-step strategy" and treats the spin-orbit interaction (SOC) on the same footing as the electron-electron interaction. As opposed to other existing approaches, the program employs an intermediate coupling scheme in which spin and space symmetry adapted configuration space functions are allowed to interact via SOC. This adds to the transparency and computational efficiency of the procedure. The approach requires the utilization of complex-valued configuration interaction coefficients, but the molecular orbital coefficients can be kept real-valued without loss of generality. Hence, expensive arithmetic associated with evaluation of complex-valued transformed molecular integrals is completely avoided. In order to investigate the quality of the calculated wave function, we extended the method to the calculation of electronic g-tensors. As the SOC is already treated to all orders in the SA-CASSCF process, first order perturbation theory with the Zeeman operator is sufficient to accomplish this task. As a test-set, we calculated g-tensors of a set of diatomics, a set of d(1) transition metal complexes MOX4 (n-), and a set of 5f(1) actinide complexes AnX6 (n-). These calculations reveal that the effect of the wavefunction relaxation due to variation inclusion of SOC is of the same order of magnitude as the effect of inclusion of dynamic correlation and hence cannot be neglected for the accurate prediction of electronic g-tensors.