The polar decomposition of an $m x n$ matrix $A$ of full rank, where $m \geq n$, can be computed using a quadratically convergent algorithm of Higham [SIAM J. Sci. Stat. Comput., 7 (1986), pp.1160-1174]. The algorithm is based on a Newton iteration involving a matrix inverse. We show how with the use of a preliminary complete orthogonal decomposition the algorithm can be extended to arbitrary $A$. We also describe how to use the algorithm to compute the positive semi-definite square root of a Hermitian positive semi-definite matrix. We formulate a hybrid algorithm which adaptively switches from the matrix inversion based iteration to a matrix multiplication based iteration due to Kovarik, and to Bjorck and Bowie. The decision when to switch is made using a condition estimator. This "matrix multiplication rich" algorithm is shown to be more efficient on machines for which matrix multiplication can be executed 1.5 times faster than matrix inversion.