annotate src/fftw-3.3.3/support/addchain.c @ 10:37bf6b4a2645

Add FFTW3
author Chris Cannam
date Wed, 20 Mar 2013 15:35:50 +0000
parents
children
rev   line source
Chris@10 1 /* addition-chain optimizer */
Chris@10 2 #include <stdio.h>
Chris@10 3 #include <stdlib.h>
Chris@10 4 #include <unistd.h>
Chris@10 5
Chris@10 6 static int verbose;
Chris@10 7 static int mulcost = 18;
Chris@10 8 static int ldcost = 2;
Chris@10 9 static int sqcost = 10;
Chris@10 10 static int reflcost = 8;
Chris@10 11 #define INFTY 100000
Chris@10 12
Chris@10 13 static int *answer;
Chris@10 14 static int best_so_far;
Chris@10 15
Chris@10 16 static void print_answer(int n, int t)
Chris@10 17 {
Chris@10 18 int i;
Chris@10 19 printf("| (%d, %d) -> [", n, t);
Chris@10 20 for (i = 0; i < t; ++i)
Chris@10 21 printf("%d;", answer[i]);
Chris@10 22 printf("] (* %d *)\n", best_so_far);
Chris@10 23 }
Chris@10 24
Chris@10 25 #define DO(i, j, k, cst) \
Chris@10 26 if (k < n) { \
Chris@10 27 int c = A[i] + A[j] + cst; \
Chris@10 28 if (c < A[k]) { \
Chris@10 29 A[k] = c; \
Chris@10 30 changed = 1; \
Chris@10 31 } \
Chris@10 32 }
Chris@10 33
Chris@10 34 #define DO3(i, j, l, k, cst) \
Chris@10 35 if (k < n) { \
Chris@10 36 int c = A[i] + A[j] + A[l] + cst; \
Chris@10 37 if (c < A[k]) { \
Chris@10 38 A[k] = c; \
Chris@10 39 changed = 1; \
Chris@10 40 } \
Chris@10 41 }
Chris@10 42
Chris@10 43 static int optimize(int n, int *A)
Chris@10 44 {
Chris@10 45 int i, j, k, changed, cst, cstmax;
Chris@10 46
Chris@10 47 do {
Chris@10 48 changed = 0;
Chris@10 49 for (i = 0; i < n; ++i) {
Chris@10 50 k = i + i;
Chris@10 51 DO(i, i, k, sqcost);
Chris@10 52 }
Chris@10 53
Chris@10 54 for (i = 0; i < n; ++i) {
Chris@10 55 for (j = 0; j <= i; ++j) {
Chris@10 56 k = i + j;
Chris@10 57 DO(i, j, k, mulcost);
Chris@10 58 k = i - j;
Chris@10 59 DO(i, j, k, mulcost);
Chris@10 60
Chris@10 61 k = i + j;
Chris@10 62 DO3(i, j, i - j, k, reflcost);
Chris@10 63 }
Chris@10 64 }
Chris@10 65
Chris@10 66 } while (changed);
Chris@10 67
Chris@10 68 cst = cstmax = 0;
Chris@10 69 for (i = 0; i < n; ++i) {
Chris@10 70 cst += A[i];
Chris@10 71 if (A[i] > cstmax) cstmax = A[i];
Chris@10 72 }
Chris@10 73 /* return cstmax; */
Chris@10 74 return cst;
Chris@10 75 }
Chris@10 76
Chris@10 77 static void search(int n, int t, int *A, int *B, int depth)
Chris@10 78 {
Chris@10 79 if (depth == 0) {
Chris@10 80 int i, tc;
Chris@10 81 for (i = 0; i < n; ++i)
Chris@10 82 A[i] = INFTY;
Chris@10 83 A[0] = 0; /* always free */
Chris@10 84 for (i = 1; i <= t; ++i)
Chris@10 85 A[B[-i]] = ldcost;
Chris@10 86
Chris@10 87 tc = optimize(n, A);
Chris@10 88 if (tc < best_so_far) {
Chris@10 89 best_so_far = tc;
Chris@10 90 for (i = 1; i <= t; ++i)
Chris@10 91 answer[t - i] = B[-i];
Chris@10 92 if (verbose)
Chris@10 93 print_answer(n, t);
Chris@10 94 }
Chris@10 95 } else {
Chris@10 96 for (B[0] = B[-1] + 1; B[0] < n; ++B[0])
Chris@10 97 search(n, t, A, B + 1, depth - 1);
Chris@10 98 }
Chris@10 99 }
Chris@10 100
Chris@10 101 static void doit(int n, int t)
Chris@10 102 {
Chris@10 103 int *A;
Chris@10 104 int *B;
Chris@10 105
Chris@10 106 A = malloc(n * sizeof(int));
Chris@10 107 B = malloc((t + 1) * sizeof(int));
Chris@10 108 answer = malloc(t * sizeof(int));
Chris@10 109
Chris@10 110 B[0] = 0;
Chris@10 111 best_so_far = INFTY;
Chris@10 112 search(n, t, A, B + 1, t);
Chris@10 113
Chris@10 114 print_answer(n, t);
Chris@10 115
Chris@10 116 free(A); free(B); free(answer);
Chris@10 117 }
Chris@10 118
Chris@10 119 int main(int argc, char *argv[])
Chris@10 120 {
Chris@10 121 int n = 32;
Chris@10 122 int t = 3;
Chris@10 123 int all;
Chris@10 124 int ch;
Chris@10 125
Chris@10 126 verbose = 0;
Chris@10 127 all = 0;
Chris@10 128 while ((ch = getopt(argc, argv, "n:t:m:l:r:s:va")) != -1) {
Chris@10 129 switch (ch) {
Chris@10 130 case 'n':
Chris@10 131 n = atoi(optarg);
Chris@10 132 break;
Chris@10 133 case 't':
Chris@10 134 t = atoi(optarg);
Chris@10 135 break;
Chris@10 136 case 'm':
Chris@10 137 mulcost = atoi(optarg);
Chris@10 138 break;
Chris@10 139 case 'l':
Chris@10 140 ldcost = atoi(optarg);
Chris@10 141 break;
Chris@10 142 case 's':
Chris@10 143 sqcost = atoi(optarg);
Chris@10 144 break;
Chris@10 145 case 'r':
Chris@10 146 reflcost = atoi(optarg);
Chris@10 147 break;
Chris@10 148 case 'v':
Chris@10 149 ++verbose;
Chris@10 150 break;
Chris@10 151 case 'a':
Chris@10 152 ++all;
Chris@10 153 break;
Chris@10 154 case '?':
Chris@10 155 fprintf(stderr, "use the source\n");
Chris@10 156 exit(1);
Chris@10 157 }
Chris@10 158 }
Chris@10 159
Chris@10 160 if (all) {
Chris@10 161 for (n = 4; n <= 64; n *= 2) {
Chris@10 162 int n1 = n - 1; if (n1 > 7) n1 = 7;
Chris@10 163 for (t = 1; t <= n1; ++t)
Chris@10 164 doit(n, t);
Chris@10 165 }
Chris@10 166 } else {
Chris@10 167 doit(n, t);
Chris@10 168 }
Chris@10 169
Chris@10 170 return 0;
Chris@10 171 }