Merge pull request #38 from McMasterU/matrix_multiplication

Matrix multiplication

Author: dandoh

algorithms/ipopt/ipopt.c

@@ -20,7 +20,9 @@ extern const int var_size[NUM_VARIABLES];
 extern const int var_offset[NUM_VARIABLES];
 extern const int partial_derivative_offset[NUM_VARIABLES];
 extern const int objective_offset;
-extern double ptr[MEM_SIZE];
+
+extern double ptr[MEMORY_NUM_DOUBLES];
+extern complex double ptr_c[MEMORY_NUM_COMPLEX_DOUBLES];
 
 extern const int bound_pos[NUM_VARIABLES];
 extern double lower_bound[NUM_ACTUAL_VARIABLES];

algorithms/lbfgs-b/lbfgs-b.c

@@ -6,6 +6,7 @@
 #include <stdlib.h>
 #include <stdbool.h>
 #include <string.h>
+#include <complex.h>
 
 #define M 10
 
@@ -18,7 +19,9 @@ extern const int var_size[NUM_VARIABLES];
 extern const int var_offset[NUM_VARIABLES];
 extern const int partial_derivative_offset[NUM_VARIABLES];
 extern const int objective_offset;
-extern double ptr[MEM_SIZE];
+
+extern double ptr[MEMORY_NUM_DOUBLES];
+extern complex double ptr_c[MEMORY_NUM_COMPLEX_DOUBLES];
 
 extern void assign_values();
 extern void evaluate_partial_derivatives_and_objective();

algorithms/lbfgs/lbfgs.c

@@ -14,7 +14,9 @@ extern const int var_size[NUM_VARIABLES];
 extern const int var_offset[NUM_VARIABLES];
 extern const int partial_derivative_offset[NUM_VARIABLES];
 extern const int objective_offset;
-extern double ptr[MEM_SIZE];
+
+extern double ptr[MEMORY_NUM_DOUBLES];
+extern complex double ptr_c[MEMORY_NUM_COMPLEX_DOUBLES];
 
 extern void assign_values();
 extern void evaluate_partial_derivatives_and_objective();

embed/fftw.c

@@ -6,72 +6,31 @@
 
 #include <fftw3.h>
 
-// For simplicity, just use the dft for complex by adding 0 to the input complex part
-// We will use the special version for computing dft real later
-void dft_1d(int N, double *in, double *out, int forward_or_backward) {
-  fftw_complex *aux;
-  fftw_plan p;
-  int i;
-  aux = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
-  p = fftw_plan_dft_1d(N, aux, aux, forward_or_backward, FFTW_MEASURE);
-  for (i = 0; i < N; i++) {
-    aux[i][0] = in[i];
-    aux[i][1] = in[i + N];
-  }
-  fftw_execute(p); /* repeat as needed */
-
-
-  if (forward_or_backward == FFTW_FORWARD) {
-    for (i = 0; i < N; i++) {
-      out[i] = aux[i][0];
-      out[i + N] = aux[i][1];
-    }
-  } else {
+void dft_1d(int N, double complex *in, double complex *out, int forward_or_backward) {
+  fftw_plan p = fftw_plan_dft_1d(N, in, out, forward_or_backward, FFTW_ESTIMATE);
+  fftw_execute(p);
+  if (forward_or_backward == FFTW_BACKWARD) {
+    int i;
     for (i = 0; i < N; i++) {
-      out[i] = aux[i][0] / N;
-      out[i + N] = aux[i][1] / N;
+      out[i] = out[i] / N;
     }
   }
 
   fftw_destroy_plan(p);
-  fftw_free(aux);
 }
 
-// For simplicity, just use the dft for complex by adding 0 to the input complex part
-// We will use the special version for computing dft real later
-void dft_2d(int ROW, int COLUMN, double *in, double *out, int forward_or_backward) {
-  fftw_complex *aux;
-  fftw_plan p;
-  int i, j;
+void dft_2d(int ROW, int COLUMN, double complex *in, double complex *out, int forward_or_backward) {
   int size = ROW * COLUMN;
-  aux = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * size);
-  p = fftw_plan_dft_2d(ROW, COLUMN, aux, aux, forward_or_backward, FFTW_MEASURE);
-  for (i = 0; i < ROW; i++) {
-    for (j = 0; j < COLUMN; j++) {
-      aux[i * COLUMN + j][0] = in[i * COLUMN + j];
-      aux[i * COLUMN + j][1] = in[i * COLUMN + j + size];
-    }
-  }
-  fftw_execute(p); /* repeat as needed */
-
-  if (forward_or_backward == FFTW_FORWARD) {
-    for (i = 0; i < ROW; i++) {
-      for (j = 0; j < COLUMN; j++) {
-        out[i * COLUMN + j] = aux[i * COLUMN + j][0];
-        out[i * COLUMN + j + size] = aux[i * COLUMN + j][1];
-      }
-    }
-  } else {
-    for (i = 0; i < ROW; i++) {
-      for (j = 0; j < COLUMN; j++) {
-        out[i * COLUMN + j] = aux[i * COLUMN + j][0] / size;
-        out[i * COLUMN + j + size] = aux[i * COLUMN + j][1] / size;
-      }
+  fftw_plan p = fftw_plan_dft_2d(ROW, COLUMN, in, out, forward_or_backward, FFTW_ESTIMATE);
+  fftw_execute(p);
+  if (forward_or_backward == FFTW_BACKWARD) {
+    int i;
+    for (i = 0; i < size; i++) {
+        out[i] = out[i] / size;
     }
   }
 
   fftw_destroy_plan(p);
-  fftw_free(aux);
 }
 
 /**