## File: matrix.c

package info (click to toggle)
gnuplot 3.7.1p1-4
• area: main
• in suites: potato
• size: 5,808 kB
• ctags: 4,418
• sloc: ansic: 43,279; lisp: 661; makefile: 647; asm: 539; sh: 416; objc: 379; csh: 297; pascal: 194; perl: 138
 file content (269 lines) | stat: -rw-r--r-- 7,152 bytes parent folder | download
 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269 #ifndef lint static char *RCSid = "$Id: matrix.c,v 1.12 1998/11/19 10:40:31 lhecking Exp$"; #endif /* * Matrix algebra, part of * * Nonlinear least squares fit according to the * Marquardt-Levenberg-algorithm * * added as Patch to Gnuplot (v3.2 and higher) * by Carsten Grammes * Experimental Physics, University of Saarbruecken, Germany * * Internet address: cagr@rz.uni-sb.de * * Copyright of this module: Carsten Grammes, 1993 * * Permission to use, copy, and distribute this software and its * documentation for any purpose with or without fee is hereby granted, * provided that the above copyright notice appear in all copies and * that both that copyright notice and this permission notice appear * in supporting documentation. * * This software is provided "as is" without express or implied warranty. */ #include "fit.h" #include "matrix.h" #include "stdfn.h" #include "alloc.h" /*****************************************************************/ #define Swap(a,b) {double temp = (a); (a) = (b); (b) = temp;} #define WINZIG 1e-30 /***************************************************************** internal prototypes *****************************************************************/ static GP_INLINE int fsign __PROTO((double x)); /***************************************************************** first straightforward vector and matrix allocation functions *****************************************************************/ double *vec (n) int n; { /* allocates a double vector with n elements */ double *dp; if( n < 1 ) return (double *) NULL; dp = (double *) gp_alloc ( n * sizeof(double), "vec"); return dp; } double **matr (rows, cols) int rows; int cols; { /* allocates a double matrix */ register int i; register double **m; if ( rows < 1 || cols < 1 ) return NULL; m = (double **) gp_alloc ( rows * sizeof(double *) , "matrix row pointers"); m[0] = (double *) gp_alloc ( rows * cols * sizeof(double), "matrix elements"); for ( i = 1; i0 ? 1 : (x < 0) ? -1 : 0) ; } /***************************************************************** Solve least squares Problem C*x+d = r, |r| = min!, by Given rotations (QR-decomposition). Direct implementation of the algorithm presented in H.R.Schwarz: Numerische Mathematik, 'equation' number (7.33) If 'd == NULL', d is not accesed: the routine just computes the QR decomposition of C and exits. If 'want_r == 0', r is not rotated back (\hat{r} is returned instead). *****************************************************************/ void Givens (C, d, x, r, N, n, want_r) double **C; double *d; double *x; double *r; int N; int n; int want_r; { int i, j, k; double w, gamma, sigma, rho, temp; double epsilon = DBL_EPSILON; /* FIXME (?)*/ /* * First, construct QR decomposition of C, by 'rotating away' * all elements of C below the diagonal. The rotations are * stored in place as Givens coefficients rho. * Vector d is also rotated in this same turn, if it exists */ for (j = 0; j= 0; i--) { /* solve R*x+d = 0, by backsubstitution */ double s = d[i]; r[i] = 0; /* ... and also set r[i] = 0 for i= 0; j--) for (i = N-1; i >= 0; i--) { if ((rho = C[i][j]) == 1) { /* reconstruct gamma, sigma from stored rho */ gamma = 0; sigma = 1; } else if (fabs(rho)<1) { sigma = rho; gamma = sqrt(1-sigma*sigma); } else { gamma = 1/fabs(rho); sigma = fsign(rho)*sqrt(1-gamma*gamma); } temp = gamma*r[j] + sigma*r[i]; /* rotate back indices (i,j) */ r[i] = -sigma*r[j] + gamma*r[i]; r[j] = temp; } } /* Given a triangular Matrix R, compute (R^T * R)^(-1), by forward * then back substitution * * R, I are n x n Matrices, I is for the result. Both must already be * allocated. * * Will only calculate the lower triangle of I, as it is symmetric */ void Invert_RtR ( R, I, n) double **R; double **I; int n; { int i, j, k; /* fill in the I matrix, and check R for regularity : */ for (i = 0; i= k; i--) { /* don't compute upper triangle of A */ double s = I[i][k]; for (j = i+1; j