We revise formal and numerical aspects of collinear and noncollinear density functional theory (DFT) in the context of a two-component self-consistent treatment of spin-orbit coupling (SOC). While the extension of the standard one-component theory to a noncollinear magnetization is formally well-defined within the local density approximation, and therefore results in a numerically stable theory, this is not the case within the generalized gradient approximation (GGA). Previously reported formulations of noncollinear DFT based on GGA exchange-correlation potentials have several limitations: (i) they fail at reducing (either formally or numerically) to the proper collinear limit (i.e., when the magnetization is parallel or antiparallel to the z axis everywhere in space); (ii) they fail at ensuring a quantitative rotational invariance of the total energy and even a qualitative rotational invariance of the spatial distribution of the magnetization when a SOC operator is included in the Hamiltonian; (iii) they are numerically very unstable in regions of small magnetization. All of the above-mentioned problems are here shown (both formally and through test examples) to be solved by using instead a new formulation of noncollinear DFT for GGA functionals, which we call the "signed canonical" theory, as combined with an effective screening algorithm for unstable terms of the exchange-correlation potential in regions of small magnetization. All methods are implemented in the CRYSTAL program and tests are performed on simple molecules to compare the different formulations of noncollinear DFT.