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