The variational formulation and computational aspects of a three-dimensional finite-strain rod model, considered in Part I, are presented. A particular parametrization is employed that bypasses the singularity typically associated with the use of Euler angles. As in the classical Kirchhoff-Love model, rotations have the standard interpretation of orthogonal, generally noncommutative, transformations. This is in contrast with alternative formulations proposed by Argyris et al. [5–8], based on the notion of semitangential rotation. Emphasis is placed on a geometric approach, which proves essential in the formulation of algorithms. In particular, the configuration update procedure becomes the algorithmic counterpart of the exponential map. The computational implementation relies on the formula for the exponential of a skew-symmetric matrix. Consistent linearization procedures are employed to obtain linearized weak forms of the balance equations. The geometric stiffness then becomes generally nonsymmetric as a result of the non-Euclidean character of the configuration space. However, complete symmetry is recovered at an equilibrium configuration, provided that the loading is conservative. An explicit condition for this to be the case is obtained. Numerical simulations including postbuckling behavior and nonconservative loading are also presented. Details pertaining to the implementation of the present formulation are also discussed.