From: Ian Jackson <ian@davenant.relativity.greenend.org.uk>
Date: Sun, 30 Dec 2007 23:07:52 +0000 (+0000)
Subject: working on rotation - this is not the decomposition I wanted
X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=commitdiff_plain;h=76f84179ceee1307c87a665a924abaeae8cb9f18

working on rotation - this is not the decomposition I wanted
---

diff --git a/.bzrignore b/.bzrignore
index edc617c..57e94d5 100644
--- a/.bzrignore
+++ b/.bzrignore
@@ -2,3 +2,4 @@ minimise
 primer
 initial
 project
+*.d
diff --git a/Makefile b/Makefile
index cf0bbd0..cd1aa6c 100644
--- a/Makefile
+++ b/Makefile
@@ -4,8 +4,8 @@ TARGETS= minimise primer initial project
 CWARNS=	-Wall -Wwrite-strings -Wpointer-arith -Werror
 
 OPTIMISE=	-O2
-CFLAGS=		$(CWARNS) $(OPTIMISE) -g
-CXXFLAGS=	$(CWARNS) $(OPTIMISE) -g
+CFLAGS=		$(CWARNS) -MMD $(OPTIMISE) -g
+CXXFLAGS=	$(CWARNS) -MMD $(OPTIMISE) -g
 
 LIBGSL= -lgsl -lgslcblas
 
@@ -26,5 +26,6 @@ initial:	generator primer sgtatham/z.typescript
 clean:
 		rm -f *.o $(TARGETS) *.new *.tmp
 		rm -f best initial
+		rm -f *.d
 
-%.o::		common.h mgraph.h bgl.h
+-include *.d
diff --git a/common.h b/common.h
index 52f2e5b..88ec865 100644
--- a/common.h
+++ b/common.h
@@ -46,8 +46,8 @@ void gsldie(int l, const char *what, int status);
 
 #define K FOR_COORD(k)
 
-#define STATIC_GSL_VECTOR(x) static gsl_vector x##_gsl= { D3,1,&x[0] };
-#define STATIC_GSL_MATRIX(x) static gsl_matrix x##_gsl= { D3,D3,D3,&x[0][0] };
+#define GSL_VECTOR(x) gsl_vector x##_gsl= { D3,1,&x[0] };
+#define GSL_MATRIX(x) gsl_matrix x##_gsl= { D3,D3,D3,&x[0][0] };
 
 #ifdef FP_FAST_FMA
 # define fma_fast fma
diff --git a/project.c b/project.c
index b68e473..0c4022e 100644
--- a/project.c
+++ b/project.c
@@ -16,7 +16,7 @@ static int ntris;
 static Vertices conformation;
 
 static double transform[D3][D3]= {{1,0,0}, {0,1,0}, {0,0,1}};
-STATIC_GSL_MATRIX(transform);
+static GSL_MATRIX(transform);
 
 const char *input_filename;
 
@@ -31,16 +31,18 @@ static void read_input(void) {
 }
 
 static void transform_coordinates(void) {
-  /*
-  static double result[D3];
-  static const gsl_vector result_gsl= { D3,1,&result[0]; }
+  double result[D3];
+  GSL_VECTOR(result);
+  gsl_vector input_gsl= { D3,1 };
 
-  int v;
+  int v, k;
   
-
   FOR_VERTEX(v) {
-    GA( gsl_blas_dgemv(CblasNoTrans,
-  }*/
+    input_gsl.data= &conformation[v][0];
+    GA( gsl_blas_dgemv(CblasNoTrans, 1.0,&transform_gsl, &input_gsl,
+		       0.0, &result_gsl) );
+    K conformation[v][k]= result[k];
+  }
 }
 
 static void addtriangle(int va, int vb, int vc) {
@@ -223,6 +225,18 @@ static void drag_none_conclude(void) { }
 static void drag_none_abandon(void) { }
 DRAG(none);
 
+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="";
+  }
+}
+
+#define PMATRIX(x) pmatrix(#x,x);
+
 static void drag_rotate_delta(double dx, double dy) {
   /* We multiple our transformation matrix by this matrix:
    *     [  1   0   0 ]
@@ -232,26 +246,31 @@ 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);
+  static GSL_MATRIX(rotateby);
 
-  static double qr[D3][D3];
-  STATIC_GSL_MATRIX(qr);
+  double qr[D3][D3];
+  gsl_matrix qr_gsl= { D3,D3,D3,&qr[0][0] };
 
-  static double tau[D3];
-  STATIC_GSL_VECTOR(tau);
+  double tau[D3];
+  GSL_VECTOR(tau);
   
-  int k;
-  
-  K rotateby[k][k]= 1;
-  rotateby[0][2]= dx;
-  rotateby[1][2]= dy;
+  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(rotateby);
+
   GA( gsl_linalg_QR_decomp(&qr_gsl, &tau_gsl) );
   GA( gsl_linalg_QR_unpack(&qr_gsl, &tau_gsl,
 			   &transform_gsl, &rotateby_gsl /*dummy*/) );
 
+  PMATRIX(transform);
+  PMATRIX(rotateby);
+
   printf("drag_rotate_delta...\n");
   show();
 }
@@ -359,9 +378,11 @@ int main(int argc, const char *const *argv) {
     case ConfigureNotify:   event_config(&event.xconfigure);  break;
 
     case MotionNotify:
+#if 0
       motion_x= event.xmotion.x;
       motion_y= event.xmotion.y;
       motion_deferred= 1;
+#endif
       continue;
       
     default: