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
