+ GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
+ &adjr_gsl,&temp_gsl,
+ 0.0, &change_gsl) );
+ PMATRIX(change);
+
+ static double skew[D3][D3];
+ static GSL_MATRIX(skew);
+
+ GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
+ &change_gsl,&transform_gsl,
+ 0.0, &skew_gsl) );
+ PMATRIX(skew);
+
+ memcpy(&transform,&skew,sizeof(transform));
+ show();
+ return;
+
+ /* Now we want to normalise skew, the result becomes new transform */
+ double svd_v[D3][D3];
+ GSL_MATRIX(svd_v);
+
+ double sigma[D3], tau[D3];
+ GSL_VECTOR(sigma);
+ GSL_VECTOR(tau);
+
+ /* 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*
+ */