From cd8d998ff14921022cecbfc2475365783f5a8090 Mon Sep 17 00:00:00 2001
From: Ian Jackson <ian@davenant.relativity.greenend.org.uk>
Date: Mon, 31 Dec 2007 00:15:25 +0000
Subject: [PATCH] before check svd

---
 project.c | 70 +++++++++++++++++++++++++++++++++++++------------------
 1 file changed, 47 insertions(+), 23 deletions(-)

diff --git a/project.c b/project.c
index 1fb8576..764a4dd 100644
--- a/project.c
+++ b/project.c
@@ -225,16 +225,21 @@ static void drag_none_conclude(void) { }
 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) {
@@ -248,31 +253,50 @@ 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();
-- 
2.30.2