- &rotateby_gsl,&transform_gsl, 0.0,&qr_gsl) );
-
- pmatrix("input", qr);
-
- GA( gsl_linalg_QR_decomp(&qr_gsl, &tau_gsl) );
-
- pmatrix("mangled", qr);
+ &rotateby_gsl,&transform_gsl,
+ 0.0, &skew_gsl) );
+
+ /* Now we want to normalise skew, the result becomes new transform */
+
+ pmatrix("skew", skew);
+
+ /* We use notation from Wikipedia Polar_decomposition
+ * Wikipedia's W is GSL's U
+ * Wikipedia's Sigma is GSL's S
+ * Wikipedia's V is GSL's V
+ * Wikipedia's U is our desired result
+ * Wikipedia which says if the SVD is A = W Sigma V*
+ * then the polar decomposition is A = U P
+ * where P = V Sigma V*
+ * and U = W V*
+ */
+
+ GA( gsl_linalg_SV_decomp(&skew_gsl, &svd_v_gsl, &sigma_gsl, &tau_gsl) );
+ pmatrix("W", skew);
+ pvector("Sigma",sigma);
+ pmatrix("V", svd_v);