update ast 8.6.2

1 files changed, 218 insertions, 0 deletions

diff --git a/ast/cminpack/qrsolv.c b/ast/cminpack/qrsolv.c
new file mode 100644
index 0000000..6ab9e98
--- /dev/null
+++ b/ast/cminpack/qrsolv.c
@@ -0,0 +1,218 @@
+#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_ */
+