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