Home | History | Annotate | Download | only in test 
      1 /*
      2  * Library:  lmfit (Levenberg-Marquardt least squares fitting)
      3  *
      4  * File:     demo/curve1.c
      5  *
      6  * Contents: Example for curve fitting with lmcurve():
      7  *           fit a data set y(x) by a curve f(x;p).
      8  *
      9  * Note:     Any modification of this example should be copied to
     10  *           the manual page source lmcurve.pod and to the wiki.
     11  *
     12  * Author:   Joachim Wuttke <j.wuttke (at) fz-juelich.de> 2004-2013
     13  *
     14  * Licence:  see ../COPYING (FreeBSD)
     15  *
     16  * Homepage: apps.jcns.fz-juelich.de/lmfit
     17  */
     18 
     19 #include "lmcurve.h"
     20 #include <stdio.h>
     21 #include <stdlib.h>
     22 
     23 void lm_qrfac(int m, int n, double *a, int *ipvt,
     24               double *rdiag, double *acnorm, double *wa);
     25 
     26 void set_orthogonal( int n, double *Q, double* v )
     27 {
     28     int i, j;
     29     double temp = 0;
     30     for (i=0; i<n; ++i)
     31         temp += v[i]*v[i];
     32     for (i=0; i<n; ++i)
     33         for (j=0; j<n; ++j)
     34             Q[j*n+i] = ( i==j ? 1 : 0 ) - 2*v[i]*v[j]/temp;
     35 }
     36 
     37 void matrix_multiplication( int n, double *A, double *Q, double *R )
     38 {
     39     int i,j,k;
     40     double temp;
     41     for (i=0; i<n; ++i) {
     42         for (j=0; j<n; ++j) {
     43             temp = 0;
     44             for (k=0; k<n; ++k) {
     45                 temp += Q[k*n+i]*R[j*n+k];
     46             }
     47             A[j*n+i] = temp;
     48         }
     49     }
     50 }
     51 
     52 void matrix_transposition( int n, double *A )
     53 {
     54     int i,j;
     55     double temp;
     56     for (i=0; i<n; ++i) {
     57         for (j=i+1; j<n; ++j) {
     58             temp = A[j*n+i];
     59             A[j*n+i] = A[i*n+j];
     60             A[i*n+j] = temp;
     61         }
     62     }
     63 }
     64 
     65 void print_matrix(int m, int n, double *a)
     66 {
     67     int i,j;
     68     for (i=0; i<m; ++i) {
     69         for (j=0; j<n; ++j) {
     70             printf( "%8.4g ", a[j*m+i] );
     71         }
     72         printf ("\n");
     73     }
     74 }
     75 
     76 int main( int argc, char *argv[] )
     77 {
     78     double A[9], rdiag[3], acnorm[3], wa[3];
     79     int i, ipvt[3];
     80 
     81     if ( argc!= 10 ) {
     82         fprintf( stderr, "bad # args\n" );
     83         exit(1);
     84     }
     85     for ( int i=0; i<9; ++i )
     86         A[i] = atof( argv[1+i] );
     87     matrix_transposition( 3, A );
     88 
     89     printf( "Input matrix A:\n" );
     90     print_matrix(3, 3, A);
     91 
     92     lm_qrfac(3, 3, A, ipvt, rdiag, acnorm, wa);
     93 
     94     printf( "Output matrix A:\n" );
     95     print_matrix(3, 3, A);
     96 
     97     printf( "Output vector rdiag:\n" );
     98     print_matrix(1, 3, rdiag);
     99 
    100     printf( "Output vector acnorm:\n" );
    101     print_matrix(1, 3, acnorm);
    102 
    103     printf( "Output vector ipvt:\n" );
    104     for (i=0; i<3; ++i)
    105         printf( "%i ", ipvt[i] );
    106     printf("\n");
    107 }
    108