This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: gsl_multifit
Hi,
> You could try computing the solution using the following
> functions instead and see if you get any speed improvement.
> gsl_linalg_QR_decomp
> gsl_linalg_QR_lssolve
Thanks for the tip. It's much faster now, but still slower than octave.
Order: (300,200) (600,400) (900,600) (1200,1000) (1500,1200)
Octave: 0.340s 3.320s 11.580s 1m9.680s 2m1.090s
gsl(QR): 0.390s 4.560s 22.690s 1m38.640s 5m1.580s
--
:wq `Kai Trukenmueller'
/*
* sinft - Sinus Fourier Transformation
*
* f(x) = sum_n^NoF sin(omega_n x) * F
* gesucht ist F.
* In Matixform:
* f = S*F, f:(NoX,1), S:(NoX, NoF), F:(1, NoF)
* mit der Matrix (S)_ij = sin(x_j * omega_j)
*
* --> F = (S^T * S)^-1 * S^T *f
* solved using `QR - decomposition' of
* gsl (GNU Scientific Library)
*/
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_multifit.h>
/* `F' Frequency array with NoF elements
* `f' spacearray mit `NoX' elements
* x : monotonic [0 bis `L'], with `NoX' elements
* `omega': frequencyarray mit NoF elements */
void sinft (double *F, double *f, double *omega, int NoX, double L, int NoF)
{
int i, j;
double xi;
gsl_matrix *A;
gsl_vector *b, *x, *tau, *residual;
A = gsl_matrix_alloc (NoX, NoF);
b = gsl_vector_alloc (NoX);
x = gsl_vector_alloc (NoF);
tau = gsl_vector_alloc (NoF);
residual = gsl_vector_alloc (NoX);
for (i=0; i<NoX; i++) {
xi = (float)i * L / ((double)NoX - 1.0);
for (j=0; j<NoF; j++) {
gsl_matrix_set (A, i, j, sin(xi * omega[j]));
}
gsl_vector_set (b, i, f[i]);
}
gsl_linalg_QR_decomp (A, tau);
gsl_linalg_QR_lssolve (A, tau, b, x, residual);
for (i=0; i<NoF; i++) {
F[i] = gsl_vector_get (x, i);
}
}
/*
* arguments:
* argv[1] : Number of Space - Elements
* argv[2] : Number of Frequency - Elements
*/
int main (int args, char **argv)
{
int i, j;
int N, M;
double L;
double xi;
double f[2000];
double F[1000], omega[1000];
double *fp, *Fp, *Omega;
L = 1.0;
if (args <2) {
printf("argv[1] : Number of space elements\n");
printf("argv[2] : Number of frequency elements (< argv[1])\n");
exit(1);
}
N = atoi (argv[1]);
M = atoi (argv[2]);
fp = f;
Fp = F;
Omega = omega;
for (i=0; i<M; i++) {
omega[i] = ((double)i+1.0) * M_PI / L;
}
/* Spacefunctions */
for (i=0; i<N; i++) {
xi = (double)i * L / ((double)N - 1.0);
f[i] = pow(xi, 2);
//printf("x%d = %f, f = %f\n", i, xi, f[i]);
}
sinft (F, f, Omega, N, L, M);
for (i=0; i<M; i++) {
printf("%f\n",F[i]);
}
}