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@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
set_orthogonal(int n,double * Q,double * v)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
matrix_multiplication(int n,double * A,double * Q,double * R)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
matrix_transposition(int n,double * A)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
print_matrix(int m,int n,double * a)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
main(int argc,char * argv[])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