diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..d68b314 --- /dev/null +++ b/.clang-format @@ -0,0 +1,24 @@ +BasedOnStyle: LLVM +IndentWidth: 2 +UseTab: Never +IndentCaseLabels: true +AllowShortBlocksOnASingleLine: true +AllowShortCaseLabelsOnASingleLine: true +AllowShortIfStatementsOnASingleLine: true +AllowShortFunctionsOnASingleLine: true +ColumnLimit: 100 +Standard: Cpp03 +SortIncludes: false +BreakBeforeBraces: Custom +BraceWrapping: + AfterClass: false + AfterControlStatement: false + AfterEnum: false + AfterFunction: false + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false + BeforeCatch: true + BeforeElse: true + IndentBraces: false diff --git a/check.h b/check.h index 99425cf..6cdc5fa 100644 --- a/check.h +++ b/check.h @@ -20,17 +20,19 @@ /* some useful error-checking macros: */ -#define CHECK(condition, message) do { \ - if (!(condition)) { \ - fprintf(stderr, "harminv: failure on line %d of " __FILE__ ": " \ - message "\n", __LINE__); exit(EXIT_FAILURE); \ - } \ -} while (0) +#define CHECK(condition, message) \ + do { \ + if (!(condition)) { \ + fprintf(stderr, "harminv: failure on line %d of " __FILE__ ": " message "\n", __LINE__); \ + exit(EXIT_FAILURE); \ + } \ + } while (0) -#define CHK_MALLOC(p, t, n) do { \ - size_t CHK_MALLOC_n_tmp = (n); \ - (p) = (t *) malloc(sizeof(t) * CHK_MALLOC_n_tmp); \ - CHECK((p) || CHK_MALLOC_n_tmp == 0, "out of memory!"); \ -} while (0) +#define CHK_MALLOC(p, t, n) \ + do { \ + size_t CHK_MALLOC_n_tmp = (n); \ + (p) = (t *)malloc(sizeof(t) * CHK_MALLOC_n_tmp); \ + CHECK((p) || CHK_MALLOC_n_tmp == 0, "out of memory!"); \ + } while (0) #endif /* CHECK_H */ diff --git a/harminv-int.h b/harminv-int.h index b61c32e..2030b8b 100644 --- a/harminv-int.h +++ b/harminv-int.h @@ -25,7 +25,7 @@ without it. Alternatively, use the complex STL class in C++. */ -# include +#include #endif #include "harminv.h" @@ -34,45 +34,44 @@ typedef harminv_complex cmplx; /* shortcut */ /* define extended precision (long double) cmplxl type */ #if defined(HAVE_LONG_DOUBLE) && !defined(__cplusplus) -# if defined(__cplusplus) -# include +#if defined(__cplusplus) +#include typedef std::complex cmplxl; /* problematic w/some C++ comp's */ -# elif defined(_Complex_I) && defined(complex) && defined(I) +#elif defined(_Complex_I) && defined(complex) && defined(I) /* C99 header was included before harminv.h */ typedef long double _Complex cmplxl; -# endif +#endif #else typedef cmplx cmplxl; #endif - #ifdef __cplusplus -# define I cmplx(0,1) -# define creal(c) real(c) -# define cimag(c) imag(c) -# define cabs(c) abs(c) -# define carg(c) arg(c) -# define cexp(c) exp(c) -# define csqrt(c) sqrt(c) -# define clog(c) log(c) +#define I cmplx(0, 1) +#define creal(c) real(c) +#define cimag(c) imag(c) +#define cabs(c) abs(c) +#define carg(c) arg(c) +#define cexp(c) exp(c) +#define csqrt(c) sqrt(c) +#define clog(c) log(c) #else -# ifndef HAVE_CARG /* Cray doesn't have this for some reason */ -# define carg(c) atan2(cimag(c), creal(c)) -# endif +#ifndef HAVE_CARG /* Cray doesn't have this for some reason */ +#define carg(c) atan2(cimag(c), creal(c)) +#endif #endif struct harminv_data_struct { - const cmplx *c; - int n, K, J, nfreqs; - double fmin, fmax; - cmplx *z; - cmplx *U0, *U1; - cmplxl *G0, *G0_M, *D0; /* cached G, G_M, and D arrays for U0 */ - cmplx *B, *u; /* eigen-solutions of U1*B = u*U0*B */ - cmplx *amps; /* mode amplitudes */ - double *errs; /* relative "error" estimates */ + const cmplx *c; + int n, K, J, nfreqs; + double fmin, fmax; + cmplx *z; + cmplx *U0, *U1; + cmplxl *G0, *G0_M, *D0; /* cached G, G_M, and D arrays for U0 */ + cmplx *B, *u; /* eigen-solutions of U1*B = u*U0*B */ + cmplx *amps; /* mode amplitudes */ + double *errs; /* relative "error" estimates */ }; #endif /* HARMINV_INT_H */ diff --git a/harminv-main.c b/harminv-main.c index e7929e8..2b40e4b 100644 --- a/harminv-main.c +++ b/harminv-main.c @@ -25,100 +25,92 @@ #include "copyright.h" #ifdef HAVE_UNISTD_H -# include +#include #endif #ifdef HAVE_GETOPT_H -# include +#include #endif /* eat whitespace, including #... comments, from the file. Returns the number of newlines read (so that a line count can be maintained). If echo_comments != 0, then echo #... comments to stdout. Commas count as whitespace, so that we can read comma-delimited text. */ -static int eat_whitespace(FILE *f, int echo_comments) -{ - int c, newlines = 0; - do { - do { - c = getc(f); - newlines += c == '\n'; - } while (isspace(c) || c == ','); - - if (c == '#') { /* # begins comments that extend to the newline */ - if (echo_comments) - putc(c, stdout); - do { - c = getc(f); - if (echo_comments) { - if (c != EOF) - putc(c, stdout); - else /* terminate line if we hit EOF */ - putc('\n', stdout); - } - newlines += c == '\n'; - } while (c != EOF && c != '\n'); - } - } while (isspace (c)); - - ungetc(c, f); /* put back the last character read */ - newlines -= c == '\n'; - - return newlines; +static int eat_whitespace(FILE *f, int echo_comments) { + int c, newlines = 0; + do { + do { + c = getc(f); + newlines += c == '\n'; + } while (isspace(c) || c == ','); + + if (c == '#') { /* # begins comments that extend to the newline */ + if (echo_comments) putc(c, stdout); + do { + c = getc(f); + if (echo_comments) { + if (c != EOF) + putc(c, stdout); + else /* terminate line if we hit EOF */ + putc('\n', stdout); + } + newlines += c == '\n'; + } while (c != EOF && c != '\n'); + } + } while (isspace(c)); + + ungetc(c, f); /* put back the last character read */ + newlines -= c == '\n'; + + return newlines; } -static int eat_plus(FILE *f) -{ - int c = getc(f); - if (c != EOF && c != '+') - ungetc(c, f); - return (c == '+' || c == '-'); +static int eat_plus(FILE *f) { + int c = getc(f); + if (c != EOF && c != '+') ungetc(c, f); + return (c == '+' || c == '-'); } -static int eat_i(FILE *f) -{ - int c = getc(f); - if (c != EOF && tolower(c) != 'i') - ungetc(c, f); - return (tolower(c) == 'i'); +static int eat_i(FILE *f) { + int c = getc(f); + if (c != EOF && tolower(c) != 'i') ungetc(c, f); + return (tolower(c) == 'i'); } -static cmplx *read_input_data(FILE *f, int *n, int verbose) -{ - cmplx *data = NULL; - int line = 1, n_alloc = 0; - *n = 0; - - do { - double re=0.0, im=0.0; - int nread; - line += eat_whitespace(f, verbose); - nread = fscanf(f, "%lg", &re); - if (nread == 1 && eat_plus(f)) { - nread = fscanf(f, "%lg", &im); - if (nread == 1) nread = eat_i(f); - } - else - im = 0.0; - if (nread != EOF) { - if (nread < 1) { - fprintf(stderr, "harminv: invalid input on line %d.\n", - line); - free(data); - *n = 0; - return NULL; - } - if (*n >= n_alloc) { - n_alloc = (n_alloc + 1) * 2; - data = (cmplx*) realloc(data, sizeof(cmplx) * n_alloc); - CHECK(data, "out of memory"); - } - data[*n] = re + I*im; - ++*n; - } - } while (!feof(f)); - - data = (cmplx*) realloc(data, sizeof(cmplx) * *n); - return data; +static cmplx *read_input_data(FILE *f, int *n, int verbose) { + cmplx *data = NULL; + int line = 1, n_alloc = 0; + *n = 0; + + do { + double re = 0.0, im = 0.0; + int nread; + line += eat_whitespace(f, verbose); + nread = fscanf(f, "%lg", &re); + if (nread == 1 && eat_plus(f)) { + nread = fscanf(f, "%lg", &im); + if (nread == 1) nread = eat_i(f); + } + else + im = 0.0; + if (nread != EOF) { + if (nread < 1) { + fprintf(stderr, "harminv: invalid input on line %d.\n", line); + free(data); + *n = 0; + return NULL; + } + if (*n >= n_alloc) { + n_alloc = (n_alloc + 1) * 2; + data = (cmplx *)realloc(data, sizeof(cmplx) * n_alloc); + CHECK(data, "out of memory"); + } + data[*n] = re + I * im; + ++*n; + } + } while (!feof(f)); + + data = (cmplx *)realloc(data, sizeof(cmplx) * *n); + return data; } static const double inf = HUGE_VAL; @@ -132,364 +124,303 @@ static const double inf = HUGE_VAL; #define REL_AMP_THRESH -1.0 #define Q_THRESH 10.0 -static void usage(FILE *f) -{ - fprintf(f, "Usage: harminv [options] -...\n" - "Options: \n" - " -h : this help message\n" - " -V : print version number and copyright\n" - " -v : verbose output\n" - " -T : specify periods instead of frequencies\n" - " -w : specify/output angular frequency, not frequency\n" - " -n : flip the sign of the frequency convention\n" - " -t
: specify sampling interval dt [default: 1]\n" - " -d : specify spectral density [default: %g]\n" - " -f : use at least basis functions [default: %d]\n" - " -s : sort by = freq/err/decay/amp [default: freq]\n" - " -F : discard frequencies outside of specified range\n" - " -a : discard amplitudes < max * [default: %e]\n" - " -A : discard amplitudes < [default: %g]\n" - " -e : discard relative errors > min * [default: %e]\n" - " -E : discard relative errors > [default: %e]\n" - " -Q : discard Q > [default: %g]\n", - DENSITY, - NFMIN, - AMP_THRESH, REL_AMP_THRESH, - ERR_THRESH, REL_ERR_THRESH, - Q_THRESH); +static void usage(FILE *f) { + fprintf(f, + "Usage: harminv [options] -...\n" + "Options: \n" + " -h : this help message\n" + " -V : print version number and copyright\n" + " -v : verbose output\n" + " -T : specify periods instead of frequencies\n" + " -w : specify/output angular frequency, not frequency\n" + " -n : flip the sign of the frequency convention\n" + " -t
: specify sampling interval dt [default: 1]\n" + " -d : specify spectral density [default: %g]\n" + " -f : use at least basis functions [default: %d]\n" + " -s : sort by = freq/err/decay/amp [default: freq]\n" + " -F : discard frequencies outside of specified range\n" + " -a : discard amplitudes < max * [default: %e]\n" + " -A : discard amplitudes < [default: %g]\n" + " -e : discard relative errors > min * [default: %e]\n" + " -E : discard relative errors > [default: %e]\n" + " -Q : discard Q > [default: %g]\n", + DENSITY, NFMIN, AMP_THRESH, REL_AMP_THRESH, ERR_THRESH, REL_ERR_THRESH, Q_THRESH); } #define TWOPI 6.2831853071795864769252867665590057683943388 harminv_data hd; -enum { - SORT_FREQUENCY, SORT_DECAY, SORT_ERROR, SORT_AMPLITUDE, SORT_Q -} sortby = SORT_FREQUENCY; - -static int cmp(double a, double b) -{ - return a > b ? 1 : (a < b ? -1 : 0); -} - -static int compar(const void *a, const void *b) -{ - const int *ia = (const int *) a; - const int *ib = (const int *) b; - cmplx aa, ab; - - switch (sortby) { - case SORT_FREQUENCY: - return cmp(harminv_get_freq(hd,*ia), harminv_get_freq(hd,*ib)); - case SORT_DECAY: - return cmp(harminv_get_decay(hd,*ia), harminv_get_decay(hd,*ib)); - case SORT_ERROR: - return cmp(harminv_get_freq_error(hd, *ia), - harminv_get_freq_error(hd, *ib)); - case SORT_AMPLITUDE: - harminv_get_amplitude(&aa, hd, *ia); - harminv_get_amplitude(&ab, hd, *ia); - return cmp(cabs(aa), cabs(ab)); - case SORT_Q: - return cmp(harminv_get_freq(hd,*ia) / harminv_get_decay(hd,*ia), - harminv_get_freq(hd,*ib) / harminv_get_decay(hd,*ib)); - } - return 0; +enum { SORT_FREQUENCY, SORT_DECAY, SORT_ERROR, SORT_AMPLITUDE, SORT_Q } sortby = SORT_FREQUENCY; + +static int cmp(double a, double b) { return a > b ? 1 : (a < b ? -1 : 0); } + +static int compar(const void *a, const void *b) { + const int *ia = (const int *)a; + const int *ib = (const int *)b; + cmplx aa, ab; + + switch (sortby) { + case SORT_FREQUENCY: return cmp(harminv_get_freq(hd, *ia), harminv_get_freq(hd, *ib)); + case SORT_DECAY: return cmp(harminv_get_decay(hd, *ia), harminv_get_decay(hd, *ib)); + case SORT_ERROR: return cmp(harminv_get_freq_error(hd, *ia), harminv_get_freq_error(hd, *ib)); + case SORT_AMPLITUDE: + harminv_get_amplitude(&aa, hd, *ia); + harminv_get_amplitude(&ab, hd, *ia); + return cmp(cabs(aa), cabs(ab)); + case SORT_Q: + return cmp(harminv_get_freq(hd, *ia) / harminv_get_decay(hd, *ia), + harminv_get_freq(hd, *ib) / harminv_get_decay(hd, *ib)); + } + return 0; } typedef struct { - int verbose; - double fmin, fmax; - int only_f_inrange; - double err_thresh, rel_err_thresh, amp_thresh, rel_amp_thresh, Q_thresh; - double min_err, max_amp; - int num_ok; + int verbose; + double fmin, fmax; + int only_f_inrange; + double err_thresh, rel_err_thresh, amp_thresh, rel_amp_thresh, Q_thresh; + double min_err, max_amp; + int num_ok; } mode_ok_data; -static int mode_ok(harminv_data d, int k, void *ok_d_) -{ - mode_ok_data *ok_d = (mode_ok_data *) ok_d_; - double errk, ampk, f; - int ok; - cmplx aa; - - if (k == -1) { /* initialize */ - int i; - ok_d->num_ok = 0; - if (!harminv_get_num_freqs(d)) - return 0; - ok_d->min_err = harminv_get_freq_error(d, 0);; - harminv_get_amplitude(&aa, d, 0); - ok_d->max_amp = cabs(aa); - for (i = 1; i < harminv_get_num_freqs(d); ++i) { - double err, amp; - if ((err = harminv_get_freq_error(d, i)) < ok_d->min_err) - ok_d->min_err = err; - harminv_get_amplitude(&aa, d, i); - if ((amp = cabs(aa)) > ok_d->max_amp) - ok_d->max_amp = amp; - - } - return 0; - } - else if (k == -2) { /* finish */ - if (ok_d->verbose && harminv_get_num_freqs(d)) - printf("# harminv: %d/%d modes are ok: " - "errs <= %e and %e * %e\n, " - "amps >= %g, %e * %g, " - "|Q| >= %g\n", - ok_d->num_ok, harminv_get_num_freqs(d), - ok_d->err_thresh, ok_d->rel_err_thresh, ok_d->min_err, - ok_d->amp_thresh, ok_d->rel_amp_thresh, ok_d->max_amp, - ok_d->Q_thresh); - return 0; - } - - f = fabs(harminv_get_freq(d, k)); - errk = harminv_get_freq_error(d, k); - harminv_get_amplitude(&aa, d, k); - ampk = cabs(aa); - - ok = ((!ok_d->only_f_inrange || (f >= ok_d->fmin && f <= ok_d->fmax)) - && errk <= ok_d->err_thresh - && (!isfinite(ok_d->rel_err_thresh) || errk <= ok_d->min_err * ok_d->rel_err_thresh) - && ampk >= ok_d->amp_thresh - && ampk >= ok_d->rel_amp_thresh * ok_d->max_amp - && fabs(harminv_get_Q(d,k)) >= ok_d->Q_thresh); - - ok_d->num_ok += ok; - - return ok; +static int mode_ok(harminv_data d, int k, void *ok_d_) { + mode_ok_data *ok_d = (mode_ok_data *)ok_d_; + double errk, ampk, f; + int ok; + cmplx aa; + + if (k == -1) { /* initialize */ + int i; + ok_d->num_ok = 0; + if (!harminv_get_num_freqs(d)) return 0; + ok_d->min_err = harminv_get_freq_error(d, 0); + ; + harminv_get_amplitude(&aa, d, 0); + ok_d->max_amp = cabs(aa); + for (i = 1; i < harminv_get_num_freqs(d); ++i) { + double err, amp; + if ((err = harminv_get_freq_error(d, i)) < ok_d->min_err) ok_d->min_err = err; + harminv_get_amplitude(&aa, d, i); + if ((amp = cabs(aa)) > ok_d->max_amp) ok_d->max_amp = amp; + } + return 0; + } + else if (k == -2) { /* finish */ + if (ok_d->verbose && harminv_get_num_freqs(d)) + printf("# harminv: %d/%d modes are ok: " + "errs <= %e and %e * %e\n, " + "amps >= %g, %e * %g, " + "|Q| >= %g\n", + ok_d->num_ok, harminv_get_num_freqs(d), ok_d->err_thresh, ok_d->rel_err_thresh, + ok_d->min_err, ok_d->amp_thresh, ok_d->rel_amp_thresh, ok_d->max_amp, ok_d->Q_thresh); + return 0; + } + + f = fabs(harminv_get_freq(d, k)); + errk = harminv_get_freq_error(d, k); + harminv_get_amplitude(&aa, d, k); + ampk = cabs(aa); + + ok = ((!ok_d->only_f_inrange || (f >= ok_d->fmin && f <= ok_d->fmax)) && + errk <= ok_d->err_thresh && + (!isfinite(ok_d->rel_err_thresh) || errk <= ok_d->min_err * ok_d->rel_err_thresh) && + ampk >= ok_d->amp_thresh && ampk >= ok_d->rel_amp_thresh * ok_d->max_amp && + fabs(harminv_get_Q(d, k)) >= ok_d->Q_thresh); + + ok_d->num_ok += ok; + + return ok; } #define SOLVE_ONCE_ONLY 0 /* 1 to use harminv_solve_once */ -#define SOLVE_OK_ONLY 0 /* 1 for experimental solver */ - -int main(int argc, char **argv) -{ - int verbose = 0; - int c; - extern char *optarg; - extern int optind; - int specify_periods = 0; - int specify_omega = 0; - int negate_omega = 0; - double dt = 1.0; - mode_ok_data ok_d; - int n, nf, nfmin = NFMIN; - double density = DENSITY; - int iarg; - cmplx *data; - - ok_d.only_f_inrange = 0; - ok_d.err_thresh = ERR_THRESH; - ok_d.rel_err_thresh = REL_ERR_THRESH; - ok_d.amp_thresh = AMP_THRESH; - ok_d.rel_amp_thresh = REL_AMP_THRESH; - ok_d.Q_thresh = Q_THRESH; - - while ((c = getopt(argc, argv, "hvVTFwnt:d:f:s:e:E:a:A:Q:")) != -1) - switch (c) { - case 'h': - usage(stdout); - return EXIT_SUCCESS; - case 'V': - printf("harminv " PACKAGE_VERSION " by Steven G. Johnson\n" - COPYRIGHT); - return EXIT_SUCCESS; - case 'v': - verbose = 1; - break; - case 'T': - specify_periods = 1; - break; - case 'w': - specify_omega = 1; - break; - case 'n': - negate_omega = 1; - break; - case 'F': - ok_d.only_f_inrange = 1; - break; - case 'a': - ok_d.rel_amp_thresh = atof(optarg); - break; - case 'A': - ok_d.amp_thresh = atof(optarg); - break; - case 'E': - ok_d.err_thresh = atof(optarg); - break; - case 'e': - ok_d.rel_err_thresh = atof(optarg); - break; - case 'Q': - ok_d.Q_thresh = atof(optarg); - break; - case 't': - dt = atof(optarg); - break; - case 'f': - nfmin = atoi(optarg); - break; - case 'd': - density = atof(optarg); - if (density < 0) { - fprintf(stderr, "harminv: -d argument must be >= 0\n"); - return EXIT_FAILURE; - } - break; - case 's': - switch (tolower(optarg[0])) { - case 'f': - sortby = SORT_FREQUENCY; - break; - case 'd': - sortby = SORT_DECAY; - break; - case 'e': - sortby = SORT_ERROR; - break; - case 'a': - sortby = SORT_AMPLITUDE; - break; - case 'q': - sortby = SORT_Q; - break; - default: - fprintf(stderr, "harminv: invalid sort type -s %c\n", tolower(optarg[0])); - usage(stderr); - return EXIT_FAILURE; - } - break; - default: - fprintf(stderr, "harminv: invalid argument -%c\n", c); - usage(stderr); - return EXIT_FAILURE; - } - if (optind == argc) { /* no parameters left */ - fprintf(stderr, "harminv: missing required frequency range(s)\n"); - usage(stderr); +#define SOLVE_OK_ONLY 0 /* 1 for experimental solver */ + +int main(int argc, char **argv) { + int verbose = 0; + int c; + extern char *optarg; + extern int optind; + int specify_periods = 0; + int specify_omega = 0; + int negate_omega = 0; + double dt = 1.0; + mode_ok_data ok_d; + int n, nf, nfmin = NFMIN; + double density = DENSITY; + int iarg; + cmplx *data; + + ok_d.only_f_inrange = 0; + ok_d.err_thresh = ERR_THRESH; + ok_d.rel_err_thresh = REL_ERR_THRESH; + ok_d.amp_thresh = AMP_THRESH; + ok_d.rel_amp_thresh = REL_AMP_THRESH; + ok_d.Q_thresh = Q_THRESH; + + while ((c = getopt(argc, argv, "hvVTFwnt:d:f:s:e:E:a:A:Q:")) != -1) + switch (c) { + case 'h': usage(stdout); return EXIT_SUCCESS; + case 'V': + printf("harminv " PACKAGE_VERSION " by Steven G. Johnson\n" COPYRIGHT); + return EXIT_SUCCESS; + case 'v': verbose = 1; break; + case 'T': specify_periods = 1; break; + case 'w': specify_omega = 1; break; + case 'n': negate_omega = 1; break; + case 'F': ok_d.only_f_inrange = 1; break; + case 'a': ok_d.rel_amp_thresh = atof(optarg); break; + case 'A': ok_d.amp_thresh = atof(optarg); break; + case 'E': ok_d.err_thresh = atof(optarg); break; + case 'e': ok_d.rel_err_thresh = atof(optarg); break; + case 'Q': ok_d.Q_thresh = atof(optarg); break; + case 't': dt = atof(optarg); break; + case 'f': nfmin = atoi(optarg); break; + case 'd': + density = atof(optarg); + if (density < 0) { + fprintf(stderr, "harminv: -d argument must be >= 0\n"); return EXIT_FAILURE; - } - - /* harminv requires nf > 1 */ - if (nfmin < 2) nfmin = 2; - - data = read_input_data(stdin, &n, verbose); - - if (n < 1) { - fprintf(stderr, "harminv: no data read\n"); - return EXIT_FAILURE; - } - - if (verbose) - printf("# harminv: %d inputs, dt = %g\n", n, dt); - - printf("frequency, decay constant, Q, amplitude, phase, error\n"); - - ok_d.verbose = verbose; - - for (iarg = optind; iarg < argc; ++iarg) { - double fmin, fmax; - int i; - int *isort = NULL; - - if (sscanf(argv[iarg], "%lf-%lf", &fmin, &fmax) != 2) { - fprintf(stderr, "harminv: invalid argument \"%s\"\n", - argv[iarg]); - return EXIT_FAILURE; - } - if (specify_periods) { - if (fmin == 0 || fmax == 0) { - fprintf(stderr, "harminv: invalid argument \"%s\"" - ": 0 not a valid period\n", argv[iarg]); - return EXIT_FAILURE; - } - fmin = 1/fmin; - fmax = 1/fmax; - } - if (specify_omega) { - fmin /= TWOPI; - fmax /= TWOPI; - } - if (negate_omega) - dt *= -1; - if ((fmin > fmax && dt > 0) || (fmin < fmax && dt < 0)) { - double dummy = fmin; - fmin = fmax; - fmax = dummy; - } - if (verbose) - printf("# searching frequency range %g - %g\n", fmin, fmax); - - ok_d.fmin = fabs(fmin*dt); - ok_d.fmax = fabs(fmax*dt); - if (ok_d.fmin > ok_d.fmax) { - double dummy = ok_d.fmin; - ok_d.fmin = ok_d.fmax; - ok_d.fmax = dummy; - } - - nf = (fmax - fmin) * dt * n * density; - if (nf > NFMAX) nf = NFMAX; - if (nf < nfmin) nf = nfmin; - if (verbose) - printf("# using %d spectral basis functions, density %g\n", - nf, nf / ((fmax - fmin) * dt * n)); - - hd = harminv_data_create(n, data, fmin*dt, fmax*dt, nf); + } + break; + case 's': + switch (tolower(optarg[0])) { + case 'f': sortby = SORT_FREQUENCY; break; + case 'd': sortby = SORT_DECAY; break; + case 'e': sortby = SORT_ERROR; break; + case 'a': sortby = SORT_AMPLITUDE; break; + case 'q': sortby = SORT_Q; break; + default: + fprintf(stderr, "harminv: invalid sort type -s %c\n", tolower(optarg[0])); + usage(stderr); + return EXIT_FAILURE; + } + break; + default: + fprintf(stderr, "harminv: invalid argument -%c\n", c); + usage(stderr); + return EXIT_FAILURE; + } + if (optind == argc) { /* no parameters left */ + fprintf(stderr, "harminv: missing required frequency range(s)\n"); + usage(stderr); + return EXIT_FAILURE; + } + + /* harminv requires nf > 1 */ + if (nfmin < 2) nfmin = 2; + + data = read_input_data(stdin, &n, verbose); + + if (n < 1) { + fprintf(stderr, "harminv: no data read\n"); + return EXIT_FAILURE; + } + + if (verbose) printf("# harminv: %d inputs, dt = %g\n", n, dt); + + printf("frequency, decay constant, Q, amplitude, phase, error\n"); + + ok_d.verbose = verbose; + + for (iarg = optind; iarg < argc; ++iarg) { + double fmin, fmax; + int i; + int *isort = NULL; + + if (sscanf(argv[iarg], "%lf-%lf", &fmin, &fmax) != 2) { + fprintf(stderr, "harminv: invalid argument \"%s\"\n", argv[iarg]); + return EXIT_FAILURE; + } + if (specify_periods) { + if (fmin == 0 || fmax == 0) { + fprintf(stderr, + "harminv: invalid argument \"%s\"" + ": 0 not a valid period\n", + argv[iarg]); + return EXIT_FAILURE; + } + fmin = 1 / fmin; + fmax = 1 / fmax; + } + if (specify_omega) { + fmin /= TWOPI; + fmax /= TWOPI; + } + if (negate_omega) dt *= -1; + if ((fmin > fmax && dt > 0) || (fmin < fmax && dt < 0)) { + double dummy = fmin; + fmin = fmax; + fmax = dummy; + } + if (verbose) printf("# searching frequency range %g - %g\n", fmin, fmax); + + ok_d.fmin = fabs(fmin * dt); + ok_d.fmax = fabs(fmax * dt); + if (ok_d.fmin > ok_d.fmax) { + double dummy = ok_d.fmin; + ok_d.fmin = ok_d.fmax; + ok_d.fmax = dummy; + } + + nf = (fmax - fmin) * dt * n * density; + if (nf > NFMAX) nf = NFMAX; + if (nf < nfmin) nf = nfmin; + if (verbose) + printf("# using %d spectral basis functions, density %g\n", nf, + nf / ((fmax - fmin) * dt * n)); + + hd = harminv_data_create(n, data, fmin * dt, fmax * dt, nf); #if SOLVE_OK_ONLY - harminv_solve_ok_modes(hd, mode_ok, &ok_d); + harminv_solve_ok_modes(hd, mode_ok, &ok_d); #elif SOLVE_ONCE_ONLY - harminv_solve_once(hd); + harminv_solve_once(hd); #else - harminv_solve(hd); + harminv_solve(hd); #endif - mode_ok(hd, -1, &ok_d); /* initialize ok_d */ + mode_ok(hd, -1, &ok_d); /* initialize ok_d */ - CHK_MALLOC(isort, int, harminv_get_num_freqs(hd)); - for (i = 0; i < harminv_get_num_freqs(hd); ++i) - isort[i] = i; - qsort(isort, harminv_get_num_freqs(hd), sizeof(int), compar); + CHK_MALLOC(isort, int, harminv_get_num_freqs(hd)); + for (i = 0; i < harminv_get_num_freqs(hd); ++i) + isort[i] = i; + qsort(isort, harminv_get_num_freqs(hd), sizeof(int), compar); - for (i = 0; i < harminv_get_num_freqs(hd); ++i) { - double freq, decay, err; - cmplx amp; - int j = isort[i]; + for (i = 0; i < harminv_get_num_freqs(hd); ++i) { + double freq, decay, err; + cmplx amp; + int j = isort[i]; #if SOLVE_OK_ONLY - CHECK(mode_ok(hd, j, &ok_d), "bug: invalid mode"); + CHECK(mode_ok(hd, j, &ok_d), "bug: invalid mode"); #else - if (!mode_ok(hd, j, &ok_d)) - continue; + if (!mode_ok(hd, j, &ok_d)) continue; #endif - freq = harminv_get_freq(hd, j) / dt; - decay = harminv_get_decay(hd, j) / fabs(dt); - harminv_get_amplitude(&, hd, j); - err = harminv_get_freq_error(hd, j); - printf("%g, %e, %g, %g, %g, %e\n", - freq * (specify_omega ? TWOPI : 1.0), decay, - harminv_get_Q(hd, j), - cabs(amp), carg(amp) * (dt < 0 ? -1 : 1), err); - } + freq = harminv_get_freq(hd, j) / dt; + decay = harminv_get_decay(hd, j) / fabs(dt); + harminv_get_amplitude(&, hd, j); + err = harminv_get_freq_error(hd, j); + printf("%g, %e, %g, %g, %g, %e\n", freq * (specify_omega ? TWOPI : 1.0), decay, + harminv_get_Q(hd, j), cabs(amp), carg(amp) * (dt < 0 ? -1 : 1), err); + } #if !SOLVE_OK_ONLY - mode_ok(hd, -2, &ok_d); + mode_ok(hd, -2, &ok_d); #endif - harminv_data_destroy(hd); - } + harminv_data_destroy(hd); + } - free(data); - return EXIT_SUCCESS; + free(data); + return EXIT_SUCCESS; } #ifdef F77_DUMMY_MAIN -# ifdef __cplusplus +#ifdef __cplusplus extern "C" -# endif -int F77_DUMMY_MAIN() { return 1; } +#endif + int + F77_DUMMY_MAIN() { + return 1; +} #endif diff --git a/harminv.c b/harminv.c index f115923..1b298e8 100644 --- a/harminv.c +++ b/harminv.c @@ -35,10 +35,7 @@ #undef clog #define clog my_clog -static cmplx my_clog(cmplx z) -{ - return (log(cabs(z)) + I * carg(z)); -} +static cmplx my_clog(cmplx z) { return (log(cabs(z)) + I * carg(z)); } /**************************************************************************/ @@ -64,7 +61,7 @@ static cmplx my_clog(cmplx z) Rongqing Chen and Hua Guo, "Efficient calculation of matrix elements in low storate filter diagonalization," J. Chem. Phys., - vol. 111, no. 2, p. 464-471(Jul. 8 1999). + vol. 111, no. 2, p. 464-471(Jul. 8 1999). The seminal work (though less practical than the M&T algorithm) on this class of methods was done by: @@ -74,7 +71,7 @@ static cmplx my_clog(cmplx z) normal mode frequencies from a small number of residues or a short-time segment of a signal. I. Theory and application to a quantum-dynamics model," J. Chem. Phys., 102, no. 20, p. 8011-8022 - (May 22 1995). + (May 22 1995). A more recent reference is: @@ -93,19 +90,19 @@ static cmplx my_clog(cmplx z) /* Crays have float == double, and don't have the Z* functions in LAPACK/BLAS...we have to use C*. Sigh. */ #if defined(CRAY) || defined(_UNICOS) || defined(_CRAYMPP) -# define BLAS_FUNC(x,X) F77_FUNC(c##x,C##X) -#else /* ! CRAY */ -# define BLAS_FUNC(x,X) F77_FUNC(z##x,Z##X) +#define BLAS_FUNC(x, X) F77_FUNC(c##x, C##X) +#else /* ! CRAY */ +#define BLAS_FUNC(x, X) F77_FUNC(z##x, Z##X) #endif /* ! CRAY */ -#define ZGEEV BLAS_FUNC(geev,GEEV) -#define ZGGEVX BLAS_FUNC(ggevx,GGEVX) -#define ZGGEV BLAS_FUNC(ggev,GGEV) -#define ZGEMM BLAS_FUNC(gemm,GEMM) -#define ZCOPY BLAS_FUNC(copy,COPY) -#define ZAXPY BLAS_FUNC(axpy,AXPY) -#define ZGEMV BLAS_FUNC(gemv,GEMV) -#define ZSCAL BLAS_FUNC(scal,SCAL) +#define ZGEEV BLAS_FUNC(geev, GEEV) +#define ZGGEVX BLAS_FUNC(ggevx, GGEVX) +#define ZGGEV BLAS_FUNC(ggev, GGEV) +#define ZGEMM BLAS_FUNC(gemm, GEMM) +#define ZCOPY BLAS_FUNC(copy, COPY) +#define ZAXPY BLAS_FUNC(axpy, AXPY) +#define ZGEMV BLAS_FUNC(gemv, GEMV) +#define ZSCAL BLAS_FUNC(scal, SCAL) #ifdef __cplusplus extern "C" { @@ -114,37 +111,31 @@ extern "C" { /* We have to pass strings in special ways on Crays, even for passing a single character as with LAPACK. Sigh. */ #if defined(CRAY) || defined(_UNICOS) || defined(_CRAYMPP) -# include -# define FCHARP _fcd -# define F_(s) _cptofcd(s,1) /* second argument is the string length */ -#else /* ! CRAY */ -# define FCHARP const char* -# define F_(s) (s) +#include +#define FCHARP _fcd +#define F_(s) _cptofcd(s, 1) /* second argument is the string length */ +#else /* ! CRAY */ +#define FCHARP const char * +#define F_(s) (s) #endif -extern void ZGEEV(FCHARP,FCHARP, int*, cmplx*,int*, cmplx*, cmplx*,int*, - cmplx*,int*, cmplx*,int*, double*, int*); -extern void ZGGEVX(FCHARP balanc, FCHARP jobvl, FCHARP jobvr, FCHARP sense, - int* n, cmplx* a, int* lda, cmplx* b, int* ldb, - cmplx* alpha, cmplx* beta, - cmplx* vl, int*, cmplx* vr, int*, - int *ilo, int *ihi, double *lscale, double *rscale, - double *abnrm, double *bbnrm, - double *rconde, double *rcondv, - cmplx *wrk, int *lwork, double *rwork, int *iwork, - int *bwork, int *info); -extern void ZGGEV(FCHARP jobvl, FCHARP jobvr, - int* n, cmplx* a, int* lda, cmplx* b, int* ldb, - cmplx* alpha, cmplx* beta, - cmplx* vl, int*, cmplx* vr, int*, - cmplx *wrk, int *lwork, double *rwork, int *info); -extern void ZGEMM(FCHARP,FCHARP, int*,int*,int*, cmplx*, - cmplx*,int*, cmplx*,int*, cmplx*, cmplx*,int*); -extern void ZCOPY(int*, const cmplx*,int*, cmplx*,int*); -extern void ZAXPY(int*, cmplx*, cmplx*,int*, cmplx*,int*); -extern void ZGEMV(FCHARP, int*,int*, cmplx*, cmplx*,int*, cmplx*,int*, - cmplx*, cmplx*,int*); -extern void ZSCAL(int*, cmplx*, cmplx*,int*); +extern void ZGEEV(FCHARP, FCHARP, int *, cmplx *, int *, cmplx *, cmplx *, int *, cmplx *, int *, + cmplx *, int *, double *, int *); +extern void ZGGEVX(FCHARP balanc, FCHARP jobvl, FCHARP jobvr, FCHARP sense, int *n, cmplx *a, + int *lda, cmplx *b, int *ldb, cmplx *alpha, cmplx *beta, cmplx *vl, int *, + cmplx *vr, int *, int *ilo, int *ihi, double *lscale, double *rscale, + double *abnrm, double *bbnrm, double *rconde, double *rcondv, cmplx *wrk, + int *lwork, double *rwork, int *iwork, int *bwork, int *info); +extern void ZGGEV(FCHARP jobvl, FCHARP jobvr, int *n, cmplx *a, int *lda, cmplx *b, int *ldb, + cmplx *alpha, cmplx *beta, cmplx *vl, int *, cmplx *vr, int *, cmplx *wrk, + int *lwork, double *rwork, int *info); +extern void ZGEMM(FCHARP, FCHARP, int *, int *, int *, cmplx *, cmplx *, int *, cmplx *, int *, + cmplx *, cmplx *, int *); +extern void ZCOPY(int *, const cmplx *, int *, cmplx *, int *); +extern void ZAXPY(int *, cmplx *, cmplx *, int *, cmplx *, int *); +extern void ZGEMV(FCHARP, int *, int *, cmplx *, cmplx *, int *, cmplx *, int *, cmplx *, cmplx *, + int *); +extern void ZSCAL(int *, cmplx *, cmplx *, int *); extern void HARMINV_ZDOTU(cmplx *, int *, cmplx *, int *, cmplx *, int *); #ifdef __cplusplus @@ -154,22 +145,19 @@ extern void HARMINV_ZDOTU(cmplx *, int *, cmplx *, int *, cmplx *, int *); /**************************************************************************/ /* compute c^n, where n is an integer: */ -static cmplx cpow_i(cmplx c, int n) -{ - if (n < 0) - return (1.0 / cpow_i(c, -n)); - else { - cmplx result = 1; - while (n > 1) { - if (n % 2 == 1) - result *= c; - c *= c; - n /= 2; - } - if (n > 0) - result *= c; - return result; - } +static cmplx cpow_i(cmplx c, int n) { + if (n < 0) + return (1.0 / cpow_i(c, -n)); + else { + cmplx result = 1; + while (n > 1) { + if (n % 2 == 1) result *= c; + c *= c; + n /= 2; + } + if (n > 0) result *= c; + return result; + } } /* Computing powers by cumulative multiplication, below, is faster @@ -181,7 +169,7 @@ static cmplx cpow_i(cmplx c, int n) /* FIXME: instead of this, we should really do the first-order expansion of the U matrix in |z - z2|, below. */ #define SMALL (1e-12) -#define C_CLOSE(c1,c2) (cabs((c1) - (c2)) < SMALL) +#define C_CLOSE(c1, c2) (cabs((c1) - (c2)) < SMALL) /**************************************************************************/ @@ -190,264 +178,249 @@ static cmplx cpow_i(cmplx c, int n) must be the case that no two elements of z are the same and that J2 == J1; in this case the matrix U will be symmetric. Note that c must be an array whose size n, is at least 2*K+p elements. */ -static void generate_U(cmplx *U, cmplx *U1, - int p, - const cmplx *c, int n, - int K, - int J, int J2, const cmplx *z, const cmplx *z2, - cmplxl **G0, cmplxl **G0_M, cmplxl **D0) -{ - int M = K - 1; - int i, j, m; - /* temp. arrays for 1/z, z^(-m), z^(-M), the G function of C&G, - and the diagonal elements D[i] = U(z[i],z[i]): */ - cmplx *z_inv, *z_m, *z_M; - cmplxl *G, *G_M, *D; - cmplx *z2_inv, *z2_m, *z2_M; - cmplxl *G2, *G2_M; - - CHECK(U && c && z && z2, "invalid arguments to generate_U"); - CHECK(n >= 2*K + p, "too few coefficients in generate_U"); - CHECK(z != z2 || J == J2, "invalid sizes passed to generate_U"); - CHECK((!G0 && !G0_M && !D0) || (G0 && G0_M && D0), - "G0/G0_M/D0 must be all-non-NULL/all-NULL"); - CHECK(!G0 || (!*G0 && !*G0_M && !*D0) || (*G0 && *G0_M && *D0), - "*G0/*G0_M/*D0 must be all-non-NULL/all-NULL"); - - /* Now, compute U according to eqs. 25-27 of Chen & Guo, but - using the notation of eq. 25 of M&T. This operation has - complexity O(N*J + J*J). At the same time, we can compute the - matrix U1 as well by eqs. 29-30 of C&G, saving an extra pass - over the input array. */ - - /* first, set up some temporary arrays for caching things like - z^m and 1/z, so we don't need to recompute them all the time. */ - CHK_MALLOC(z_inv, cmplx, J); - CHK_MALLOC(z_m, cmplx, J); - CHK_MALLOC(z_M, cmplx, J); - if (G0 && *G0) { - G = *G0; - G_M = *G0_M; - D = *D0; - } - else { - CHK_MALLOC(G, cmplxl, J); - CHK_MALLOC(G_M, cmplxl, J); - CHK_MALLOC(D, cmplxl, J); - for (i = 0; i < J; ++i) { - D[i] = G[i] = G_M[i] = 0; - } - } - for (i = 0; i < J; ++i) { - z_inv[i] = 1.0 / z[i]; - z_m[i] = 1; - z_M[i] = cpow_i(z[i], -M); - } - if (z2 != z) { - CHK_MALLOC(z2_inv, cmplx, J2); - CHK_MALLOC(z2_m, cmplx, J2); - CHK_MALLOC(z2_M, cmplx, J2); - CHK_MALLOC(G2, cmplxl, J2); - CHK_MALLOC(G2_M, cmplxl, J2); - for (i = 0; i < J2; ++i) { - z2_inv[i] = 1.0 / z2[i]; - z2_m[i] = 1; - z2_M[i] = cpow_i(z2[i], -M); - G2[i] = G2_M[i] = 0; - } - } - else { - z2_inv = z2_m = z2_M = NULL; - G2 = G2_M = NULL; - } - - /* First, loop over the signal array (c), building up the - spectral functions G and G_M (corresponding to G_p and - G_{p+M+1} in C&G), as well as the diagonal matrix entries: */ - for (m = 0; m <= M; ++m) { - cmplx c1 = c[m + p], c2 = c[m + p + M + 1]; - double d = m + 1; /* M - fabs(M - m) + 1 */ - double d2 = M - m; /* M - fabs(M - (m + M + 1)) + 1 */ - - if (!G0 || !*G0) { - for (i = 0; i < J; ++i) { - cmplx x1 = z_m[i] * c1; - cmplx x2 = z_m[i] * c2; - G[i] += x1; - G_M[i] += x2; - D[i] += x1 * d + x2 * d2 * z_M[i] * z_inv[i]; - if (m % NPOW == NPOW - 1) - z_m[i] = cpow_i(z_inv[i], m + 1); - else - z_m[i] *= z_inv[i]; - } - } - if (z2 != z) - for (i = 0; i < J2; ++i) { - G2[i] += z2_m[i] * c1; - G2_M[i] += z2_m[i] * c2; - if (m % NPOW == NPOW - 1) - z2_m[i] = cpow_i(z2_inv[i], m + 1); - else - z2_m[i] *= z2_inv[i]; - } - } - - /* Compute U (or just the upper part if U is symmetric), via the - formula from C&G; compute U1 at the same time as in C&G. */ - if (z2 != z) { - for (i = 0; i < J; ++i) - for (j = 0; j < J2; ++j) { - if (C_CLOSE(z[i], z2[j])) - U[i*J2 + j] = D[i]; - else - U[i*J2 + j] = (z[i] * G2[j] - z2[j] * G[i] + - z2_M[j] * G_M[i] - z_M[i] * G2_M[j]) - / (z[i] - z2[j]); - } - - if (U1) - for (i = 0; i < J; ++i) - for (j = 0; j < J2; ++j) { - if (C_CLOSE(z[i], z2[j])) - U1[i*J2 + j] = z[i] * (D[i] - G[i]) + - z_M[i] * G_M[i]; - else - U1[i*J2 + j] = (z[i] * z2[j] * (G2[j] - G[i]) - + z2_M[j] * z[i] * G_M[i] - - z_M[i] * z2[j] * G2_M[j]) - / (z[i] - z2[j]); - } - } - else { /* z == z2 */ - for (i = 0; i < J; ++i) { - U[i*J + i] = D[i]; - for (j = i + 1; j < J; ++j) { - U[i*J + j] = (z[i] * G[j] - z[j] * G[i] + - z_M[j] * G_M[i] - z_M[i] * G_M[j]) - / (z[i] - z[j]); - } - } - - if (U1) - for (i = 0; i < J; ++i) { - U1[i*J + i] = z[i] * (D[i] - G[i]) + z_M[i] * G_M[i]; - for (j = i + 1; j < J; ++j) { - U1[i*J + j] = (z[i] * z[j] * (G[j] - G[i]) - + z_M[j] * z[i] * G_M[i] - - z_M[i] * z[j] * G_M[j]) - / (z[i] - z[j]); - } - } - } - - - /* finally, copy the upper to the lower triangle if U is symmetric: */ - if (z == z2) { - for (i = 0; i < J; ++i) - for (j = i + 1; j < J; ++j) - U[j*J + i] = U[i*J + j]; - if (U1) - for (i = 0; i < J; ++i) - for (j = i + 1; j < J; ++j) - U1[j*J + i] = U1[i*J + j]; - } - - free(G2_M); - free(G2); - free(z2_M); - free(z2_m); - free(z2_inv); - - if (G0 && !*G0) { - *G0 = G; - *G0_M = G_M; - *D0 = D; - } - else if (!G0) { - free(D); - free(G_M); - free(G); - } - free(z_M); - free(z_m); - free(z_inv); +static void generate_U(cmplx *U, cmplx *U1, int p, const cmplx *c, int n, int K, int J, int J2, + const cmplx *z, const cmplx *z2, cmplxl **G0, cmplxl **G0_M, cmplxl **D0) { + int M = K - 1; + int i, j, m; + /* temp. arrays for 1/z, z^(-m), z^(-M), the G function of C&G, + and the diagonal elements D[i] = U(z[i],z[i]): */ + cmplx *z_inv, *z_m, *z_M; + cmplxl *G, *G_M, *D; + cmplx *z2_inv, *z2_m, *z2_M; + cmplxl *G2, *G2_M; + + CHECK(U && c && z && z2, "invalid arguments to generate_U"); + CHECK(n >= 2 * K + p, "too few coefficients in generate_U"); + CHECK(z != z2 || J == J2, "invalid sizes passed to generate_U"); + CHECK((!G0 && !G0_M && !D0) || (G0 && G0_M && D0), "G0/G0_M/D0 must be all-non-NULL/all-NULL"); + CHECK(!G0 || (!*G0 && !*G0_M && !*D0) || (*G0 && *G0_M && *D0), + "*G0/*G0_M/*D0 must be all-non-NULL/all-NULL"); + + /* Now, compute U according to eqs. 25-27 of Chen & Guo, but + using the notation of eq. 25 of M&T. This operation has + complexity O(N*J + J*J). At the same time, we can compute the + matrix U1 as well by eqs. 29-30 of C&G, saving an extra pass + over the input array. */ + + /* first, set up some temporary arrays for caching things like + z^m and 1/z, so we don't need to recompute them all the time. */ + CHK_MALLOC(z_inv, cmplx, J); + CHK_MALLOC(z_m, cmplx, J); + CHK_MALLOC(z_M, cmplx, J); + if (G0 && *G0) { + G = *G0; + G_M = *G0_M; + D = *D0; + } + else { + CHK_MALLOC(G, cmplxl, J); + CHK_MALLOC(G_M, cmplxl, J); + CHK_MALLOC(D, cmplxl, J); + for (i = 0; i < J; ++i) { + D[i] = G[i] = G_M[i] = 0; + } + } + for (i = 0; i < J; ++i) { + z_inv[i] = 1.0 / z[i]; + z_m[i] = 1; + z_M[i] = cpow_i(z[i], -M); + } + if (z2 != z) { + CHK_MALLOC(z2_inv, cmplx, J2); + CHK_MALLOC(z2_m, cmplx, J2); + CHK_MALLOC(z2_M, cmplx, J2); + CHK_MALLOC(G2, cmplxl, J2); + CHK_MALLOC(G2_M, cmplxl, J2); + for (i = 0; i < J2; ++i) { + z2_inv[i] = 1.0 / z2[i]; + z2_m[i] = 1; + z2_M[i] = cpow_i(z2[i], -M); + G2[i] = G2_M[i] = 0; + } + } + else { + z2_inv = z2_m = z2_M = NULL; + G2 = G2_M = NULL; + } + + /* First, loop over the signal array (c), building up the + spectral functions G and G_M (corresponding to G_p and + G_{p+M+1} in C&G), as well as the diagonal matrix entries: */ + for (m = 0; m <= M; ++m) { + cmplx c1 = c[m + p], c2 = c[m + p + M + 1]; + double d = m + 1; /* M - fabs(M - m) + 1 */ + double d2 = M - m; /* M - fabs(M - (m + M + 1)) + 1 */ + + if (!G0 || !*G0) { + for (i = 0; i < J; ++i) { + cmplx x1 = z_m[i] * c1; + cmplx x2 = z_m[i] * c2; + G[i] += x1; + G_M[i] += x2; + D[i] += x1 * d + x2 * d2 * z_M[i] * z_inv[i]; + if (m % NPOW == NPOW - 1) + z_m[i] = cpow_i(z_inv[i], m + 1); + else + z_m[i] *= z_inv[i]; + } + } + if (z2 != z) + for (i = 0; i < J2; ++i) { + G2[i] += z2_m[i] * c1; + G2_M[i] += z2_m[i] * c2; + if (m % NPOW == NPOW - 1) + z2_m[i] = cpow_i(z2_inv[i], m + 1); + else + z2_m[i] *= z2_inv[i]; + } + } + + /* Compute U (or just the upper part if U is symmetric), via the + formula from C&G; compute U1 at the same time as in C&G. */ + if (z2 != z) { + for (i = 0; i < J; ++i) + for (j = 0; j < J2; ++j) { + if (C_CLOSE(z[i], z2[j])) + U[i * J2 + j] = D[i]; + else + U[i * J2 + j] = + (z[i] * G2[j] - z2[j] * G[i] + z2_M[j] * G_M[i] - z_M[i] * G2_M[j]) / (z[i] - z2[j]); + } + + if (U1) + for (i = 0; i < J; ++i) + for (j = 0; j < J2; ++j) { + if (C_CLOSE(z[i], z2[j])) + U1[i * J2 + j] = z[i] * (D[i] - G[i]) + z_M[i] * G_M[i]; + else + U1[i * J2 + j] = (z[i] * z2[j] * (G2[j] - G[i]) + z2_M[j] * z[i] * G_M[i] - + z_M[i] * z2[j] * G2_M[j]) / + (z[i] - z2[j]); + } + } + else { /* z == z2 */ + for (i = 0; i < J; ++i) { + U[i * J + i] = D[i]; + for (j = i + 1; j < J; ++j) { + U[i * J + j] = + (z[i] * G[j] - z[j] * G[i] + z_M[j] * G_M[i] - z_M[i] * G_M[j]) / (z[i] - z[j]); + } + } + + if (U1) + for (i = 0; i < J; ++i) { + U1[i * J + i] = z[i] * (D[i] - G[i]) + z_M[i] * G_M[i]; + for (j = i + 1; j < J; ++j) { + U1[i * J + j] = + (z[i] * z[j] * (G[j] - G[i]) + z_M[j] * z[i] * G_M[i] - z_M[i] * z[j] * G_M[j]) / + (z[i] - z[j]); + } + } + } + + /* finally, copy the upper to the lower triangle if U is symmetric: */ + if (z == z2) { + for (i = 0; i < J; ++i) + for (j = i + 1; j < J; ++j) + U[j * J + i] = U[i * J + j]; + if (U1) + for (i = 0; i < J; ++i) + for (j = i + 1; j < J; ++j) + U1[j * J + i] = U1[i * J + j]; + } + + free(G2_M); + free(G2); + free(z2_M); + free(z2_m); + free(z2_inv); + + if (G0 && !*G0) { + *G0 = G; + *G0_M = G_M; + *D0 = D; + } + else if (!G0) { + free(D); + free(G_M); + free(G); + } + free(z_M); + free(z_m); + free(z_inv); } /**************************************************************************/ -static void init_z(harminv_data d, int J, cmplx *z) -{ - d->J = J; - d->z = z; - CHK_MALLOC(d->U0, cmplx, J*J); - CHK_MALLOC(d->U1, cmplx, J*J); - generate_U(d->U0, d->U1, 0, d->c, d->n, d->K, d->J, d->J, d->z, d->z, - &d->G0, &d->G0_M, &d->D0); +static void init_z(harminv_data d, int J, cmplx *z) { + d->J = J; + d->z = z; + CHK_MALLOC(d->U0, cmplx, J * J); + CHK_MALLOC(d->U1, cmplx, J * J); + generate_U(d->U0, d->U1, 0, d->c, d->n, d->K, d->J, d->J, d->z, d->z, &d->G0, &d->G0_M, &d->D0); } /**************************************************************************/ -harminv_data harminv_data_create(int n, - const cmplx *signal, - double fmin, double fmax, int nf) -{ - int i; - harminv_data d; - - CHECK(nf > 1, "# frequencies must > 1"); - CHECK(n > 0, "invalid number of data points"); - CHECK(signal, "invalid NULL signal array"); - CHECK(fmin < fmax, "should have fmin < fmax"); - - CHK_MALLOC(d, struct harminv_data_struct, 1); - d->c = signal; - d->n = n; - d->K = n/2 - 1; - d->fmin = fmin; - d->fmax = fmax; - d->nfreqs = -1; /* we haven't computed eigen-solutions yet */ - d->B = d->u = d->amps = d->U0 = d->U1 = (cmplx *) NULL; - d->G0 = d->G0_M = d->D0 = (cmplxl *) NULL; - d->errs = (double *) NULL; - - CHK_MALLOC(d->z, cmplx, nf); - for (i = 0; i < nf; ++i) - d->z[i] = cexp(-I * TWOPI * (fmin + i * ((fmax - fmin) / (nf - 1)))); - - init_z(d, nf, d->z); - - return d; +harminv_data harminv_data_create(int n, const cmplx *signal, double fmin, double fmax, int nf) { + int i; + harminv_data d; + + CHECK(nf > 1, "# frequencies must > 1"); + CHECK(n > 0, "invalid number of data points"); + CHECK(signal, "invalid NULL signal array"); + CHECK(fmin < fmax, "should have fmin < fmax"); + + CHK_MALLOC(d, struct harminv_data_struct, 1); + d->c = signal; + d->n = n; + d->K = n / 2 - 1; + d->fmin = fmin; + d->fmax = fmax; + d->nfreqs = -1; /* we haven't computed eigen-solutions yet */ + d->B = d->u = d->amps = d->U0 = d->U1 = (cmplx *)NULL; + d->G0 = d->G0_M = d->D0 = (cmplxl *)NULL; + d->errs = (double *)NULL; + + CHK_MALLOC(d->z, cmplx, nf); + for (i = 0; i < nf; ++i) + d->z[i] = cexp(-I * TWOPI * (fmin + i * ((fmax - fmin) / (nf - 1)))); + + init_z(d, nf, d->z); + + return d; } /**************************************************************************/ -void harminv_data_destroy(harminv_data d) -{ - if (d) { - free(d->u); free(d->B); - free(d->U1); free(d->U0); - free(d->G0); free(d->G0_M); free(d->D0); - free(d->z); - free(d->amps); - free(d->errs); - free(d); - } +void harminv_data_destroy(harminv_data d) { + if (d) { + free(d->u); + free(d->B); + free(d->U1); + free(d->U0); + free(d->G0); + free(d->G0_M); + free(d->D0); + free(d->z); + free(d->amps); + free(d->errs); + free(d); + } } /**************************************************************************/ /* Compute the symmetric dot product of x and y, both vectors of - length n. If they are column-vectors, this is: transpose(x) * y. + length n. If they are column-vectors, this is: transpose(x) * y. (We could use the BLAS ZDOTU function for this, but calling Fortran functions, as opposed to subroutines, from C is problematic.) */ -static cmplx symmetric_dot(int n, cmplx *x, cmplx *y) -{ - cmplxl dot = 0; - int i; - for (i = 0; i < n; ++i) - dot += x[i] * y[i]; - return dot; +static cmplx symmetric_dot(int n, cmplx *x, cmplx *y) { + cmplxl dot = 0; + int i; + for (i = 0; i < n; ++i) + dot += x[i] * y[i]; + return dot; } /**************************************************************************/ @@ -458,27 +431,26 @@ static cmplx symmetric_dot(int n, cmplx *x, cmplx *y) complex-symmetric n x n matrix A. The eigenvectors are normalized to 1 according to the symmetric dot product (i.e. no complex conjugation). */ -static void solve_eigenvects(int n, const cmplx *A0, cmplx *V, cmplx *v) -{ - int lwork, info; - cmplx *work; - double *rwork; - cmplx *A; - - /* according to the ZGEEV documentation, the matrix A is overwritten, - and we don't want to overwrite our input A0 */ - CHK_MALLOC(A, cmplx, n*n); - { - int n2 = n*n, one = 1; - ZCOPY(&n2, A0, &one, A, &one); - } - - /* Unfortunately, LAPACK doesn't have a special solver for the - complex-symmetric eigenproblem. For now, just use the general - non-symmetric solver, and realize that the left eigenvectors - are the complex-conjugates of the right eigenvectors. */ - -#if 0 /* LAPACK seems to be buggy here, returning ridiculous sizes at times */ +static void solve_eigenvects(int n, const cmplx *A0, cmplx *V, cmplx *v) { + int lwork, info; + cmplx *work; + double *rwork; + cmplx *A; + + /* according to the ZGEEV documentation, the matrix A is overwritten, + and we don't want to overwrite our input A0 */ + CHK_MALLOC(A, cmplx, n * n); + { + int n2 = n * n, one = 1; + ZCOPY(&n2, A0, &one, A, &one); + } + + /* Unfortunately, LAPACK doesn't have a special solver for the + complex-symmetric eigenproblem. For now, just use the general + non-symmetric solver, and realize that the left eigenvectors + are the complex-conjugates of the right eigenvectors. */ + +#if 0 /* LAPACK seems to be buggy here, returning ridiculous sizes at times */ cmplx wsize; lwork = -1; /* compute optimal workspace size */ ZGEEV(F_("N"), F_("V"), &n, A, &n, v, V, &n, V, &n, &wsize, &lwork, rwork, &info); @@ -488,32 +460,32 @@ static void solve_eigenvects(int n, const cmplx *A0, cmplx *V, cmplx *v) lwork = 2*n; CHECK(lwork > 0, "zgeev is not returning a positive work size!"); #else - lwork = 4*n; /* minimum is 2*n; we'll be generous. */ + lwork = 4 * n; /* minimum is 2*n; we'll be generous. */ #endif - CHK_MALLOC(rwork, double, 2*n); - CHK_MALLOC(work, cmplx, lwork); - - ZGEEV(F_("N"), F_("V"), &n, A, &n, v, V, &n, V, &n, work, &lwork, rwork, &info); - - free(work); - free(rwork); - free(A); - - CHECK(info >= 0, "invalid argument to ZGEEV"); - CHECK(info <= 0, "failed convergence in ZGEEV"); - - /* Finally, we need to fix the normalization of the eigenvectors, - since LAPACK normalizes them under the ordinary dot product, - i.e. with complex conjugation. (In principle, do we also need - to re-orthogonalize, for the case of degenerate eigenvalues?) */ - { - int i, one = 1; - for (i = 0; i < n; ++i) { - cmplx norm = 1.0 / csqrt(symmetric_dot(n, V+i*n, V+i*n)); - ZSCAL(&n, &norm, V+i*n, &one); - } - } + CHK_MALLOC(rwork, double, 2 * n); + CHK_MALLOC(work, cmplx, lwork); + + ZGEEV(F_("N"), F_("V"), &n, A, &n, v, V, &n, V, &n, work, &lwork, rwork, &info); + + free(work); + free(rwork); + free(A); + + CHECK(info >= 0, "invalid argument to ZGEEV"); + CHECK(info <= 0, "failed convergence in ZGEEV"); + + /* Finally, we need to fix the normalization of the eigenvectors, + since LAPACK normalizes them under the ordinary dot product, + i.e. with complex conjugation. (In principle, do we also need + to re-orthogonalize, for the case of degenerate eigenvalues?) */ + { + int i, one = 1; + for (i = 0; i < n; ++i) { + cmplx norm = 1.0 / csqrt(symmetric_dot(n, V + i * n, V + i * n)); + ZSCAL(&n, &norm, V + i * n, &one); + } + } } /**************************************************************************/ @@ -524,135 +496,134 @@ static void solve_eigenvects(int n, const cmplx *A0, cmplx *V, cmplx *v) /* Solve the eigenvalue problem U1 b = u U0 b, where b is the eigenvector and u is the eigenvalue. u = exp(iwt - at) then contains both the frequency and the decay constant. */ -void harminv_solve_once(harminv_data d) -{ - int J, i, one=1; - cmplx zone = 1.0, zzero = 0.0; - cmplx *V0, *v0, *H1, *V1; /* for eigensolutions of U0 and U1 */ - double max_v0 = 0.0; - - J = d->J; - CHK_MALLOC(V0, cmplx, J*J); - CHK_MALLOC(v0, cmplx, J); - - /* Unfortunately, U0 is very likely to be singular, so we must - first extract the non-singular eigenvectors and only work in - that sub-space. See the Wall & Neuhauser paper. */ - - solve_eigenvects(J, d->U0, V0, v0); - - /* find maximum |eigenvalue| */ - for (i = 0; i < J; ++i) { - double v = cabs(v0[i]); - if (v > max_v0) - max_v0 = v; - } - - /* we must remove the singular components of U0, those - that are less than some threshold times the maximum eigenvalue. - Also, we need to scale the eigenvectors by 1/sqrt(eigenval). */ - d->nfreqs = J; - for (i = 0; i < J; ++i) { - if (cabs(v0[i]) < SINGULAR_THRESHOLD * max_v0) { - v0[i] = 0; /* tag as a "hole" */ - d->nfreqs -= 1; - } - else { /* not singular */ - cmplx s; - int j; - /* move the eigenvector to the first "hole" left by - deleting singular eigenvalues: */ - for (j = 0; j < i && v0[j] != 0.0; ++j) - ; - if (j < i) { - ZCOPY(&J, V0 + i*J, &one, V0 + j*J, &one); - v0[j] = v0[i]; - v0[i] = 0; /* tag as a "hole" */ - } - s = 1.0 / csqrt(v0[j]); - ZSCAL(&J, &s, V0 + j*J, &one); - } - } - - CHK_MALLOC(d->B, cmplx, d->nfreqs * J); - CHK_MALLOC(d->u, cmplx, d->nfreqs); - CHK_MALLOC(V1, cmplx, d->nfreqs * d->nfreqs); - CHK_MALLOC(H1, cmplx, d->nfreqs * d->nfreqs); - - /* compute H1 = V0 * U1 * V0': */ - - /* B = V0 * U1: */ - ZGEMM(F_("N"), F_("N"), &J, &d->nfreqs, &J, - &zone, d->U1, &J, V0, &J, &zzero, d->B, &J); - /* H1 = B * transpose(V0) */ - ZGEMM(F_("T"), F_("N"), &d->nfreqs, &d->nfreqs, &J, - &zone, V0, &J, d->B, &J, &zzero, H1, &d->nfreqs); - - /* Finally, we can find the eigenvalues and eigenvectors: */ - solve_eigenvects(d->nfreqs, H1, V1, d->u); - /* B = V1 * V0: */ - ZGEMM(F_("N"), F_("N"), &J, &d->nfreqs, &d->nfreqs, - &zone, V0, &J, V1, &d->nfreqs, &zzero, d->B, &J); - - free(H1); - free(V1); - free(v0); - free(V0); +void harminv_solve_once(harminv_data d) { + int J, i, one = 1; + cmplx zone = 1.0, zzero = 0.0; + cmplx *V0, *v0, *H1, *V1; /* for eigensolutions of U0 and U1 */ + double max_v0 = 0.0; + + J = d->J; + CHK_MALLOC(V0, cmplx, J * J); + CHK_MALLOC(v0, cmplx, J); + + /* Unfortunately, U0 is very likely to be singular, so we must + first extract the non-singular eigenvectors and only work in + that sub-space. See the Wall & Neuhauser paper. */ + + solve_eigenvects(J, d->U0, V0, v0); + + /* find maximum |eigenvalue| */ + for (i = 0; i < J; ++i) { + double v = cabs(v0[i]); + if (v > max_v0) max_v0 = v; + } + + /* we must remove the singular components of U0, those + that are less than some threshold times the maximum eigenvalue. + Also, we need to scale the eigenvectors by 1/sqrt(eigenval). */ + d->nfreqs = J; + for (i = 0; i < J; ++i) { + if (cabs(v0[i]) < SINGULAR_THRESHOLD * max_v0) { + v0[i] = 0; /* tag as a "hole" */ + d->nfreqs -= 1; + } + else { /* not singular */ + cmplx s; + int j; + /* move the eigenvector to the first "hole" left by + deleting singular eigenvalues: */ + for (j = 0; j < i && v0[j] != 0.0; ++j) + ; + if (j < i) { + ZCOPY(&J, V0 + i * J, &one, V0 + j * J, &one); + v0[j] = v0[i]; + v0[i] = 0; /* tag as a "hole" */ + } + s = 1.0 / csqrt(v0[j]); + ZSCAL(&J, &s, V0 + j * J, &one); + } + } + + CHK_MALLOC(d->B, cmplx, d->nfreqs * J); + CHK_MALLOC(d->u, cmplx, d->nfreqs); + CHK_MALLOC(V1, cmplx, d->nfreqs * d->nfreqs); + CHK_MALLOC(H1, cmplx, d->nfreqs * d->nfreqs); + + /* compute H1 = V0 * U1 * V0': */ + + /* B = V0 * U1: */ + ZGEMM(F_("N"), F_("N"), &J, &d->nfreqs, &J, &zone, d->U1, &J, V0, &J, &zzero, d->B, &J); + /* H1 = B * transpose(V0) */ + ZGEMM(F_("T"), F_("N"), &d->nfreqs, &d->nfreqs, &J, &zone, V0, &J, d->B, &J, &zzero, H1, + &d->nfreqs); + + /* Finally, we can find the eigenvalues and eigenvectors: */ + solve_eigenvects(d->nfreqs, H1, V1, d->u); + /* B = V1 * V0: */ + ZGEMM(F_("N"), F_("N"), &J, &d->nfreqs, &d->nfreqs, &zone, V0, &J, V1, &d->nfreqs, &zzero, d->B, + &J); + + free(H1); + free(V1); + free(v0); + free(V0); } /**************************************************************************/ /* After solving once, solve again using the solutions from last - time as the input to the spectra estimator this time. + time as the input to the spectra estimator this time. Optionally, if mode_ok is not NULL, we can only use the solutions k from last time that pass ok(d, k, ok_d). Currently, this is not recommended as it seems to make things worse. */ -void harminv_solve_again(harminv_data d, harminv_mode_ok_func ok, void *ok_d) -{ - int i, j; - char *mode_ok = 0; - CHECK(d->nfreqs >= 0, "haven't computed eigensolutions yet"); - - if (!d->nfreqs) return; /* no eigensolutions to work with */ - - if (ok) { - CHK_MALLOC(mode_ok, char, d->nfreqs); - ok(d, -1, ok_d); /* initialize */ - for (i = 0; i < d->nfreqs; ++i) - mode_ok[i] = ok(d, i, ok_d); - } - - free(d->B); - free(d->U1); free(d->U0); free(d->G0); free(d->G0_M); free(d->D0); - free(d->z); - free(d->amps); - free(d->errs); - d->B = d->U1 = d->U0 = d->z = d->amps = (cmplx *) NULL; - d->G0 = d->G0_M = d->D0 = (cmplxl *) NULL; - d->errs = (double *) NULL; - - /* Spectral grid needs to be on the unit circle or system is unstable: */ - for (i = j = 0; i < d->nfreqs; ++i) - if (!ok || mode_ok[i]) - d->u[j++] = d->u[i] / cabs(d->u[i]); - d->nfreqs = j; - - if (ok) { - ok(d, -2, ok_d); /* finish */ - free(mode_ok); - } - - d->u = (cmplx *) realloc(d->u, sizeof(cmplx) * d->nfreqs); - - if (!d->nfreqs) return; /* no eigensolutions to work with */ - - init_z(d, d->nfreqs, d->u); - - d->nfreqs = 0; - d->B = d->u = NULL; - - harminv_solve_once(d); +void harminv_solve_again(harminv_data d, harminv_mode_ok_func ok, void *ok_d) { + int i, j; + char *mode_ok = 0; + CHECK(d->nfreqs >= 0, "haven't computed eigensolutions yet"); + + if (!d->nfreqs) return; /* no eigensolutions to work with */ + + if (ok) { + CHK_MALLOC(mode_ok, char, d->nfreqs); + ok(d, -1, ok_d); /* initialize */ + for (i = 0; i < d->nfreqs; ++i) + mode_ok[i] = ok(d, i, ok_d); + } + + free(d->B); + free(d->U1); + free(d->U0); + free(d->G0); + free(d->G0_M); + free(d->D0); + free(d->z); + free(d->amps); + free(d->errs); + d->B = d->U1 = d->U0 = d->z = d->amps = (cmplx *)NULL; + d->G0 = d->G0_M = d->D0 = (cmplxl *)NULL; + d->errs = (double *)NULL; + + /* Spectral grid needs to be on the unit circle or system is unstable: */ + for (i = j = 0; i < d->nfreqs; ++i) + if (!ok || mode_ok[i]) d->u[j++] = d->u[i] / cabs(d->u[i]); + d->nfreqs = j; + + if (ok) { + ok(d, -2, ok_d); /* finish */ + free(mode_ok); + } + + d->u = (cmplx *)realloc(d->u, sizeof(cmplx) * d->nfreqs); + + if (!d->nfreqs) return; /* no eigensolutions to work with */ + + init_z(d, d->nfreqs, d->u); + + d->nfreqs = 0; + d->B = d->u = NULL; + + harminv_solve_once(d); } /**************************************************************************/ @@ -661,42 +632,37 @@ void harminv_solve_again(harminv_data d, harminv_mode_ok_func ok, void *ok_d) Currently, it is recommended that you use harminv_solve (i.e. pass ok = 0); see harminv_solve_again. */ -void harminv_solve_ok_modes(harminv_data d, harminv_mode_ok_func ok,void *ok_d) -{ - - harminv_solve_once(d); - - /* This is not in the papers, but seems to be a good idea: - plug the u's back in as z's for another pass, and repeat - as long as the number of eigenvalues decreases. Effectively, - this gives us more basis functions where the modes are. */ - { - int prev_nf, cur_nf, nf_ok; - cur_nf = harminv_get_num_freqs(d); - do { - prev_nf = cur_nf; - harminv_solve_again(d, ok, ok_d); - cur_nf = harminv_get_num_freqs(d); - if (ok) { - ok(d, -1, ok_d); /* initialize */ - for (nf_ok = 0; nf_ok < cur_nf - && ok(d, nf_ok, ok_d); ++nf_ok) - ; - ok(d, -2, ok_d); /* finish */ - } - else - nf_ok = cur_nf; - } while (cur_nf < prev_nf || nf_ok < cur_nf); - /* FIXME: solve one more time for good measure? */ - } +void harminv_solve_ok_modes(harminv_data d, harminv_mode_ok_func ok, void *ok_d) { + + harminv_solve_once(d); + + /* This is not in the papers, but seems to be a good idea: + plug the u's back in as z's for another pass, and repeat + as long as the number of eigenvalues decreases. Effectively, + this gives us more basis functions where the modes are. */ + { + int prev_nf, cur_nf, nf_ok; + cur_nf = harminv_get_num_freqs(d); + do { + prev_nf = cur_nf; + harminv_solve_again(d, ok, ok_d); + cur_nf = harminv_get_num_freqs(d); + if (ok) { + ok(d, -1, ok_d); /* initialize */ + for (nf_ok = 0; nf_ok < cur_nf && ok(d, nf_ok, ok_d); ++nf_ok) + ; + ok(d, -2, ok_d); /* finish */ + } + else + nf_ok = cur_nf; + } while (cur_nf < prev_nf || nf_ok < cur_nf); + /* FIXME: solve one more time for good measure? */ + } } /**************************************************************************/ -void harminv_solve(harminv_data d) -{ - harminv_solve_ok_modes(d, NULL, NULL); -} +void harminv_solve(harminv_data d) { harminv_solve_ok_modes(d, NULL, NULL); } /**************************************************************************/ @@ -705,47 +671,42 @@ void harminv_solve(harminv_data d) /* Returns an array (of size harminv_get_num_freqs(d)) of estimates for the |error| in the solution frequencies. Solutions with errors much larger than the smallest error are likely to be spurious. */ -double *harminv_compute_frequency_errors(harminv_data d) -{ - int i, J2, one = 1; - cmplx *U2, *U2b; - double *freq_err; - - CHECK(d->nfreqs >= 0, "haven't computed eigensolutions yet"); - if (!d->nfreqs) return NULL; - - CHK_MALLOC(freq_err, double, d->nfreqs); - - J2 = d->J*d->J; - CHK_MALLOC(U2, cmplx, J2); - generate_U(U2, NULL, 2, d->c, d->n, d->K, d->J, d->J, d->z, d->z, - NULL, NULL, NULL); - CHK_MALLOC(U2b, cmplx, d->J); - - /* For each eigenstate, compute an estimate of the error, roughly - as suggested in W&N, eq. (2.19). */ - - for (i = 0; i < d->nfreqs; ++i) { - cmplx zone = 1.0, zzero = 0.0; - - /* compute U2b = U2 * B[i] */ - ZGEMV(F_("T"), &d->J, &d->J, - &zone, U2, &d->J, d->B + i * d->J, &one, - &zzero, U2b, &one); - - /* ideally, B[i] should satisfy U2 B[i] = u^2 U0 B[i]. - since B U0 B = 1, then we can get a second estimate - for u by sqrt(B[i] U2 B[i]), and from this we compute - the relative error in the (complex) frequency. */ - - freq_err[i] = - cabs(clog(csqrt(symmetric_dot(d->J, d->B + i * d->J, U2b)) - / d->u[i])) / cabs(clog(d->u[i])); - } - - free(U2b); - free(U2); - return freq_err; +double *harminv_compute_frequency_errors(harminv_data d) { + int i, J2, one = 1; + cmplx *U2, *U2b; + double *freq_err; + + CHECK(d->nfreqs >= 0, "haven't computed eigensolutions yet"); + if (!d->nfreqs) return NULL; + + CHK_MALLOC(freq_err, double, d->nfreqs); + + J2 = d->J * d->J; + CHK_MALLOC(U2, cmplx, J2); + generate_U(U2, NULL, 2, d->c, d->n, d->K, d->J, d->J, d->z, d->z, NULL, NULL, NULL); + CHK_MALLOC(U2b, cmplx, d->J); + + /* For each eigenstate, compute an estimate of the error, roughly + as suggested in W&N, eq. (2.19). */ + + for (i = 0; i < d->nfreqs; ++i) { + cmplx zone = 1.0, zzero = 0.0; + + /* compute U2b = U2 * B[i] */ + ZGEMV(F_("T"), &d->J, &d->J, &zone, U2, &d->J, d->B + i * d->J, &one, &zzero, U2b, &one); + + /* ideally, B[i] should satisfy U2 B[i] = u^2 U0 B[i]. + since B U0 B = 1, then we can get a second estimate + for u by sqrt(B[i] U2 B[i]), and from this we compute + the relative error in the (complex) frequency. */ + + freq_err[i] = cabs(clog(csqrt(symmetric_dot(d->J, d->B + i * d->J, U2b)) / d->u[i])) / + cabs(clog(d->u[i])); + } + + free(U2b); + free(U2); + return freq_err; } /**************************************************************************/ @@ -753,119 +714,96 @@ double *harminv_compute_frequency_errors(harminv_data d) #define UNITY_THRESH 1e-4 /* FIXME? */ /* true if UNITY_THRESH < |u|^n < 1/UNITY_THRESH. */ -static int u_near_unity(cmplx u, int n) -{ - double nlgabsu = n * log(cabs(u)); - return (log(UNITY_THRESH) < nlgabsu && nlgabsu < -log(UNITY_THRESH)); +static int u_near_unity(cmplx u, int n) { + double nlgabsu = n * log(cabs(u)); + return (log(UNITY_THRESH) < nlgabsu && nlgabsu < -log(UNITY_THRESH)); } /* Return an array (of size harminv_get_num_freqs(d)) of complex amplitudes of each sinusoid in the solution. */ -cmplx *harminv_compute_amplitudes(harminv_data d) -{ - int k, j; - cmplx *u; - cmplx *Uu; - cmplx *a; /* the amplitudes of the eigenfrequencies */ - int ku, nu; - - CHECK(d->nfreqs >= 0, "haven't computed eigensolutions yet"); - if (!d->nfreqs) return NULL; - - CHK_MALLOC(a, cmplx, d->nfreqs); - CHK_MALLOC(u, cmplx, d->nfreqs); - - for (k = ku = 0; k < d->nfreqs; ++k) - if (u_near_unity(d->u[k], d->n)) - u[ku++] = d->u[k]; - nu = ku; - - CHK_MALLOC(Uu, cmplx, d->J * nu); - generate_U(Uu, NULL, 0, d->c, d->n, d->K, d->J, nu, d->z, u, - &d->G0, &d->G0_M, &d->D0); - - /* compute the amplitudes via eq. 27 of M&T, except when |u| is - too small..in that case, the computation of Uu is unstable, - and we use eq. 26 instead (which doesn't use half of the data, - but doesn't blow up either): */ - for (k = ku = 0; k < d->nfreqs; ++k) { - cmplxl asum = 0; - if (u_near_unity(d->u[k], d->n)) { /* eq. 27 */ - for (j = 0; j < d->J; ++j) - asum += d->B[k * d->J + j] * Uu[j * nu + ku]; - asum /= d->K; - ku++; - } - else { /* eq. 26 */ - for (j = 0; j < d->J; ++j) - asum += d->B[k * d->J + j] * d->G0[j]; - } - a[k] = asum * asum; - } - - free(Uu); - free(u); - return a; +cmplx *harminv_compute_amplitudes(harminv_data d) { + int k, j; + cmplx *u; + cmplx *Uu; + cmplx *a; /* the amplitudes of the eigenfrequencies */ + int ku, nu; + + CHECK(d->nfreqs >= 0, "haven't computed eigensolutions yet"); + if (!d->nfreqs) return NULL; + + CHK_MALLOC(a, cmplx, d->nfreqs); + CHK_MALLOC(u, cmplx, d->nfreqs); + + for (k = ku = 0; k < d->nfreqs; ++k) + if (u_near_unity(d->u[k], d->n)) u[ku++] = d->u[k]; + nu = ku; + + CHK_MALLOC(Uu, cmplx, d->J * nu); + generate_U(Uu, NULL, 0, d->c, d->n, d->K, d->J, nu, d->z, u, &d->G0, &d->G0_M, &d->D0); + + /* compute the amplitudes via eq. 27 of M&T, except when |u| is + too small..in that case, the computation of Uu is unstable, + and we use eq. 26 instead (which doesn't use half of the data, + but doesn't blow up either): */ + for (k = ku = 0; k < d->nfreqs; ++k) { + cmplxl asum = 0; + if (u_near_unity(d->u[k], d->n)) { /* eq. 27 */ + for (j = 0; j < d->J; ++j) + asum += d->B[k * d->J + j] * Uu[j * nu + ku]; + asum /= d->K; + ku++; + } + else { /* eq. 26 */ + for (j = 0; j < d->J; ++j) + asum += d->B[k * d->J + j] * d->G0[j]; + } + a[k] = asum * asum; + } + + free(Uu); + free(u); + return a; } /**************************************************************************/ -int harminv_get_num_freqs(const harminv_data d) -{ - return d->nfreqs; -} +int harminv_get_num_freqs(const harminv_data d) { return d->nfreqs; } -double harminv_get_freq(harminv_data d, int k) -{ - CHECK(d->nfreqs >= 0, "haven't computed eigensolutions yet"); - CHECK(k >= 0 && k < d->nfreqs, - "argument out of range in harminv_get_freq"); - return(-carg(d->u[k]) / TWOPI); +double harminv_get_freq(harminv_data d, int k) { + CHECK(d->nfreqs >= 0, "haven't computed eigensolutions yet"); + CHECK(k >= 0 && k < d->nfreqs, "argument out of range in harminv_get_freq"); + return (-carg(d->u[k]) / TWOPI); } -double harminv_get_decay(harminv_data d, int k) -{ - CHECK(d->nfreqs >= 0, "haven't computed eigensolutions yet"); - CHECK(k >= 0 && k < d->nfreqs, - "argument out of range in harminv_get_decay"); - return(-log(cabs(d->u[k]))); +double harminv_get_decay(harminv_data d, int k) { + CHECK(d->nfreqs >= 0, "haven't computed eigensolutions yet"); + CHECK(k >= 0 && k < d->nfreqs, "argument out of range in harminv_get_decay"); + return (-log(cabs(d->u[k]))); } -double harminv_get_Q(harminv_data d, int k) -{ - CHECK(k >= 0 && k < d->nfreqs, - "argument out of range in harminv_get_Q"); - return(TWOPI * fabs(harminv_get_freq(d, k)) - / (2 * harminv_get_decay(d, k))); +double harminv_get_Q(harminv_data d, int k) { + CHECK(k >= 0 && k < d->nfreqs, "argument out of range in harminv_get_Q"); + return (TWOPI * fabs(harminv_get_freq(d, k)) / (2 * harminv_get_decay(d, k))); } -void harminv_get_omega(cmplx *omega, harminv_data d, int k) -{ - CHECK(d->nfreqs >= 0, "haven't computed eigensolutions yet"); - CHECK(k >= 0 && k < d->nfreqs, - "argument out of range in harminv_get_omega"); - *omega = (I * clog(d->u[k])); - return; +void harminv_get_omega(cmplx *omega, harminv_data d, int k) { + CHECK(d->nfreqs >= 0, "haven't computed eigensolutions yet"); + CHECK(k >= 0 && k < d->nfreqs, "argument out of range in harminv_get_omega"); + *omega = (I * clog(d->u[k])); + return; } -void harminv_get_amplitude(cmplx *amplitude, harminv_data d, int k) -{ - CHECK(k >= 0 && k < d->nfreqs, - "argument out of range in harminv_get_amplitude"); - if (!d->amps) - d->amps = harminv_compute_amplitudes(d); - *amplitude = d->amps[k]; - return; +void harminv_get_amplitude(cmplx *amplitude, harminv_data d, int k) { + CHECK(k >= 0 && k < d->nfreqs, "argument out of range in harminv_get_amplitude"); + if (!d->amps) d->amps = harminv_compute_amplitudes(d); + *amplitude = d->amps[k]; + return; } -double harminv_get_freq_error(harminv_data d, int k) -{ - CHECK(k >= 0 && k < d->nfreqs, - "argument out of range in harminv_get_freq_error"); - if (!d->errs) - d->errs = harminv_compute_frequency_errors(d); - return d->errs[k]; +double harminv_get_freq_error(harminv_data d, int k) { + CHECK(k >= 0 && k < d->nfreqs, "argument out of range in harminv_get_freq_error"); + if (!d->errs) d->errs = harminv_compute_frequency_errors(d); + return d->errs[k]; } /**************************************************************************/ - diff --git a/harminv.h b/harminv.h index cdf0825..fec9959 100644 --- a/harminv.h +++ b/harminv.h @@ -26,7 +26,7 @@ /**************************************************************************/ #if defined(__cplusplus) -# include +#include extern "C" { typedef std::complex harminv_complex; @@ -44,9 +44,8 @@ typedef struct harminv_data_struct *harminv_data; typedef int (*harminv_mode_ok_func)(harminv_data d, int k, void *); -extern harminv_data harminv_data_create(int n, - const harminv_complex *signal, - double fmin, double fmax, int nf); +extern harminv_data harminv_data_create(int n, const harminv_complex *signal, double fmin, + double fmax, int nf); extern void harminv_data_destroy(harminv_data d); extern void harminv_solve(harminv_data d); @@ -61,10 +60,8 @@ extern double harminv_get_freq_error(harminv_data d, int k); /* the following routines are undocumented and not recommended */ extern void harminv_solve_once(harminv_data d); -extern void harminv_solve_again(harminv_data d, - harminv_mode_ok_func ok, void *ok_d); -extern void harminv_solve_ok_modes(harminv_data d, - harminv_mode_ok_func ok, void *ok_d); +extern void harminv_solve_again(harminv_data d, harminv_mode_ok_func ok, void *ok_d); +extern void harminv_solve_ok_modes(harminv_data d, harminv_mode_ok_func ok, void *ok_d); extern double *harminv_compute_freq_errors(harminv_data d); extern harminv_complex *harminv_compute_amplitudes(harminv_data d); diff --git a/sines.c b/sines.c index 49837ff..b583376 100644 --- a/sines.c +++ b/sines.c @@ -26,10 +26,10 @@ #include "copyright.h" #ifdef HAVE_UNISTD_H -# include +#include #endif #ifdef HAVE_GETOPT_H -# include +#include #endif int verbose = 0; @@ -37,172 +37,159 @@ int verbose = 0; #define TWOPI 6.2831853071795864769252867665590057683943388 #define NPERIODS 10 -void usage(FILE *f) -{ - fprintf(f, "Usage: sines [options] ...\n" - "Note that Im[freq] is the decay rate, or its inverse for -T.\n" - "Options: \n" - " -h : this help message\n" - " -V : print version number and copyright\n" - " -v : verbose output\n" - " -T : specify periods instead of frequencies\n" - " -R : real output\n" - " -r : random amplitudes\n" - " -N : add white noise with amplitude \n" - " -s : use seed for random #s (default =time)\n" - " -n : output points (default %d * max period)\n" - " -t
: time step
(default 1.0)\n", - NPERIODS); +void usage(FILE *f) { + fprintf(f, + "Usage: sines [options] ...\n" + "Note that Im[freq] is the decay rate, or its inverse for -T.\n" + "Options: \n" + " -h : this help message\n" + " -V : print version number and copyright\n" + " -v : verbose output\n" + " -T : specify periods instead of frequencies\n" + " -R : real output\n" + " -r : random amplitudes\n" + " -N : add white noise with amplitude \n" + " -s : use seed for random #s (default =time)\n" + " -n : output points (default %d * max period)\n" + " -t
: time step
(default 1.0)\n", + NPERIODS); } typedef struct { - double freq, decay, amplitude, phase; + double freq, decay, amplitude, phase; } sinusoid; -#define MAX2(a,b) ((a) > (b) ? (a) : (b)) - -int main(int argc, char **argv) -{ - int c; - extern char *optarg; - extern int optind; - int iarg; - int specify_periods = 0; - int random_amplitudes = 0; - int real_output = 0; - sinusoid *sines = NULL; - int nsines = 0, nalloc = 0, n = 0; - double max_period = 0; - double dt = 1.0; - double noise = 0.0; - int i, is; - - srand(time(NULL)); - - while ((c = getopt(argc, argv, "hVvTRrn:t:N:s:")) != -1) - switch (c) { - case 'h': - usage(stdout); - return EXIT_SUCCESS; - case 'V': - printf("sines " PACKAGE_VERSION " by Steven G. Johnson.\n" - "Test program for harminv.\n" - COPYRIGHT); - return EXIT_SUCCESS; - case 'v': - verbose = 1; - break; - case 'T': - specify_periods = 1; - break; - case 'R': - real_output = 1; - break; - case 'r': - random_amplitudes = 1; - break; - case 'N': - noise = atof(optarg); - break; - case 's': - srand(atoi(optarg)); - break; - case 'n': - n = atoi(optarg); - if (n < 1) { - fprintf(stderr, "sines: " - "invalid non-positive -n argument %d\n", n); - return EXIT_FAILURE; - } - break; - case 't': - dt = atof(optarg); - break; - default: - fprintf(stderr, "sines: invalid argument -%c\n", c); - usage(stderr); - return EXIT_FAILURE; - } - - if (optind == argc) { /* no parameters left */ - fprintf(stderr, "sines: missing required frequency(ies)\n"); - usage(stderr); +#define MAX2(a, b) ((a) > (b) ? (a) : (b)) + +int main(int argc, char **argv) { + int c; + extern char *optarg; + extern int optind; + int iarg; + int specify_periods = 0; + int random_amplitudes = 0; + int real_output = 0; + sinusoid *sines = NULL; + int nsines = 0, nalloc = 0, n = 0; + double max_period = 0; + double dt = 1.0; + double noise = 0.0; + int i, is; + + srand(time(NULL)); + + while ((c = getopt(argc, argv, "hVvTRrn:t:N:s:")) != -1) + switch (c) { + case 'h': usage(stdout); return EXIT_SUCCESS; + case 'V': + printf("sines " PACKAGE_VERSION " by Steven G. Johnson.\n" + "Test program for harminv.\n" COPYRIGHT); + return EXIT_SUCCESS; + case 'v': verbose = 1; break; + case 'T': specify_periods = 1; break; + case 'R': real_output = 1; break; + case 'r': random_amplitudes = 1; break; + case 'N': noise = atof(optarg); break; + case 's': srand(atoi(optarg)); break; + case 'n': + n = atoi(optarg); + if (n < 1) { + fprintf(stderr, + "sines: " + "invalid non-positive -n argument %d\n", + n); return EXIT_FAILURE; - } - - for (iarg = optind; iarg < argc; ++iarg) { - sinusoid s = { 0, 0, 0, 0 }; - - if (sscanf(argv[iarg], "%lf+%lfi", &s.freq, &s.decay) < 1) { - fprintf(stderr, "sines: invalid argument \"%s\"\n", - argv[iarg]); - return EXIT_FAILURE; - } - if (specify_periods) { - if (s.freq == 0) { - fprintf(stderr, "sines: invalid argument \"%s\"" - ": 0 not a valid period\n", argv[iarg]); - return EXIT_FAILURE; - } - s.freq = 1/s.freq; - if (s.decay != 0) - s.decay = 1/s.decay; - } - - if (s.decay == 0 || fabs(1/s.freq) > 1/s.decay) - max_period = MAX2(max_period, fabs(1/s.freq)); - else - max_period = MAX2(max_period, 1/s.decay); - - if (random_amplitudes) { - s.amplitude = rand() * 1.0/RAND_MAX; - s.phase = rand() * TWOPI/RAND_MAX - TWOPI/2; - } - else { - s.amplitude = (iarg - optind + 1); - s.phase = (iarg - optind + 1) * TWOPI / (argc-optind) - TWOPI/2; - } - - if (verbose) - printf("# mode: frequency = %g (period %g), decay = %g (lifetime %g), amplitude = %g, phase = %g\n", s.freq, 1/s.freq, s.decay, s.decay != 0 ? 1/s.decay : 0, s.amplitude, s.phase); - - if (nsines >= nalloc) { - nalloc = (nalloc + 1) * 2; - sines = (sinusoid*) realloc(sines, sizeof(sinusoid) * nalloc); - CHECK(sines, "out of memory"); - } - sines[nsines++] = s; - } - - if (n == 0 && max_period == 0) { - fprintf(stderr, "sines: must specify -n or non-zero period\n"); - return EXIT_FAILURE; - } - - if (n == 0) - n = (int) (max_period / fabs(dt) + 0.5) * NPERIODS; - - for (i = 0; i < n; ++i) { - cmplx output = 0; - - for (is = 0; is < nsines; ++is) - output += sines[is].amplitude * - cexp(I * (sines[is].phase - TWOPI*sines[is].freq * i*dt) - - sines[is].decay * i*dt); - output += noise * (rand() * 2.0/RAND_MAX - 1); - if (real_output) - printf("%0.17e\n", creal(output)); - else - printf("%0.17e%+0.17ei\n", creal(output), cimag(output)); - } - - free(sines); - return EXIT_SUCCESS; + } + break; + case 't': dt = atof(optarg); break; + default: + fprintf(stderr, "sines: invalid argument -%c\n", c); + usage(stderr); + return EXIT_FAILURE; + } + + if (optind == argc) { /* no parameters left */ + fprintf(stderr, "sines: missing required frequency(ies)\n"); + usage(stderr); + return EXIT_FAILURE; + } + + for (iarg = optind; iarg < argc; ++iarg) { + sinusoid s = {0, 0, 0, 0}; + + if (sscanf(argv[iarg], "%lf+%lfi", &s.freq, &s.decay) < 1) { + fprintf(stderr, "sines: invalid argument \"%s\"\n", argv[iarg]); + return EXIT_FAILURE; + } + if (specify_periods) { + if (s.freq == 0) { + fprintf(stderr, + "sines: invalid argument \"%s\"" + ": 0 not a valid period\n", + argv[iarg]); + return EXIT_FAILURE; + } + s.freq = 1 / s.freq; + if (s.decay != 0) s.decay = 1 / s.decay; + } + + if (s.decay == 0 || fabs(1 / s.freq) > 1 / s.decay) + max_period = MAX2(max_period, fabs(1 / s.freq)); + else + max_period = MAX2(max_period, 1 / s.decay); + + if (random_amplitudes) { + s.amplitude = rand() * 1.0 / RAND_MAX; + s.phase = rand() * TWOPI / RAND_MAX - TWOPI / 2; + } + else { + s.amplitude = (iarg - optind + 1); + s.phase = (iarg - optind + 1) * TWOPI / (argc - optind) - TWOPI / 2; + } + + if (verbose) + printf("# mode: frequency = %g (period %g), decay = %g (lifetime %g), amplitude = %g, phase " + "= %g\n", + s.freq, 1 / s.freq, s.decay, s.decay != 0 ? 1 / s.decay : 0, s.amplitude, s.phase); + + if (nsines >= nalloc) { + nalloc = (nalloc + 1) * 2; + sines = (sinusoid *)realloc(sines, sizeof(sinusoid) * nalloc); + CHECK(sines, "out of memory"); + } + sines[nsines++] = s; + } + + if (n == 0 && max_period == 0) { + fprintf(stderr, "sines: must specify -n or non-zero period\n"); + return EXIT_FAILURE; + } + + if (n == 0) n = (int)(max_period / fabs(dt) + 0.5) * NPERIODS; + + for (i = 0; i < n; ++i) { + cmplx output = 0; + + for (is = 0; is < nsines; ++is) + output += sines[is].amplitude * cexp(I * (sines[is].phase - TWOPI * sines[is].freq * i * dt) - + sines[is].decay * i * dt); + output += noise * (rand() * 2.0 / RAND_MAX - 1); + if (real_output) + printf("%0.17e\n", creal(output)); + else + printf("%0.17e%+0.17ei\n", creal(output), cimag(output)); + } + + free(sines); + return EXIT_SUCCESS; } #ifdef F77_DUMMY_MAIN -# ifdef __cplusplus +#ifdef __cplusplus extern "C" -# endif -int F77_DUMMY_MAIN() { return 1; } +#endif + int + F77_DUMMY_MAIN() { + return 1; +} #endif