diff options
Diffstat (limited to 'math/curfit/cvsolve.gx')
| -rw-r--r-- | math/curfit/cvsolve.gx | 51 |
1 files changed, 51 insertions, 0 deletions
diff --git a/math/curfit/cvsolve.gx b/math/curfit/cvsolve.gx new file mode 100644 index 00000000..08e0bf92 --- /dev/null +++ b/math/curfit/cvsolve.gx @@ -0,0 +1,51 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math/curfit.h> + +$if (datatype == r) +include "curfitdef.h" +$else +include "dcurfitdef.h" +$endif + + +# CVSOLVE -- Solve the matrix normal equations of the form ca = b for a, +# where c is a symmetric, positive semi-definite, banded matrix with +# CV_NCOEFF(cv) rows and a and b are CV_NCOEFF(cv)-vectors. +# Initially c is stored in the CV_ORDER(cv) by CV_NCOEFF(cv) matrix MATRIX +# and b is stored in VECTOR. +# The Cholesky factorization of MATRIX is calculated and stored in CHOFAC. +# Finally the coefficients are calculated by forward and back substitution +# and stored in COEFF. + +$if (datatype == r) +procedure cvsolve (cv, ier) +$else +procedure dcvsolve (cv, ier) +$endif + + +pointer cv # curve descriptor +int ier # ier = OK, everything OK + # ier = SINGULAR, matrix is singular, 1 or more + # coefficients are 0. + # ier = NO_DEG_FREEDOM, too few points to solve matrix +int nfree + +begin + ier = OK + nfree = CV_NPTS(cv) - CV_NCOEFF(cv) + + if (nfree < 0) { + ier = NO_DEG_FREEDOM + return + } + + # calculate the Cholesky factorization of the data matrix + call $tcvchofac (MATRIX(CV_MATRIX(cv)), CV_ORDER(cv), CV_NCOEFF(cv), + CHOFAC(CV_CHOFAC(cv)), ier) + + # solve for the coefficients by forward and back substitution + call $tcvchoslv (CHOFAC(CV_CHOFAC(cv)), CV_ORDER(cv), CV_NCOEFF(cv), + VECTOR(CV_VECTOR(cv)), COEFF(CV_COEFF(cv))) +end |