From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- math/curfit/cvsolver.x | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 math/curfit/cvsolver.x (limited to 'math/curfit/cvsolver.x') diff --git a/math/curfit/cvsolver.x b/math/curfit/cvsolver.x new file mode 100644 index 00000000..b52f012e --- /dev/null +++ b/math/curfit/cvsolver.x @@ -0,0 +1,43 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include + +include "curfitdef.h" + + +# 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. + +procedure cvsolve (cv, ier) + + +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 rcvchofac (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 rcvchoslv (CHOFAC(CV_CHOFAC(cv)), CV_ORDER(cv), CV_NCOEFF(cv), + VECTOR(CV_VECTOR(cv)), COEFF(CV_COEFF(cv))) +end -- cgit