static void drag_none_abandon(void) { }
DRAG(none);
+static void pvectorcore(const char *n, double v[D3]) {
+ int k;
+ printf("%10s [ ",n);
+ K printf("%# 10.10f ",v[k]);
+ printf("]\n");
+}
+static void pvector(const char *n, double v[D3]) {
+ pvectorcore(m,v);
+ putchar('\n');
+}
static void pmatrix(const char *n, double m[D3][D3]) {
- int j,k;
- for (j=0; j<D3; j++) {
- printf("%10s [ ",n);
- K printf("%# f ",m[j][k]);
- printf("]\n");
- n="";
- }
+ int j;
+ for (j=0; j<D3; j++) { pvector(n,m[j]); n=""; }
+ putchar('\n');
}
-
#define PMATRIX(x) pmatrix(#x,x);
static void drag_rotate_delta(double dx, double dy) {
static double rotateby[D3][D3]= {{1,0,0},{0,1,0},{-20,-20,1}};
static GSL_MATRIX(rotateby);
- double qr[D3][D3];
- gsl_matrix qr_gsl= { D3,D3,D3,&qr[0][0] };
+ static double skew[D3][D3];
+ static GSL_MATRIX(skew);
+
+ double svd_v[D3][D3];
+ GSL_MATRIX(svd_v);
- double tau[D3];
+ double sigma[D3], tau[D3];
+ GSL_VECTOR(sigma);
GSL_VECTOR(tau);
rotateby[2][0]= dx;
rotateby[2][1]= dy;
-
PMATRIX(rotateby);
GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
- &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);
- GA( gsl_linalg_QR_unpack(&qr_gsl, &tau_gsl,
- &transform_gsl, &rotateby_gsl /*dummy*/) );
+ /* We only need U, not P. */
+ GA( gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0,
+ &skew_gsl,&svd_v_gsl,
+ 0.0,&transform_gsl) );
- pmatrix("Q", transform);
- pmatrix("R", rotateby);
+ pmatrix("U", transform);
printf("drag_rotate_delta...\n");
show();