This actually simplifies the math somewhat because in the stress calculation you can just use the (stretch tensor)^2 which is available as an argument in the subroutine definition. If you do it this way, there is no need to rotate the stresses at all...they will automatically be expressed in the reference configuration instead of the deformed configuration.
You actually can use the regular Cauchy stress equation and then rotate it back to the reference configuration with the rotation matrix:
stress_corotational = R^T(stress)R
where R is the rotation tensor and T is transpose operator. You can calculate R and R^T from the deformation gradient and stretch tensor passed in by the subroutine arguments:
F = RU
R = FU^-1
Shawn Chester took this more complicated approach in his Neo-Hookean subroutine at: