luke
/
ngspice
mirror of https://git.code.sf.net/p/ngspice/ngspice
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

fft calculation with type double (instead of float)

This commit is contained in:
h_vogt 2011-08-19 20:46:23 +00:00
parent 56229ce907
commit 220cc83458
7 changed files with 376 additions and 466 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
ChangeLog
View File

@ -1,6 +1,10 @@
2011-08-19 Holger Vogt 2011-08-19 Holger Vogt
* rawfile.c: code beautified, 'strange line' error improved * rawfile.c: code beautified, 'strange line' error improved
* maths/fft/fftext.h removed
* include/fftext.h, 1-f-code.c, 1-f-code.h, com_fft.c:
type float --> double
* control.c: code beautified
2011-08-14 Robert Larice 2011-08-14 Robert Larice
* src/frontend/postsc.c : * src/frontend/postsc.c :
postsc.c, swallow warnings postsc.c, swallow warnings

2
src/frontend/com_fft.c
View File

@ -16,7 +16,7 @@ Author: 2008 Holger Vogt
#include "variable.h" #include "variable.h"
#include "parse.h" #include "parse.h"
#include "../misc/misc_time.h" #include "../misc/misc_time.h"
#include "../maths/fft/fftext.h" #include "fftext.h"
static void fftext(double*, double*, long int, long int, int); static void fftext(double*, double*, long int, long int, int);

636
src/frontend/control.c
View File

@ -26,7 +26,7 @@ Author: 1985 Wayne A. Christopher, U. C. Berkeley CAD Group
/* Are we waiting for a command? This lets signal handling be /* Are we waiting for a command? This lets signal handling be
* more clever. */ * more clever. */
bool cp_cwait = FALSE; bool cp_cwait = FALSE;
char *cp_csep = ";"; char *cp_csep = ";";
bool cp_dounixcom = FALSE; bool cp_dounixcom = FALSE;
@ -47,22 +47,21 @@ int stackp = 0;
* blown away every time we return -- probably every time we type * blown away every time we return -- probably every time we type
* source at the keyboard and every time a source returns to keyboard * source at the keyboard and every time a source returns to keyboard
* input is ok though -- use ft_controlreset. */ * input is ok though -- use ft_controlreset. */
/* Notes by CDHW: /* Notes by CDHW:
* This routine leaked like a sieve because each getcommand() created a * This routine leaked like a sieve because each getcommand() created a
* wordlist that was never freed because it might have been added into * wordlist that was never freed because it might have been added into
* the control structure. I've tackled this by making sure that everything * the control structure. I've tackled this by making sure that everything
* put into the cend[stackp] is a copy. This means that wlist can be * put into the cend[stackp] is a copy. This means that wlist can be
* destroyed safely * destroyed safely
*/ */
/* no redirection after the following commands (we may need more to add here!) */ /* no redirection after the following commands (we may need more to add here!) */
static char *noredirect[] = { "stop", "define", NULL } ; static char *noredirect[] = { "stop", "define", NULL } ;
static struct control * static struct control *
findlabel(char *s, struct control *ct) findlabel(char *s, struct control *ct) {
{
while (ct) { while (ct) {
if ((ct->co_type == CO_LABEL) && eq(s, ct->co_text->wl_word)) if ((ct->co_type == CO_LABEL) && eq(s, ct->co_text->wl_word))
break; break;
@ -106,15 +105,23 @@ pwlist_echo(wordlist *wlist, char *name) /*CDHW used to perform function of se
/*CDHW Remove control structure and free the memory its hogging CDHW*/ /*CDHW Remove control structure and free the memory its hogging CDHW*/
static void ctl_free(struct control *ctrl) { static void ctl_free(struct control *ctrl)
if (!ctrl) return; {
wl_free(ctrl->co_cond); ctrl->co_cond = NULL; if (!ctrl) return;
tfree(ctrl->co_foreachvar); ctrl->co_foreachvar = NULL; wl_free(ctrl->co_cond);
wl_free(ctrl->co_text); ctrl->co_text = NULL; ctrl->co_cond = NULL;
ctl_free(ctrl->co_children); ctrl->co_children = NULL; tfree(ctrl->co_foreachvar);
ctl_free(ctrl->co_elseblock); ctrl->co_elseblock = NULL; ctrl->co_foreachvar = NULL;
ctl_free(ctrl->co_next); ctrl->co_next = NULL; wl_free(ctrl->co_text);
tfree(ctrl); ctrl = NULL; ctrl->co_text = NULL;
ctl_free(ctrl->co_children);
ctrl->co_children = NULL;
ctl_free(ctrl->co_elseblock);
ctrl->co_elseblock = NULL;
ctl_free(ctrl->co_next);
ctrl->co_next = NULL;
tfree(ctrl);
ctrl = NULL;
} }
@ -145,7 +152,7 @@ docommand(wordlist *wlist)
wlist = cp_doglob(wlist); wlist = cp_doglob(wlist);
pwlist(wlist, "After globbing"); pwlist(wlist, "After globbing");
pwlist_echo(wlist, "Becomes >"); pwlist_echo(wlist, "Becomes >");
if (!wlist || !wlist->wl_word) /*CDHW need to free wlist in second case? CDHW*/ if (!wlist || !wlist->wl_word) /*CDHW need to free wlist in second case? CDHW*/
@ -191,13 +198,13 @@ docommand(wordlist *wlist)
for (i = 0; cp_coms[i].co_comname; i++) { for (i = 0; cp_coms[i].co_comname; i++) {
/* strcmp(cp_coms[i].co_comname, s) ... */ /* strcmp(cp_coms[i].co_comname, s) ... */
for (t = cp_coms[i].co_comname, r = s; *t && *r; for (t = cp_coms[i].co_comname, r = s; *t && *r;
t++, r++) t++, r++)
if (*t != *r) if (*t != *r)
break; break;
if (!*t && !*r) if (!*t && !*r)
break; break;
} }
/* Now give the user-supplied command routine a try... */ /* Now give the user-supplied command routine a try... */
if (!cp_coms[i].co_func && cp_oddcomm(s, wlist->wl_next)) if (!cp_coms[i].co_func && cp_oddcomm(s, wlist->wl_next))
goto out; goto out;
@ -207,14 +214,14 @@ docommand(wordlist *wlist)
if (cp_dounixcom && cp_unixcom(wlist)) if (cp_dounixcom && cp_unixcom(wlist))
goto out; goto out;
fprintf(cp_err,"%s: no such command available in %s\n", fprintf(cp_err,"%s: no such command available in %s\n",
s, cp_program); s, cp_program);
goto out; goto out;
/* If it hasn't been implemented */ /* If it hasn't been implemented */
} else if (!cp_coms[i].co_func) { } else if (!cp_coms[i].co_func) {
fprintf(cp_err,"%s: command is not implemented\n", s); fprintf(cp_err,"%s: command is not implemented\n", s);
goto out; goto out;
/* If it's there but spiceonly, and this is nutmeg, error. */ /* If it's there but spiceonly, and this is nutmeg, error. */
} else if (ft_nutmeg && cp_coms[i].co_spiceonly) { } else if (ft_nutmeg && cp_coms[i].co_spiceonly) {
fprintf(cp_err,"%s: command available only in spice\n", s); fprintf(cp_err,"%s: command available only in spice\n", s);
goto out; goto out;
@ -227,11 +234,11 @@ docommand(wordlist *wlist)
nargs++; nargs++;
{ {
if (nargs < command->co_minargs) { if (nargs < command->co_minargs) {
if (command->co_argfn) { if (command->co_argfn) {
command->co_argfn (wlist->wl_next, command); command->co_argfn (wlist->wl_next, command);
} else { } else {
fprintf(cp_err, "%s: too few args.\n", s); fprintf(cp_err, "%s: too few args.\n", s);
} }
} else if (nargs > command->co_maxargs) { } else if (nargs > command->co_maxargs) {
fprintf(cp_err, "%s: too many args.\n", s); fprintf(cp_err, "%s: too many args.\n", s);
} else { } else {
@ -240,7 +247,8 @@ docommand(wordlist *wlist)
} }
/* Now fix the pointers and advance wlist. */ /* Now fix the pointers and advance wlist. */
out: wlist->wl_prev = ee; out:
wlist->wl_prev = ee;
if (nextc) { if (nextc) {
for(wl=wlist; wl->wl_next; wl=wl->wl_next) for(wl=wlist; wl->wl_next; wl=wl->wl_next)
; ;
@ -282,14 +290,14 @@ doblock(struct control *bl, int *num)
wordlist *wl; wordlist *wl;
char *i; char *i;
int nn; int nn;
nn = *num + 1 ; /*CDHW this is a guess... CDHW*/ nn = *num + 1 ; /*CDHW this is a guess... CDHW*/
switch (bl->co_type) { switch (bl->co_type) {
case CO_WHILE: case CO_WHILE:
if (!bl->co_children) { if (!bl->co_children) {
fprintf(cp_err, "Warning: Executing empty 'while' block.\n"); fprintf(cp_err, "Warning: Executing empty 'while' block.\n");
fprintf(cp_err, " (Use a label statement as a no-op to suppress this warning.)\n"); fprintf(cp_err, " (Use a label statement as a no-op to suppress this warning.)\n");
} }
while (bl->co_cond && cp_istrue(bl->co_cond)) { while (bl->co_cond && cp_istrue(bl->co_cond)) {
if (!bl->co_children) cp_periodic(); /*CDHW*/ if (!bl->co_children) cp_periodic(); /*CDHW*/
@ -298,10 +306,10 @@ doblock(struct control *bl, int *num)
i = doblock(ch, &nn); i = doblock(ch, &nn);
switch (*i) { switch (*i) {
case NORMAL: case NORMAL:
break; break;
case BROKEN: /* Break. */ case BROKEN: /* Break. */
if (nn < 2) if (nn < 2)
return (NORMAL_STR); return (NORMAL_STR);
else { else {
@ -309,7 +317,7 @@ doblock(struct control *bl, int *num)
return (BROKEN_STR); return (BROKEN_STR);
} }
case CONTINUED: /* Continue. */ case CONTINUED: /* Continue. */
if (nn < 2) { if (nn < 2) {
cn = NULL; cn = NULL;
break; break;
@ -318,7 +326,7 @@ doblock(struct control *bl, int *num)
return (CONTINUED_STR); return (CONTINUED_STR);
} }
default: default:
cn = findlabel(i, bl->co_children); cn = findlabel(i, bl->co_children);
if (!cn) if (!cn)
return (i); return (i);
@ -334,10 +342,10 @@ doblock(struct control *bl, int *num)
i = doblock(ch, &nn); i = doblock(ch, &nn);
switch (*i) { switch (*i) {
case NORMAL: case NORMAL:
break; break;
case BROKEN: /* Break. */ case BROKEN: /* Break. */
if (nn < 2) if (nn < 2)
return (NORMAL_STR); return (NORMAL_STR);
else { else {
@ -345,7 +353,7 @@ doblock(struct control *bl, int *num)
return (BROKEN_STR); return (BROKEN_STR);
} }
case CONTINUED: /* Continue. */ case CONTINUED: /* Continue. */
if (nn < 2) { if (nn < 2) {
cn = NULL; cn = NULL;
break; break;
@ -354,7 +362,7 @@ doblock(struct control *bl, int *num)
return (CONTINUED_STR); return (CONTINUED_STR);
} }
default: default:
cn = findlabel(i, bl->co_children); cn = findlabel(i, bl->co_children);
if (!cn) if (!cn)
return (i); return (i);
@ -365,15 +373,15 @@ doblock(struct control *bl, int *num)
case CO_REPEAT: case CO_REPEAT:
if (!bl->co_children) { if (!bl->co_children) {
fprintf(cp_err, "Warning: Executing empty 'repeat' block.\n"); fprintf(cp_err, "Warning: Executing empty 'repeat' block.\n");
fprintf(cp_err, " (Use a label statement as a no-op to suppress this warning.)\n"); fprintf(cp_err, " (Use a label statement as a no-op to suppress this warning.)\n");
} }
if (!bl->co_timestodo) bl->co_timestodo = bl->co_numtimes; if (!bl->co_timestodo) bl->co_timestodo = bl->co_numtimes;
/*bl->co_numtimes: total repeat count /*bl->co_numtimes: total repeat count
bl->co_numtimes = -1: repeat forever bl->co_numtimes = -1: repeat forever
bl->co_timestodo: remaining repeats*/ bl->co_timestodo: remaining repeats*/
while ((bl->co_timestodo > 0) || while ((bl->co_timestodo > 0) ||
(bl->co_timestodo == -1)) { (bl->co_timestodo == -1)) {
if (!bl->co_children) cp_periodic(); /*CDHW*/ if (!bl->co_children) cp_periodic(); /*CDHW*/
if (bl->co_timestodo != -1) bl->co_timestodo--; if (bl->co_timestodo != -1) bl->co_timestodo--;
/* loop through all stements inside rpeat ... end */ /* loop through all stements inside rpeat ... end */
@ -382,10 +390,10 @@ doblock(struct control *bl, int *num)
i = doblock(ch, &nn); i = doblock(ch, &nn);
switch (*i) { switch (*i) {
case NORMAL: case NORMAL:
break; break;
case BROKEN: /* Break. */ case BROKEN: /* Break. */
/* before leaving repeat loop set remaining timestodo to 0 */ /* before leaving repeat loop set remaining timestodo to 0 */
bl->co_timestodo = 0; bl->co_timestodo = 0;
if (nn < 2) if (nn < 2)
@ -395,7 +403,7 @@ doblock(struct control *bl, int *num)
return (BROKEN_STR); return (BROKEN_STR);
} }
case CONTINUED: /* Continue. */ case CONTINUED: /* Continue. */
if (nn < 2) { if (nn < 2) {
cn = NULL; cn = NULL;
break; break;
@ -406,12 +414,12 @@ doblock(struct control *bl, int *num)
return (CONTINUED_STR); return (CONTINUED_STR);
} }
default: default:
cn = findlabel(i, bl->co_children); cn = findlabel(i, bl->co_children);
if (!cn) { if (!cn) {
/* no label found inside repeat loop: /* no label found inside repeat loop:
before leaving loop set remaining timestodo to 0 */ before leaving loop set remaining timestodo to 0 */
bl->co_timestodo = 0; bl->co_timestodo = 0;
return (i); return (i);
} }
@ -427,7 +435,7 @@ doblock(struct control *bl, int *num)
i = doblock(ch, &nn); i = doblock(ch, &nn);
if (*i > 2) { if (*i > 2) {
cn = findlabel(i, cn = findlabel(i,
bl->co_children); bl->co_children);
if (!cn) if (!cn)
return (i); return (i);
} else if (*i != NORMAL) { } else if (*i != NORMAL) {
@ -441,7 +449,7 @@ doblock(struct control *bl, int *num)
i = doblock(ch, &nn); i = doblock(ch, &nn);
if (*i > 2) { if (*i > 2) {
cn = findlabel(i, cn = findlabel(i,
bl->co_elseblock); bl->co_elseblock);
if (!cn) if (!cn)
return (i); return (i);
} else if (*i != NORMAL) { } else if (*i != NORMAL) {
@ -454,18 +462,18 @@ doblock(struct control *bl, int *num)
case CO_FOREACH: case CO_FOREACH:
for (wl = cp_variablesubst(cp_bquote(cp_doglob(wl_copy(bl->co_text)))); for (wl = cp_variablesubst(cp_bquote(cp_doglob(wl_copy(bl->co_text))));
wl; wl;
wl = wl->wl_next) { wl = wl->wl_next) {
cp_vset(bl->co_foreachvar, CP_STRING, wl->wl_word); cp_vset(bl->co_foreachvar, CP_STRING, wl->wl_word);
for (ch = bl->co_children; ch; ch = cn) { for (ch = bl->co_children; ch; ch = cn) {
cn = ch->co_next; cn = ch->co_next;
i = doblock(ch, &nn); i = doblock(ch, &nn);
switch (*i) { switch (*i) {
case NORMAL: case NORMAL:
break; break;
case BROKEN: /* Break. */ case BROKEN: /* Break. */
if (nn < 2) if (nn < 2)
return (NORMAL_STR); return (NORMAL_STR);
else { else {
@ -473,7 +481,7 @@ doblock(struct control *bl, int *num)
return (BROKEN_STR); return (BROKEN_STR);
} }
case CONTINUED: /* Continue. */ case CONTINUED: /* Continue. */
if (nn < 2) { if (nn < 2) {
cn = NULL; cn = NULL;
break; break;
@ -482,7 +490,7 @@ doblock(struct control *bl, int *num)
return (CONTINUED_STR); return (CONTINUED_STR);
} }
default: default:
cn = findlabel(i, bl->co_children); cn = findlabel(i, bl->co_children);
if (!cn) if (!cn)
return (i); return (i);
@ -513,13 +521,13 @@ doblock(struct control *bl, int *num)
case CO_GOTO: case CO_GOTO:
wl = cp_variablesubst(cp_bquote(cp_doglob( wl = cp_variablesubst(cp_bquote(cp_doglob(
wl_copy(bl->co_text)))); /*CDHW Leak ? CDHW*/ wl_copy(bl->co_text)))); /*CDHW Leak ? CDHW*/
return (wl->wl_word); return (wl->wl_word);
case CO_LABEL: case CO_LABEL:
/* Do nothing. */ /* Do nothing. */
cp_periodic(); /*CDHW needed to avoid lock-ups when loop contains only a label CDHW*/ cp_periodic(); /*CDHW needed to avoid lock-ups when loop contains only a label CDHW*/
break; break;
case CO_STATEMENT: case CO_STATEMENT:
docommand(wl_copy(bl->co_text)); docommand(wl_copy(bl->co_text));
@ -531,9 +539,9 @@ doblock(struct control *bl, int *num)
break; break;
default: default:
fprintf(cp_err, fprintf(cp_err,
"doblock: Internal Error: bad block type %d\n", "doblock: Internal Error: bad block type %d\n",
bl->co_type); bl->co_type);
return (NORMAL_STR); return (NORMAL_STR);
} }
return (NORMAL_STR); return (NORMAL_STR);
@ -546,41 +554,41 @@ doblock(struct control *bl, int *num)
Number of chevrons indicates stack depth. Number of chevrons indicates stack depth.
Returns NULL when there is no alternate prompt. Returns NULL when there is no alternate prompt.
SJB 28th April 2005 */ SJB 28th April 2005 */
char * char *
get_alt_prompt(void) get_alt_prompt(void)
{ {
int i = 0, j; int i = 0, j;
static char buf[MAX_CHEVRONS + 2]; /* includes terminating space & null */ static char buf[MAX_CHEVRONS + 2]; /* includes terminating space & null */
struct control *c; struct control *c;
/* if nothing on the command stack return NULL */ /* if nothing on the command stack return NULL */
if (cend[stackp] == 0) if (cend[stackp] == 0)
return NULL; return NULL;
/* measure stack depth */ /* measure stack depth */
for (c = cend[stackp]->co_parent; c; c = c->co_parent) for (c = cend[stackp]->co_parent; c; c = c->co_parent)
i++; i++;
if (i <= 0) if (i <= 0)
return NULL; return NULL;
/* Avoid overflow of buffer and /* Avoid overflow of buffer and
indicate when we've limited the chevrons by starting with a '+' */ indicate when we've limited the chevrons by starting with a '+' */
if(i > MAX_CHEVRONS) { if(i > MAX_CHEVRONS) {
i = MAX_CHEVRONS; i = MAX_CHEVRONS;
buf[0]='+'; buf[0]='+';
} else { } else {
buf[0]='>'; buf[0]='>';
} }
/* return one chevron per command stack depth */ /* return one chevron per command stack depth */
for (j = 1; j < i; j++) for (j = 1; j < i; j++)
buf[j] = '>'; buf[j] = '>';
/* Add space and terminate */ /* Add space and terminate */
buf[j] = ' '; buf[j] = ' ';
buf[j + 1] = '\0'; buf[j + 1] = '\0';
return buf; return buf;
} }
@ -599,7 +607,7 @@ getcommand(char *string)
/* set cp_altprompt for use by the lexer - see parser/lexical.c */ /* set cp_altprompt for use by the lexer - see parser/lexical.c */
cp_altprompt = get_alt_prompt(); cp_altprompt = get_alt_prompt();
#endif /* !defined(HAVE_GNUREADLINE) && !defined(HAVE_BSDEDITLINE) */ #endif /* !defined(HAVE_GNUREADLINE) && !defined(HAVE_BSDEDITLINE) */
cp_cwait = TRUE; cp_cwait = TRUE;
wlist = cp_parse(string); wlist = cp_parse(string);
cp_cwait = FALSE; cp_cwait = FALSE;
@ -627,232 +635,238 @@ cp_evloop(char *string)
cend[stackp]->co_type = CO_UNFILLED; cend[stackp]->co_type = CO_UNFILLED;
for (;;) { for (;;) {
wlist = getcommand(string); wlist = getcommand(string);
if (wlist == NULL) { /* End of file or end of user input. */ if (wlist == NULL) { /* End of file or end of user input. */
if (cend[stackp]->co_parent && !string) { if (cend[stackp]->co_parent && !string) {
cp_resetcontrol(); cp_resetcontrol();
continue; continue;
} else } else
return (0); return (0);
} }
if ((wlist->wl_word == NULL) || (*wlist->wl_word == '\0')) { if ((wlist->wl_word == NULL) || (*wlist->wl_word == '\0')) {
/* User just typed return. */ /* User just typed return. */
wl_free(wlist); /* va, avoid memory leak */ wl_free(wlist); /* va, avoid memory leak */
if (string) if (string)
return (1); return (1);
else { else {
cp_event--; cp_event--;
continue; continue;
} }
} }
/* Just a check... */ /* Just a check... */
for (ww = wlist; ww; ww = ww->wl_next) for (ww = wlist; ww; ww = ww->wl_next)
if (!ww->wl_word) { if (!ww->wl_word) {
fprintf(cp_err, fprintf(cp_err,
"cp_evloop: Internal Error: NULL word pointer\n"); "cp_evloop: Internal Error: NULL word pointer\n");
continue; continue;
} }
/* Add this to the control structure list. If cend->co_type is /* Add this to the control structure list. If cend->co_type is
* CO_UNFILLED, the last line was the beginning of a block, * CO_UNFILLED, the last line was the beginning of a block,
* and this is the unfilled first statement. * and this is the unfilled first statement.
*/ */
/* va: TODO: free old structure and its content, before overwriting */ /* va: TODO: free old structure and its content, before overwriting */
if (cend[stackp] && (cend[stackp]->co_type != CO_UNFILLED)) { if (cend[stackp] && (cend[stackp]->co_type != CO_UNFILLED)) {
cend[stackp]->co_next = alloc(struct control); cend[stackp]->co_next = alloc(struct control);
ZERO(cend[stackp]->co_next, struct control); ZERO(cend[stackp]->co_next, struct control);
cend[stackp]->co_next->co_prev = cend[stackp]; cend[stackp]->co_next->co_prev = cend[stackp];
cend[stackp]->co_next->co_parent = cend[stackp]->co_parent; cend[stackp]->co_next->co_parent = cend[stackp]->co_parent;
cend[stackp] = cend[stackp]->co_next; cend[stackp] = cend[stackp]->co_next;
} else if (!cend[stackp]) { } else if (!cend[stackp]) {
control[stackp] = cend[stackp] = alloc(struct control); control[stackp] = cend[stackp] = alloc(struct control);
ZERO(cend[stackp], struct control); ZERO(cend[stackp], struct control);
} }
if (eq(wlist->wl_word, "while")) { if (eq(wlist->wl_word, "while")) {
cend[stackp]->co_type = CO_WHILE; cend[stackp]->co_type = CO_WHILE;
cend[stackp]->co_cond = wl_copy(wlist->wl_next); /* va, wl_copy */ cend[stackp]->co_cond = wl_copy(wlist->wl_next); /* va, wl_copy */
if (!cend[stackp]->co_cond) { if (!cend[stackp]->co_cond) {
fprintf(stderr, fprintf(stderr,
"Error: missing while condition, 'false' will be assumed.\n"); "Error: missing while condition, 'false' will be assumed.\n");
} }
newblock; newblock;
} else if (eq(wlist->wl_word, "dowhile")) { } else if (eq(wlist->wl_word, "dowhile")) {
cend[stackp]->co_type = CO_DOWHILE; cend[stackp]->co_type = CO_DOWHILE;
cend[stackp]->co_cond = wl_copy(wlist->wl_next); /* va, wl_copy */ cend[stackp]->co_cond = wl_copy(wlist->wl_next); /* va, wl_copy */
if (!cend[stackp]->co_cond) { if (!cend[stackp]->co_cond) {
/* va: prevent misinterpretation as trigraph sequence with \-sign */ /* va: prevent misinterpretation as trigraph sequence with \-sign */
fprintf(stderr, fprintf(stderr,
"Error: missing dowhile condition, '?\?\?' will be assumed.\n"); "Error: missing dowhile condition, '?\?\?' will be assumed.\n");
} }
newblock; newblock;
} else if (eq(wlist->wl_word, "repeat")) { } else if (eq(wlist->wl_word, "repeat")) {
cend[stackp]->co_type = CO_REPEAT; cend[stackp]->co_type = CO_REPEAT;
if (!wlist->wl_next) { if (!wlist->wl_next) {
cend[stackp]->co_numtimes = -1; cend[stackp]->co_numtimes = -1;
} else { } else {
char *s; char *s;
double *dd; double *dd;
struct wordlist *t; /*CDHW*/
/*CDHW wlist = cp_variablesubst(cp_bquote(cp_doglob(wl_copy(wlist)))); Wrong order? Leak? CDHW*/
t = cp_doglob(cp_bquote(cp_variablesubst(wl_copy(wlist)))); /*CDHW leak from cp_doglob? */
s = t->wl_next->wl_word;
dd = ft_numparse(&s, FALSE);
if (dd) {
if (*dd < 0) {
fprintf(cp_err,
"Error: can't repeat a negative number of times\n");
*dd = 0.0;
}
cend[stackp]->co_numtimes = (int) *dd;
} else
fprintf(cp_err,
"Error: bad repeat argument %s\n",
t->wl_next->wl_word); /* CDHW */
wl_free(t); t = NULL; /* CDHW */
}
newblock;
} else if (eq(wlist->wl_word, "if")) {
cend[stackp]->co_type = CO_IF;
cend[stackp]->co_cond = wl_copy(wlist->wl_next); /* va, wl_copy */
if (!cend[stackp]->co_cond) {
fprintf(stderr,
"Error: missing if condition.\n");
}
newblock;
} else if (eq(wlist->wl_word, "foreach")) {
cend[stackp]->co_type = CO_FOREACH;
if (wlist->wl_next) {
wlist = wlist->wl_next;
cend[stackp]->co_foreachvar =
copy(wlist->wl_word);
wlist = wlist->wl_next;
} else
fprintf(stderr,
"Error: missing foreach variable.\n");
wlist = cp_doglob(wlist); /*CDHW Possible leak around here? */
cend[stackp]->co_text = wl_copy(wlist);
newblock;
} else if (eq(wlist->wl_word, "label")) {
cend[stackp]->co_type = CO_LABEL;
if (wlist->wl_next) {
cend[stackp]->co_text = wl_copy(wlist->wl_next);
/* I think of everything, don't I? */
cp_addkword(CT_LABEL, wlist->wl_next->wl_word);
if (wlist->wl_next->wl_next)
fprintf(cp_err,
"Warning: ignored extra junk after label.\n");
} else
fprintf(stderr, "Error: missing label.\n");
} else if (eq(wlist->wl_word, "goto")) {
/* Incidentally, this won't work if the values 1 and 2 ever get
* to be valid character pointers -- I think it's reasonably
* safe to assume they aren't... */
cend[stackp]->co_type = CO_GOTO;
if (wlist->wl_next) {
cend[stackp]->co_text = wl_copy(wlist->wl_next);
if (wlist->wl_next->wl_next)
fprintf(cp_err,
"Warning: ignored extra junk after goto.\n");
} else
fprintf(stderr, "Error: missing label.\n");
} else if (eq(wlist->wl_word, "continue")) {
cend[stackp]->co_type = CO_CONTINUE;
if (wlist->wl_next) {
cend[stackp]->co_numtimes = scannum(wlist->
wl_next->wl_word);
if (wlist->wl_next->wl_next)
fprintf(cp_err,
"Warning: ignored extra junk after continue %d.\n",
cend[stackp]->co_numtimes);
} else
cend[stackp]->co_numtimes = 1;
} else if (eq(wlist->wl_word, "break")) {
cend[stackp]->co_type = CO_BREAK;
if (wlist->wl_next) {
cend[stackp]->co_numtimes = scannum(wlist->
wl_next->wl_word);
if (wlist->wl_next->wl_next)
fprintf(cp_err,
"Warning: ignored extra junk after break %d.\n",
cend[stackp]->co_numtimes);
} else
cend[stackp]->co_numtimes = 1;
} else if (eq(wlist->wl_word, "end")) {
/* Throw away this thing. */
if (!cend[stackp]->co_parent) {
fprintf(stderr, "Error: no block to end.\n");
cend[stackp]->co_type = CO_UNFILLED;
} else if (cend[stackp]->co_prev) {
cend[stackp]->co_prev->co_next = NULL;
x = cend[stackp];
cend[stackp] = cend[stackp]->co_parent;
tfree(x); x=NULL;
} else {
x = cend[stackp];
cend[stackp] = cend[stackp]->co_parent;
cend[stackp]->co_children = NULL;
tfree(x); x=NULL;
}
} else if (eq(wlist->wl_word, "else")) {
if (!cend[stackp]->co_parent ||
(cend[stackp]->co_parent->co_type !=
CO_IF)) {
fprintf(stderr, "Error: misplaced else.\n");
cend[stackp]->co_type = CO_UNFILLED;
} else {
if (cend[stackp]->co_prev)
cend[stackp]->co_prev->co_next = NULL;
else
cend[stackp]->co_parent->co_children = NULL;
cend[stackp]->co_parent->co_elseblock = cend[stackp];
cend[stackp]->co_prev = NULL;
}
} else {
cend[stackp]->co_type = CO_STATEMENT;
cend[stackp]->co_text = wl_copy(wlist);
}
if (!cend[stackp]->co_parent) { struct wordlist *t; /*CDHW*/
x = cend[stackp]; /*CDHW wlist = cp_variablesubst(cp_bquote(cp_doglob(wl_copy(wlist)))); Wrong order? Leak? CDHW*/
/* We have to toss this do-while loop in here so t = cp_doglob(cp_bquote(cp_variablesubst(wl_copy(wlist)))); /*CDHW leak from cp_doglob? */
* that gotos at the top level will work. s = t->wl_next->wl_word;
*/
do { nn = 0; /* CDHW */ dd = ft_numparse(&s, FALSE);
i = doblock(x, &nn); if (dd) {
switch (*i) { if (*dd < 0) {
case NORMAL: fprintf(cp_err,
break; "Error: can't repeat a negative number of times\n");
case BROKEN: *dd = 0.0;
fprintf(cp_err, }
"Error: break not in loop or too many break levels given\n"); cend[stackp]->co_numtimes = (int) *dd;
break; } else
case CONTINUED: fprintf(cp_err,
fprintf(cp_err, "Error: bad repeat argument %s\n",
"Error: continue not in loop or too many continue levels given\n"); t->wl_next->wl_word); /* CDHW */
break; wl_free(t);
default: t = NULL; /* CDHW */
x = findlabel(i, control[stackp]); }
if (!x) newblock;
fprintf(cp_err, "Error: label %s not found\n", i);
} } else if (eq(wlist->wl_word, "if")) {
if (x) cend[stackp]->co_type = CO_IF;
x = x->co_next; cend[stackp]->co_cond = wl_copy(wlist->wl_next); /* va, wl_copy */
} while (x); if (!cend[stackp]->co_cond) {
} fprintf(stderr,
wl_free(wlist); wlist = NULL; "Error: missing if condition.\n");
if (string) { }
return (1); /* The return value is irrelevant. */ newblock;
}
} else if (eq(wlist->wl_word, "foreach")) {
cend[stackp]->co_type = CO_FOREACH;
if (wlist->wl_next) {
wlist = wlist->wl_next;
cend[stackp]->co_foreachvar =
copy(wlist->wl_word);
wlist = wlist->wl_next;
} else
fprintf(stderr,
"Error: missing foreach variable.\n");
wlist = cp_doglob(wlist); /*CDHW Possible leak around here? */
cend[stackp]->co_text = wl_copy(wlist);
newblock;
} else if (eq(wlist->wl_word, "label")) {
cend[stackp]->co_type = CO_LABEL;
if (wlist->wl_next) {
cend[stackp]->co_text = wl_copy(wlist->wl_next);
/* I think of everything, don't I? */
cp_addkword(CT_LABEL, wlist->wl_next->wl_word);
if (wlist->wl_next->wl_next)
fprintf(cp_err,
"Warning: ignored extra junk after label.\n");
} else
fprintf(stderr, "Error: missing label.\n");
} else if (eq(wlist->wl_word, "goto")) {
/* Incidentally, this won't work if the values 1 and 2 ever get
* to be valid character pointers -- I think it's reasonably
* safe to assume they aren't... */
cend[stackp]->co_type = CO_GOTO;
if (wlist->wl_next) {
cend[stackp]->co_text = wl_copy(wlist->wl_next);
if (wlist->wl_next->wl_next)
fprintf(cp_err,
"Warning: ignored extra junk after goto.\n");
} else
fprintf(stderr, "Error: missing label.\n");
} else if (eq(wlist->wl_word, "continue")) {
cend[stackp]->co_type = CO_CONTINUE;
if (wlist->wl_next) {
cend[stackp]->co_numtimes = scannum(wlist->
wl_next->wl_word);
if (wlist->wl_next->wl_next)
fprintf(cp_err,
"Warning: ignored extra junk after continue %d.\n",
cend[stackp]->co_numtimes);
} else
cend[stackp]->co_numtimes = 1;
} else if (eq(wlist->wl_word, "break")) {
cend[stackp]->co_type = CO_BREAK;
if (wlist->wl_next) {
cend[stackp]->co_numtimes = scannum(wlist->
wl_next->wl_word);
if (wlist->wl_next->wl_next)
fprintf(cp_err,
"Warning: ignored extra junk after break %d.\n",
cend[stackp]->co_numtimes);
} else
cend[stackp]->co_numtimes = 1;
} else if (eq(wlist->wl_word, "end")) {
/* Throw away this thing. */
if (!cend[stackp]->co_parent) {
fprintf(stderr, "Error: no block to end.\n");
cend[stackp]->co_type = CO_UNFILLED;
} else if (cend[stackp]->co_prev) {
cend[stackp]->co_prev->co_next = NULL;
x = cend[stackp];
cend[stackp] = cend[stackp]->co_parent;
tfree(x);
x=NULL;
} else {
x = cend[stackp];
cend[stackp] = cend[stackp]->co_parent;
cend[stackp]->co_children = NULL;
tfree(x);
x=NULL;
}
} else if (eq(wlist->wl_word, "else")) {
if (!cend[stackp]->co_parent ||
(cend[stackp]->co_parent->co_type !=
CO_IF)) {
fprintf(stderr, "Error: misplaced else.\n");
cend[stackp]->co_type = CO_UNFILLED;
} else {
if (cend[stackp]->co_prev)
cend[stackp]->co_prev->co_next = NULL;
else
cend[stackp]->co_parent->co_children = NULL;
cend[stackp]->co_parent->co_elseblock = cend[stackp];
cend[stackp]->co_prev = NULL;
}
} else {
cend[stackp]->co_type = CO_STATEMENT;
cend[stackp]->co_text = wl_copy(wlist);
}
if (!cend[stackp]->co_parent) {
x = cend[stackp];
/* We have to toss this do-while loop in here so
* that gotos at the top level will work.
*/
do {
nn = 0; /* CDHW */
i = doblock(x, &nn);
switch (*i) {
case NORMAL:
break;
case BROKEN:
fprintf(cp_err,
"Error: break not in loop or too many break levels given\n");
break;
case CONTINUED:
fprintf(cp_err,
"Error: continue not in loop or too many continue levels given\n");
break;
default:
x = findlabel(i, control[stackp]);
if (!x)
fprintf(cp_err, "Error: label %s not found\n", i);
}
if (x)
x = x->co_next;
} while (x);
}
wl_free(wlist);
wlist = NULL;
if (string) {
return (1); /* The return value is irrelevant. */
}
} }
wl_free(wlist); wlist = NULL; wl_free(wlist);
wlist = NULL;
return (0); /* va: which value? */ return (0); /* va: which value? */
} }
@ -919,12 +933,12 @@ cp_toplevel(void)
/* va: This totally frees the control structures */ /* va: This totally frees the control structures */
void cp_free_control(void) void cp_free_control(void)
{ {
int i; int i;
for (i=stackp; i>=0; i--) ctl_free(control[i]); for (i=stackp; i>=0; i--) ctl_free(control[i]);
control[0] = cend[0] = NULL; control[0] = cend[0] = NULL;
stackp = 0; stackp = 0;
} }

24
src/frontend/trannoise/1-f-code.c
View File

@ -21,25 +21,25 @@
#include "wallace.h" #include "wallace.h"
void f_alpha(int n_pts, int n_exp, float X[], float Q_d, void f_alpha(int n_pts, int n_exp, double X[], double Q_d,
float alpha) double alpha)
{ {
int i; int i;
float *hfa, *wfa; double *hfa, *wfa;
float ha; double ha;
ha = alpha/2.0f ; ha = alpha/2.0f ;
// Q_d = sqrt(Q_d); /* find the deviation of the noise */ // Q_d = sqrt(Q_d); /* find the deviation of the noise */
hfa = TMALLOC(float,n_pts); hfa = TMALLOC(double,n_pts);
wfa = TMALLOC(float,n_pts); wfa = TMALLOC(double,n_pts);
hfa[0] = 1.0f; hfa[0] = 1.0f;
wfa[0] = Q_d * (float)GaussWa; wfa[0] = Q_d * GaussWa;
/* generate the coefficients hk */ /* generate the coefficients hk */
for (i=1 ; i < n_pts; i++) { for (i=1 ; i < n_pts; i++) {
/* generate the coefficients hk */ /* generate the coefficients hk */
hfa[i] = hfa[i-1] * (ha + (float)(i-1)) / ( (float)(i) ); hfa[i] = hfa[i-1] * (ha + (double)(i-1)) / ( (double)(i) );
/* fill the sequence wk with white noise */ /* fill the sequence wk with white noise */
wfa[i] = Q_d * (float)GaussWa; wfa[i] = Q_d * GaussWa;
} }
// for (i=0 ; i < n_pts; i++) // for (i=0 ; i < n_pts; i++)
@ -88,13 +88,13 @@ trnoise_state_gen(struct trnoise_state *this, CKTcircuit *ckt)
newexp++; newexp++;
} }
this->oneof = TMALLOC(float, newsteps); this->oneof = TMALLOC(double, newsteps);
this->oneof_length = newsteps; this->oneof_length = newsteps;
f_alpha((int) newsteps, newexp, f_alpha((int) newsteps, newexp,
this -> oneof, this -> oneof,
(float) this -> NAMP, this -> NAMP,
(float) this -> NALPHA); this -> NALPHA);
} }
trnoise_state_push(this, 0.0); /* first is deterministic */ trnoise_state_push(this, 0.0); /* first is deterministic */

6
src/include/1-f-code.h
View File

@ -4,8 +4,8 @@ extern double exprand(double);
extern int poisson(double); extern int poisson(double);
extern double gauss0(void); extern double gauss0(void);
void f_alpha(int n_pts, int n_exp, float X[], float Q_d, void f_alpha(int n_pts, int n_exp, double X[], double Q_d,
float alpha); double alpha);
void rvfft(float X[], unsigned long int n); void rvfft(float X[], unsigned long int n);
@ -18,7 +18,7 @@ struct trnoise_state
double NA, TS, NAMP, NALPHA, RTSAM, RTSCAPT, RTSEMT; double NA, TS, NAMP, NALPHA, RTSAM, RTSCAPT, RTSEMT;
float *oneof; double *oneof;
size_t oneof_length; size_t oneof_length;
double RTScapTime, RTSemTime; double RTScapTime, RTSemTime;

62
src/include/fftext.h
View File

@ -19,7 +19,7 @@ int fftInit(long M);
void fftFree(void); void fftFree(void);
// release storage for all private cosine and bit reversed tables // release storage for all private cosine and bit reversed tables
void ffts(float *data, long M, long Rows); void ffts(double *data, long M, long Rows);
/* Compute in-place complex fft on the rows of the input array */ /* Compute in-place complex fft on the rows of the input array */
/* INPUTS */ /* INPUTS */
/* *ioptr = input data array */ /* *ioptr = input data array */
@ -28,7 +28,7 @@ void ffts(float *data, long M, long Rows);
/* OUTPUTS */ /* OUTPUTS */
/* *ioptr = output data array */ /* *ioptr = output data array */
void iffts(float *data, long M, long Rows); void iffts(double *data, long M, long Rows);
/* Compute in-place inverse complex fft on the rows of the input array */ /* Compute in-place inverse complex fft on the rows of the input array */
/* INPUTS */ /* INPUTS */
/* *ioptr = input data array */ /* *ioptr = input data array */
@ -37,7 +37,7 @@ void iffts(float *data, long M, long Rows);
/* OUTPUTS */ /* OUTPUTS */
/* *ioptr = output data array */ /* *ioptr = output data array */
void rffts(float *data, long M, long Rows); void rffts(double *data, long M, long Rows);
/* Compute in-place real fft on the rows of the input array */ /* Compute in-place real fft on the rows of the input array */
/* The result is the complex spectra of the positive frequencies */ /* The result is the complex spectra of the positive frequencies */
/* except the location for the first complex number contains the real */ /* except the location for the first complex number contains the real */
@ -51,7 +51,7 @@ void rffts(float *data, long M, long Rows);
/* *ioptr = output data array in the following order */ /* *ioptr = output data array in the following order */
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */ /* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
void riffts(float *data, long M, long Rows); void riffts(double *data, long M, long Rows);
/* Compute in-place real ifft on the rows of the input array */ /* Compute in-place real ifft on the rows of the input array */
/* data order as from rffts */ /* data order as from rffts */
/* INPUTS */ /* INPUTS */
@ -62,7 +62,7 @@ void riffts(float *data, long M, long Rows);
/* OUTPUTS */ /* OUTPUTS */
/* *ioptr = real output data array */ /* *ioptr = real output data array */
void rspectprod(float *data1, float *data2, float *outdata, long N); void rspectprod(double *data1, double *data2, double *outdata, long N);
// When multiplying a pair of spectra from rfft care must be taken to multiply the // When multiplying a pair of spectra from rfft care must be taken to multiply the
// two real values seperately from the complex ones. This routine does it correctly. // two real values seperately from the complex ones. This routine does it correctly.
// the result can be stored in-place over one of the inputs // the result can be stored in-place over one of the inputs
@ -74,35 +74,33 @@ void rspectprod(float *data1, float *data2, float *outdata, long N);
/* *outdata = output data array spectra */ /* *outdata = output data array spectra */
/* The following is FYI // The following is FYI
Note that most of the fft routines require full matrices, ie Rsiz==Ncols //Note that most of the fft routines require full matrices, ie Rsiz==Ncols
This is how I like to define a real matrix: //This is how I like to define a real matrix:
struct matrix { // real matrix //struct matrix { // real matrix
float *d; // pointer to data // double *d; // pointer to data
long Nrows; // number of rows in the matrix // long Nrows; // number of rows in the matrix
long Ncols; // number of columns in the matrix (can be less than Rsiz) // long Ncols; // number of columns in the matrix (can be less than Rsiz)
long Rsiz; // number of floats from one row to the next // long Rsiz; // number of doubles from one row to the next
}; //};
typedef struct matrix matrix; //typedef struct matrix matrix;
CACHEFILLMALLOC and CEILCACHELINE can be used instead of malloc to make // CACHEFILLMALLOC and CEILCACHELINE can be used instead of malloc to make
arrays that start exactly on a cache line start. // arrays that start exactly on a cache line start.
First we CACHEFILLMALLOC a void * (use this void * when free'ing), // First we CACHEFILLMALLOC a void * (use this void * when free'ing),
then we set our array pointer equal to the properly cast CEILCACHELINE of this void * // then we set our array pointer equal to the properly cast CEILCACHELINE of this void *
example: // example:
aInit = CACHEFILLMALLOC( NUMFLOATS*sizeof(float) ); // aInit = CACHEFILLMALLOC( NUMFLOATS*sizeof(float) );
a = (float *) CEILCACHELINE(ainit); // a = (float *) CEILCACHELINE(ainit);
... main body of code ... // ... main body of code ...
free(aInit); // free(aInit);
//
To disable this alignment, set CACHELINESIZE to 1 // To disable this alignment, set CACHELINESIZE to 1
#define CACHELINESIZE 32 // Bytes per cache line //#define CACHELINESIZE 32 // Bytes per cache line
#define CACHELINEFILL (CACHELINESIZE-1) //#define CACHELINEFILL (CACHELINESIZE-1)
#define CEILCACHELINE(p) ((((unsigned long)p+CACHELINEFILL)/CACHELINESIZE)*CACHELINESIZE) //#define CEILCACHELINE(p) ((((unsigned long)p+CACHELINEFILL)/CACHELINESIZE)*CACHELINESIZE)
#define CACHEFILLMALLOC(n) malloc((n)+CACHELINEFILL) //#define CACHEFILLMALLOC(n) malloc((n)+CACHELINEFILL)
*/

106
src/maths/fft/fftext.h
View File

@ -1,106 +0,0 @@
/*******************************************************************
This file extends the fftlib with calls to maintain the cosine and bit reversed tables
for you (including mallocs and free's). Call the init routine for each fft size you will
be using. Then you can call the fft routines below which will make the fftlib library
call with the appropriate tables passed. When you are done with all fft's you can call
fftfree to release the storage for the tables. Note that you can call fftinit repeatedly
with the same size, the extra calls will be ignored. So, you could make a macro to call
fftInit every time you call ffts. For example you could have someting like:
#define FFT(a,n) if(!fftInit(roundtol(LOG2(n)))) ffts(a,roundtol(LOG2(n)),1);else printf("fft error\n");
*******************************************************************/
int fftInit(long M);
// malloc and init cosine and bit reversed tables for a given size fft, ifft, rfft, rifft
/* INPUTS */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* OUTPUTS */
/* private cosine and bit reversed tables */
void fftFree(void);
// release storage for all private cosine and bit reversed tables
void ffts(double *data, long M, long Rows);
/* Compute in-place complex fft on the rows of the input array */
/* INPUTS */
/* *ioptr = input data array */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* Rows = number of rows in ioptr array (use 1 for Rows for a single fft) */
/* OUTPUTS */
/* *ioptr = output data array */
void iffts(double *data, long M, long Rows);
/* Compute in-place inverse complex fft on the rows of the input array */
/* INPUTS */
/* *ioptr = input data array */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* Rows = number of rows in ioptr array (use 1 for Rows for a single fft) */
/* OUTPUTS */
/* *ioptr = output data array */
void rffts(double *data, long M, long Rows);
/* Compute in-place real fft on the rows of the input array */
/* The result is the complex spectra of the positive frequencies */
/* except the location for the first complex number contains the real */
/* values for DC and Nyquest */
/* See rspectprod for multiplying two of these spectra together- ex. for fast convolution */
/* INPUTS */
/* *ioptr = real input data array */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* Rows = number of rows in ioptr array (use 1 for Rows for a single fft) */
/* OUTPUTS */
/* *ioptr = output data array in the following order */
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
void riffts(double *data, long M, long Rows);
/* Compute in-place real ifft on the rows of the input array */
/* data order as from rffts */
/* INPUTS */
/* *ioptr = input data array in the following order */
/* M = log2 of fft size (ex M=10 for 1024 point fft) */
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
/* Rows = number of rows in ioptr array (use 1 for Rows for a single fft) */
/* OUTPUTS */
/* *ioptr = real output data array */
void rspectprod(double *data1, double *data2, double *outdata, long N);
// When multiplying a pair of spectra from rfft care must be taken to multiply the
// two real values seperately from the complex ones. This routine does it correctly.
// the result can be stored in-place over one of the inputs
/* INPUTS */
/* *data1 = input data array first spectra */
/* *data2 = input data array second spectra */
/* N = fft input size for both data1 and data2 */
/* OUTPUTS */
/* *outdata = output data array spectra */
// The following is FYI
//Note that most of the fft routines require full matrices, ie Rsiz==Ncols
//This is how I like to define a real matrix:
//struct matrix { // real matrix
// double *d; // pointer to data
// long Nrows; // number of rows in the matrix
// long Ncols; // number of columns in the matrix (can be less than Rsiz)
// long Rsiz; // number of doubles from one row to the next
//};
//typedef struct matrix matrix;
// CACHEFILLMALLOC and CEILCACHELINE can be used instead of malloc to make
// arrays that start exactly on a cache line start.
// First we CACHEFILLMALLOC a void * (use this void * when free'ing),
// then we set our array pointer equal to the properly cast CEILCACHELINE of this void *
// example:
// aInit = CACHEFILLMALLOC( NUMFLOATS*sizeof(float) );
// a = (float *) CEILCACHELINE(ainit);
// ... main body of code ...
// free(aInit);
//
// To disable this alignment, set CACHELINESIZE to 1
//#define CACHELINESIZE 32 // Bytes per cache line
//#define CACHELINEFILL (CACHELINESIZE-1)
//#define CEILCACHELINE(p) ((((unsigned long)p+CACHELINEFILL)/CACHELINESIZE)*CACHELINESIZE)
//#define CACHEFILLMALLOC(n) malloc((n)+CACHELINEFILL)