diff options
author | William Joye <wjoye@cfa.harvard.edu> | 2016-11-02 19:11:06 (GMT) |
---|---|---|
committer | William Joye <wjoye@cfa.harvard.edu> | 2016-11-02 19:11:06 (GMT) |
commit | 33c55bd916dff8c4932b01c7db58f0103ac31c31 (patch) | |
tree | a4cdca3287dd2df5247ce8079c424ffa438b4c2e /ast/cminpack/qrsolv.c | |
parent | 4121637f3d41d6dc23e6543a445b5a3aed9e6ddc (diff) | |
download | blt-33c55bd916dff8c4932b01c7db58f0103ac31c31.zip blt-33c55bd916dff8c4932b01c7db58f0103ac31c31.tar.gz blt-33c55bd916dff8c4932b01c7db58f0103ac31c31.tar.bz2 |
update ast
Diffstat (limited to 'ast/cminpack/qrsolv.c')
-rw-r--r-- | ast/cminpack/qrsolv.c | 218 |
1 files changed, 0 insertions, 218 deletions
diff --git a/ast/cminpack/qrsolv.c b/ast/cminpack/qrsolv.c deleted file mode 100644 index 6ab9e98..0000000 --- a/ast/cminpack/qrsolv.c +++ /dev/null @@ -1,218 +0,0 @@ -#include "cminpack.h" -#include <math.h> -#include "cminpackP.h" - -__cminpack_attr__ -void __cminpack_func__(qrsolv)(int n, real *r, int ldr, - const int *ipvt, const real *diag, const real *qtb, real *x, - real *sdiag, real *wa) -{ - /* Initialized data */ - -#define p5 .5 -#define p25 .25 - - /* Local variables */ - int i, j, k, l; - real cos, sin, sum, temp; - int nsing; - real qtbpj; - -/* ********** */ - -/* subroutine qrsolv */ - -/* given an m by n matrix a, an n by n diagonal matrix d, */ -/* and an m-vector b, the problem is to determine an x which */ -/* solves the system */ - -/* a*x = b , d*x = 0 , */ - -/* in the least squares sense. */ - -/* this subroutine completes the solution of the problem */ -/* if it is provided with the necessary information from the */ -/* qr factorization, with column pivoting, of a. that is, if */ -/* a*p = q*r, where p is a permutation matrix, q has orthogonal */ -/* columns, and r is an upper triangular matrix with diagonal */ -/* elements of nonincreasing magnitude, then qrsolv expects */ -/* the full upper triangle of r, the permutation matrix p, */ -/* and the first n components of (q transpose)*b. the system */ -/* a*x = b, d*x = 0, is then equivalent to */ - -/* t t */ -/* r*z = q *b , p *d*p*z = 0 , */ - -/* where x = p*z. if this system does not have full rank, */ -/* then a least squares solution is obtained. on output qrsolv */ -/* also provides an upper triangular matrix s such that */ - -/* t t t */ -/* p *(a *a + d*d)*p = s *s . */ - -/* s is computed within qrsolv and may be of separate interest. */ - -/* the subroutine statement is */ - -/* subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa) */ - -/* where */ - -/* n is a positive integer input variable set to the order of r. */ - -/* r is an n by n array. on input the full upper triangle */ -/* must contain the full upper triangle of the matrix r. */ -/* on output the full upper triangle is unaltered, and the */ -/* strict lower triangle contains the strict upper triangle */ -/* (transposed) of the upper triangular matrix s. */ - -/* ldr is a positive integer input variable not less than n */ -/* which specifies the leading dimension of the array r. */ - -/* ipvt is an integer input array of length n which defines the */ -/* permutation matrix p such that a*p = q*r. column j of p */ -/* is column ipvt(j) of the identity matrix. */ - -/* diag is an input array of length n which must contain the */ -/* diagonal elements of the matrix d. */ - -/* qtb is an input array of length n which must contain the first */ -/* n elements of the vector (q transpose)*b. */ - -/* x is an output array of length n which contains the least */ -/* squares solution of the system a*x = b, d*x = 0. */ - -/* sdiag is an output array of length n which contains the */ -/* diagonal elements of the upper triangular matrix s. */ - -/* wa is a work array of length n. */ - -/* subprograms called */ - -/* fortran-supplied ... dabs,dsqrt */ - -/* argonne national laboratory. minpack project. march 1980. */ -/* burton s. garbow, kenneth e. hillstrom, jorge j. more */ - -/* ********** */ - -/* copy r and (q transpose)*b to preserve input and initialize s. */ -/* in particular, save the diagonal elements of r in x. */ - - for (j = 0; j < n; ++j) { - for (i = j; i < n; ++i) { - r[i + j * ldr] = r[j + i * ldr]; - } - x[j] = r[j + j * ldr]; - wa[j] = qtb[j]; - } - -/* eliminate the diagonal matrix d using a givens rotation. */ - - for (j = 0; j < n; ++j) { - -/* prepare the row of d to be eliminated, locating the */ -/* diagonal element using p from the qr factorization. */ - - l = ipvt[j]-1; - if (diag[l] != 0.) { - for (k = j; k < n; ++k) { - sdiag[k] = 0.; - } - sdiag[j] = diag[l]; - -/* the transformations to eliminate the row of d */ -/* modify only a single element of (q transpose)*b */ -/* beyond the first n, which is initially zero. */ - - qtbpj = 0.; - for (k = j; k < n; ++k) { - -/* determine a givens rotation which eliminates the */ -/* appropriate element in the current row of d. */ - - if (sdiag[k] != 0.) { -# ifdef USE_LAPACK - dlartg_( &r[k + k * ldr], &sdiag[k], &cos, &sin, &temp ); -# else /* !USE_LAPACK */ - if (fabs(r[k + k * ldr]) < fabs(sdiag[k])) { - real cotan; - cotan = r[k + k * ldr] / sdiag[k]; - sin = p5 / sqrt(p25 + p25 * (cotan * cotan)); - cos = sin * cotan; - } else { - real tan; - tan = sdiag[k] / r[k + k * ldr]; - cos = p5 / sqrt(p25 + p25 * (tan * tan)); - sin = cos * tan; - } - -/* compute the modified diagonal element of r and */ -/* the modified element of ((q transpose)*b,0). */ - -# endif /* !USE_LAPACK */ - temp = cos * wa[k] + sin * qtbpj; - qtbpj = -sin * wa[k] + cos * qtbpj; - wa[k] = temp; - -/* accumulate the tranformation in the row of s. */ -# ifdef USE_CBLAS - cblas_drot( n-k, &r[k + k * ldr], 1, &sdiag[k], 1, cos, sin ); -# else /* !USE_CBLAS */ - r[k + k * ldr] = cos * r[k + k * ldr] + sin * sdiag[k]; - if (n > k+1) { - for (i = k+1; i < n; ++i) { - temp = cos * r[i + k * ldr] + sin * sdiag[i]; - sdiag[i] = -sin * r[i + k * ldr] + cos * sdiag[i]; - r[i + k * ldr] = temp; - } - } -# endif /* !USE_CBLAS */ - } - } - } - -/* store the diagonal element of s and restore */ -/* the corresponding diagonal element of r. */ - - sdiag[j] = r[j + j * ldr]; - r[j + j * ldr] = x[j]; - } - -/* solve the triangular system for z. if the system is */ -/* singular, then obtain a least squares solution. */ - - nsing = n; - for (j = 0; j < n; ++j) { - if (sdiag[j] == 0. && nsing == n) { - nsing = j; - } - if (nsing < n) { - wa[j] = 0.; - } - } - if (nsing >= 1) { - for (k = 1; k <= nsing; ++k) { - j = nsing - k; - sum = 0.; - if (nsing > j+1) { - for (i = j+1; i < nsing; ++i) { - sum += r[i + j * ldr] * wa[i]; - } - } - wa[j] = (wa[j] - sum) / sdiag[j]; - } - } - -/* permute the components of z back to components of x. */ - - for (j = 0; j < n; ++j) { - l = ipvt[j]-1; - x[l] = wa[j]; - } - return; - -/* last card of subroutine qrsolv. */ - -} /* qrsolv_ */ - |