Re: gsl_multifit

[Kai Trukenmueller <trukenm at ag2 dot mechanik dot tu-darmstadt dot de>]

Hi,

I benchmarked the programs on system with Athlon 1GHz, 640MB running Linux-2.4.10 with the following results:

Dimensions	(120,80)  (300,200)  (600,400)  (900,600)  (1200,800)
Octave-Script	0.060s    0.410s      3.390s     12.260s    39.790s
Compiled gsl    0.110s    2.240s      31.310s    2m7.590s   seg.fault

The gsl routine not just turns out to be much slower, it also seems to
be `instable' for higher orders (the results do not convege -> inft..)
Maybe sth. is wrong in my code. For low-orders (~<500) both results are
equivalent.

I'm not using the atlas-blas but gslcblas.

Code attached;
The octave-script `sinft.m' should be executable, and can be started
directely from the shell (iff /usr/bin/octave exists).
The matrix dimensions must be edited in the source code (in the c
programm, they are arguments).


-- 
: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 `multi parameter fitting Method' 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 chisq, xi;
  gsl_matrix *X, *cov;
  gsl_vector *y, *w, *c;
  gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc (NoX, NoF);

  X = gsl_matrix_alloc (NoX, NoF);
  y = gsl_vector_alloc (NoX);
  w = gsl_vector_alloc (NoX);
  c = gsl_vector_alloc (NoF);
  cov = gsl_matrix_alloc (NoF, NoF);
  for (i=0; i<NoX; i++) {
    xi = (float)i * L / ((double)NoX - 1.0);
    for (j=0; j<NoF; j++) {
      gsl_matrix_set (X, i, j, sin(xi * omega[j]));
    }
    gsl_vector_set (y, i, f[i]);
  }

  gsl_multifit_linear (X, y, c, cov, &chisq, work);
  gsl_multifit_linear_free (work);

  for (i=0; i<NoF; i++) {
    F[i] = gsl_vector_get (c, 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[1000];
  double F[100], omega[100];
  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]);
  }
}
#! /usr/bin/octave -q
% Sine - Fourier
% f(x) = sum_n^NoX F_n sin(omega_n x_j)
%
% as an example for f = x^2
%
% results in solvin f = S F
% -->  F = inv (S' * S) * S' * f
% with S(i,j) = sin (omega(j) * x(i))
NoX = 120;
NoM = 80;

j = 1: NoX;
n = 1: NoM;

x = 0 : 1/(NoX-1) : 1;
omega = n*pi/1;

S=sin(x' * omega);

%In this exaple, is f a quadratic polynom in space
f = [x.^2]';

F = S\f;
printf("%f\n",F);

%Reverse transform to compare f
%f_approx = S * F;
%printf("%f\n",f_approx);
%printf("\n");
%printf("%f\n",x);