Re: gsl_multifit

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

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]);
  }
}