chiark / gitweb /
working on rotation - this is not the decomposition I wanted
authorIan Jackson <ian@davenant.relativity.greenend.org.uk>
Sun, 30 Dec 2007 23:07:52 +0000 (23:07 +0000)
committerIan Jackson <ian@davenant.relativity.greenend.org.uk>
Sun, 30 Dec 2007 23:07:52 +0000 (23:07 +0000)
.bzrignore
Makefile
common.h
project.c

index edc617c..57e94d5 100644 (file)
@@ -2,3 +2,4 @@ minimise
 primer
 initial
 project
+*.d
index cf0bbd0..cd1aa6c 100644 (file)
--- 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
index 52f2e5b..88ec865 100644 (file)
--- 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
index b68e473..0c4022e 100644 (file)
--- 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: