summaryrefslogtreecommitdiff
path: root/gsl-1.9/linalg/svd.c
diff options
context:
space:
mode:
Diffstat (limited to 'gsl-1.9/linalg/svd.c')
-rw-r--r--gsl-1.9/linalg/svd.c620
1 files changed, 620 insertions, 0 deletions
diff --git a/gsl-1.9/linalg/svd.c b/gsl-1.9/linalg/svd.c
new file mode 100644
index 0000000..cd5588a
--- /dev/null
+++ b/gsl-1.9/linalg/svd.c
@@ -0,0 +1,620 @@
+/* linalg/svd.c
+ *
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman, Brian Gough
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <config.h>
+#include <stdlib.h>
+#include <string.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+
+#include <gsl/gsl_linalg.h>
+
+#include "givens.c"
+#include "svdstep.c"
+
+/* Factorise a general M x N matrix A into,
+ *
+ * A = U D V^T
+ *
+ * where U is a column-orthogonal M x N matrix (U^T U = I),
+ * D is a diagonal N x N matrix,
+ * and V is an N x N orthogonal matrix (V^T V = V V^T = I)
+ *
+ * U is stored in the original matrix A, which has the same size
+ *
+ * V is stored as a separate matrix (not V^T). You must take the
+ * transpose to form the product above.
+ *
+ * The diagonal matrix D is stored in the vector S, D_ii = S_i
+ */
+
+int
+gsl_linalg_SV_decomp (gsl_matrix * A, gsl_matrix * V, gsl_vector * S,
+ gsl_vector * work)
+{
+ size_t a, b, i, j;
+
+ const size_t M = A->size1;
+ const size_t N = A->size2;
+ const size_t K = GSL_MIN (M, N);
+
+ if (M < N)
+ {
+ GSL_ERROR ("svd of MxN matrix, M<N, is not implemented", GSL_EUNIMPL);
+ }
+ else if (V->size1 != N)
+ {
+ GSL_ERROR ("square matrix V must match second dimension of matrix A",
+ GSL_EBADLEN);
+ }
+ else if (V->size1 != V->size2)
+ {
+ GSL_ERROR ("matrix V must be square", GSL_ENOTSQR);
+ }
+ else if (S->size != N)
+ {
+ GSL_ERROR ("length of vector S must match second dimension of matrix A",
+ GSL_EBADLEN);
+ }
+ else if (work->size != N)
+ {
+ GSL_ERROR ("length of workspace must match second dimension of matrix A",
+ GSL_EBADLEN);
+ }
+
+ /* Handle the case of N = 1 (SVD of a column vector) */
+
+ if (N == 1)
+ {
+ gsl_vector_view column = gsl_matrix_column (A, 0);
+ double norm = gsl_blas_dnrm2 (&column.vector);
+
+ gsl_vector_set (S, 0, norm);
+ gsl_matrix_set (V, 0, 0, 1.0);
+
+ if (norm != 0.0)
+ {
+ gsl_blas_dscal (1.0/norm, &column.vector);
+ }
+
+ return GSL_SUCCESS;
+ }
+
+ {
+ gsl_vector_view f = gsl_vector_subvector (work, 0, K - 1);
+
+ /* bidiagonalize matrix A, unpack A into U S V */
+
+ gsl_linalg_bidiag_decomp (A, S, &f.vector);
+ gsl_linalg_bidiag_unpack2 (A, S, &f.vector, V);
+
+ /* apply reduction steps to B=(S,Sd) */
+
+ chop_small_elements (S, &f.vector);
+
+ /* Progressively reduce the matrix until it is diagonal */
+
+ b = N - 1;
+
+ while (b > 0)
+ {
+ double fbm1 = gsl_vector_get (&f.vector, b - 1);
+
+ if (fbm1 == 0.0 || gsl_isnan (fbm1))
+ {
+ b--;
+ continue;
+ }
+
+ /* Find the largest unreduced block (a,b) starting from b
+ and working backwards */
+
+ a = b - 1;
+
+ while (a > 0)
+ {
+ double fam1 = gsl_vector_get (&f.vector, a - 1);
+
+ if (fam1 == 0.0 || gsl_isnan (fam1))
+ {
+ break;
+ }
+
+ a--;
+ }
+
+ {
+ const size_t n_block = b - a + 1;
+ gsl_vector_view S_block = gsl_vector_subvector (S, a, n_block);
+ gsl_vector_view f_block = gsl_vector_subvector (&f.vector, a, n_block - 1);
+
+ gsl_matrix_view U_block =
+ gsl_matrix_submatrix (A, 0, a, A->size1, n_block);
+ gsl_matrix_view V_block =
+ gsl_matrix_submatrix (V, 0, a, V->size1, n_block);
+
+ qrstep (&S_block.vector, &f_block.vector, &U_block.matrix, &V_block.matrix);
+
+ /* remove any small off-diagonal elements */
+
+ chop_small_elements (&S_block.vector, &f_block.vector);
+ }
+ }
+ }
+ /* Make singular values positive by reflections if necessary */
+
+ for (j = 0; j < K; j++)
+ {
+ double Sj = gsl_vector_get (S, j);
+
+ if (Sj < 0.0)
+ {
+ for (i = 0; i < N; i++)
+ {
+ double Vij = gsl_matrix_get (V, i, j);
+ gsl_matrix_set (V, i, j, -Vij);
+ }
+
+ gsl_vector_set (S, j, -Sj);
+ }
+ }
+
+ /* Sort singular values into decreasing order */
+
+ for (i = 0; i < K; i++)
+ {
+ double S_max = gsl_vector_get (S, i);
+ size_t i_max = i;
+
+ for (j = i + 1; j < K; j++)
+ {
+ double Sj = gsl_vector_get (S, j);
+
+ if (Sj > S_max)
+ {
+ S_max = Sj;
+ i_max = j;
+ }
+ }
+
+ if (i_max != i)
+ {
+ /* swap eigenvalues */
+ gsl_vector_swap_elements (S, i, i_max);
+
+ /* swap eigenvectors */
+ gsl_matrix_swap_columns (A, i, i_max);
+ gsl_matrix_swap_columns (V, i, i_max);
+ }
+ }
+
+ return GSL_SUCCESS;
+}
+
+
+/* Modified algorithm which is better for M>>N */
+
+int
+gsl_linalg_SV_decomp_mod (gsl_matrix * A,
+ gsl_matrix * X,
+ gsl_matrix * V, gsl_vector * S, gsl_vector * work)
+{
+ size_t i, j;
+
+ const size_t M = A->size1;
+ const size_t N = A->size2;
+
+ if (M < N)
+ {
+ GSL_ERROR ("svd of MxN matrix, M<N, is not implemented", GSL_EUNIMPL);
+ }
+ else if (V->size1 != N)
+ {
+ GSL_ERROR ("square matrix V must match second dimension of matrix A",
+ GSL_EBADLEN);
+ }
+ else if (V->size1 != V->size2)
+ {
+ GSL_ERROR ("matrix V must be square", GSL_ENOTSQR);
+ }
+ else if (X->size1 != N)
+ {
+ GSL_ERROR ("square matrix X must match second dimension of matrix A",
+ GSL_EBADLEN);
+ }
+ else if (X->size1 != X->size2)
+ {
+ GSL_ERROR ("matrix X must be square", GSL_ENOTSQR);
+ }
+ else if (S->size != N)
+ {
+ GSL_ERROR ("length of vector S must match second dimension of matrix A",
+ GSL_EBADLEN);
+ }
+ else if (work->size != N)
+ {
+ GSL_ERROR ("length of workspace must match second dimension of matrix A",
+ GSL_EBADLEN);
+ }
+
+ if (N == 1)
+ {
+ gsl_vector_view column = gsl_matrix_column (A, 0);
+ double norm = gsl_blas_dnrm2 (&column.vector);
+
+ gsl_vector_set (S, 0, norm);
+ gsl_matrix_set (V, 0, 0, 1.0);
+
+ if (norm != 0.0)
+ {
+ gsl_blas_dscal (1.0/norm, &column.vector);
+ }
+
+ return GSL_SUCCESS;
+ }
+
+ /* Convert A into an upper triangular matrix R */
+
+ for (i = 0; i < N; i++)
+ {
+ gsl_vector_view c = gsl_matrix_column (A, i);
+ gsl_vector_view v = gsl_vector_subvector (&c.vector, i, M - i);
+ double tau_i = gsl_linalg_householder_transform (&v.vector);
+
+ /* Apply the transformation to the remaining columns */
+
+ if (i + 1 < N)
+ {
+ gsl_matrix_view m =
+ gsl_matrix_submatrix (A, i, i + 1, M - i, N - (i + 1));
+ gsl_linalg_householder_hm (tau_i, &v.vector, &m.matrix);
+ }
+
+ gsl_vector_set (S, i, tau_i);
+ }
+
+ /* Copy the upper triangular part of A into X */
+
+ for (i = 0; i < N; i++)
+ {
+ for (j = 0; j < i; j++)
+ {
+ gsl_matrix_set (X, i, j, 0.0);
+ }
+
+ {
+ double Aii = gsl_matrix_get (A, i, i);
+ gsl_matrix_set (X, i, i, Aii);
+ }
+
+ for (j = i + 1; j < N; j++)
+ {
+ double Aij = gsl_matrix_get (A, i, j);
+ gsl_matrix_set (X, i, j, Aij);
+ }
+ }
+
+ /* Convert A into an orthogonal matrix L */
+
+ for (j = N; j > 0 && j--;)
+ {
+ /* Householder column transformation to accumulate L */
+ double tj = gsl_vector_get (S, j);
+ gsl_matrix_view m = gsl_matrix_submatrix (A, j, j, M - j, N - j);
+ gsl_linalg_householder_hm1 (tj, &m.matrix);
+ }
+
+ /* unpack R into X V S */
+
+ gsl_linalg_SV_decomp (X, V, S, work);
+
+ /* Multiply L by X, to obtain U = L X, stored in U */
+
+ {
+ gsl_vector_view sum = gsl_vector_subvector (work, 0, N);
+
+ for (i = 0; i < M; i++)
+ {
+ gsl_vector_view L_i = gsl_matrix_row (A, i);
+ gsl_vector_set_zero (&sum.vector);
+
+ for (j = 0; j < N; j++)
+ {
+ double Lij = gsl_vector_get (&L_i.vector, j);
+ gsl_vector_view X_j = gsl_matrix_row (X, j);
+ gsl_blas_daxpy (Lij, &X_j.vector, &sum.vector);
+ }
+
+ gsl_vector_memcpy (&L_i.vector, &sum.vector);
+ }
+ }
+
+ return GSL_SUCCESS;
+}
+
+
+/* Solves the system A x = b using the SVD factorization
+ *
+ * A = U S V^T
+ *
+ * to obtain x. For M x N systems it finds the solution in the least
+ * squares sense.
+ */
+
+int
+gsl_linalg_SV_solve (const gsl_matrix * U,
+ const gsl_matrix * V,
+ const gsl_vector * S,
+ const gsl_vector * b, gsl_vector * x)
+{
+ if (U->size1 != b->size)
+ {
+ GSL_ERROR ("first dimension of matrix U must size of vector b",
+ GSL_EBADLEN);
+ }
+ else if (U->size2 != S->size)
+ {
+ GSL_ERROR ("length of vector S must match second dimension of matrix U",
+ GSL_EBADLEN);
+ }
+ else if (V->size1 != V->size2)
+ {
+ GSL_ERROR ("matrix V must be square", GSL_ENOTSQR);
+ }
+ else if (S->size != V->size1)
+ {
+ GSL_ERROR ("length of vector S must match size of matrix V",
+ GSL_EBADLEN);
+ }
+ else if (V->size2 != x->size)
+ {
+ GSL_ERROR ("size of matrix V must match size of vector x", GSL_EBADLEN);
+ }
+ else
+ {
+ const size_t N = U->size2;
+ size_t i;
+
+ gsl_vector *w = gsl_vector_calloc (N);
+
+ gsl_blas_dgemv (CblasTrans, 1.0, U, b, 0.0, w);
+
+ for (i = 0; i < N; i++)
+ {
+ double wi = gsl_vector_get (w, i);
+ double alpha = gsl_vector_get (S, i);
+ if (alpha != 0)
+ alpha = 1.0 / alpha;
+ gsl_vector_set (w, i, alpha * wi);
+ }
+
+ gsl_blas_dgemv (CblasNoTrans, 1.0, V, w, 0.0, x);
+
+ gsl_vector_free (w);
+
+ return GSL_SUCCESS;
+ }
+}
+
+/* This is a the jacobi version */
+/* Author: G. Jungman */
+
+/*
+ * Algorithm due to J.C. Nash, Compact Numerical Methods for
+ * Computers (New York: Wiley and Sons, 1979), chapter 3.
+ * See also Algorithm 4.1 in
+ * James Demmel, Kresimir Veselic, "Jacobi's Method is more
+ * accurate than QR", Lapack Working Note 15 (LAWN15), October 1989.
+ * Available from netlib.
+ *
+ * Based on code by Arthur Kosowsky, Rutgers University
+ * kosowsky@physics.rutgers.edu
+ *
+ * Another relevant paper is, P.P.M. De Rijk, "A One-Sided Jacobi
+ * Algorithm for computing the singular value decomposition on a
+ * vector computer", SIAM Journal of Scientific and Statistical
+ * Computing, Vol 10, No 2, pp 359-371, March 1989.
+ *
+ */
+
+int
+gsl_linalg_SV_decomp_jacobi (gsl_matrix * A, gsl_matrix * Q, gsl_vector * S)
+{
+ if (A->size1 < A->size2)
+ {
+ /* FIXME: only implemented M>=N case so far */
+
+ GSL_ERROR ("svd of MxN matrix, M<N, is not implemented", GSL_EUNIMPL);
+ }
+ else if (Q->size1 != A->size2)
+ {
+ GSL_ERROR ("square matrix Q must match second dimension of matrix A",
+ GSL_EBADLEN);
+ }
+ else if (Q->size1 != Q->size2)
+ {
+ GSL_ERROR ("matrix Q must be square", GSL_ENOTSQR);
+ }
+ else if (S->size != A->size2)
+ {
+ GSL_ERROR ("length of vector S must match second dimension of matrix A",
+ GSL_EBADLEN);
+ }
+ else
+ {
+ const size_t M = A->size1;
+ const size_t N = A->size2;
+ size_t i, j, k;
+
+ /* Initialize the rotation counter and the sweep counter. */
+ int count = 1;
+ int sweep = 0;
+ int sweepmax = 5*N;
+
+ double tolerance = 10 * M * GSL_DBL_EPSILON;
+
+ /* Always do at least 12 sweeps. */
+ sweepmax = GSL_MAX (sweepmax, 12);
+
+ /* Set Q to the identity matrix. */
+ gsl_matrix_set_identity (Q);
+
+ /* Store the column error estimates in S, for use during the
+ orthogonalization */
+
+ for (j = 0; j < N; j++)
+ {
+ gsl_vector_view cj = gsl_matrix_column (A, j);
+ double sj = gsl_blas_dnrm2 (&cj.vector);
+ gsl_vector_set(S, j, GSL_DBL_EPSILON * sj);
+ }
+
+ /* Orthogonalize A by plane rotations. */
+
+ while (count > 0 && sweep <= sweepmax)
+ {
+ /* Initialize rotation counter. */
+ count = N * (N - 1) / 2;
+
+ for (j = 0; j < N - 1; j++)
+ {
+ for (k = j + 1; k < N; k++)
+ {
+ double a = 0.0;
+ double b = 0.0;
+ double p = 0.0;
+ double q = 0.0;
+ double cosine, sine;
+ double v;
+ double abserr_a, abserr_b;
+ int sorted, orthog, noisya, noisyb;
+
+ gsl_vector_view cj = gsl_matrix_column (A, j);
+ gsl_vector_view ck = gsl_matrix_column (A, k);
+
+ gsl_blas_ddot (&cj.vector, &ck.vector, &p);
+ p *= 2.0 ; /* equation 9a: p = 2 x.y */
+
+ a = gsl_blas_dnrm2 (&cj.vector);
+ b = gsl_blas_dnrm2 (&ck.vector);
+
+ q = a * a - b * b;
+ v = hypot(p, q);
+
+ /* test for columns j,k orthogonal, or dominant errors */
+
+ abserr_a = gsl_vector_get(S,j);
+ abserr_b = gsl_vector_get(S,k);
+
+ sorted = (gsl_coerce_double(a) >= gsl_coerce_double(b));
+ orthog = (fabs (p) <= tolerance * gsl_coerce_double(a * b));
+ noisya = (a < abserr_a);
+ noisyb = (b < abserr_b);
+
+ if (sorted && (orthog || noisya || noisyb))
+ {
+ count--;
+ continue;
+ }
+
+ /* calculate rotation angles */
+ if (v == 0 || !sorted)
+ {
+ cosine = 0.0;
+ sine = 1.0;
+ }
+ else
+ {
+ cosine = sqrt((v + q) / (2.0 * v));
+ sine = p / (2.0 * v * cosine);
+ }
+
+ /* apply rotation to A */
+ for (i = 0; i < M; i++)
+ {
+ const double Aik = gsl_matrix_get (A, i, k);
+ const double Aij = gsl_matrix_get (A, i, j);
+ gsl_matrix_set (A, i, j, Aij * cosine + Aik * sine);
+ gsl_matrix_set (A, i, k, -Aij * sine + Aik * cosine);
+ }
+
+ gsl_vector_set(S, j, fabs(cosine) * abserr_a + fabs(sine) * abserr_b);
+ gsl_vector_set(S, k, fabs(sine) * abserr_a + fabs(cosine) * abserr_b);
+
+ /* apply rotation to Q */
+ for (i = 0; i < N; i++)
+ {
+ const double Qij = gsl_matrix_get (Q, i, j);
+ const double Qik = gsl_matrix_get (Q, i, k);
+ gsl_matrix_set (Q, i, j, Qij * cosine + Qik * sine);
+ gsl_matrix_set (Q, i, k, -Qij * sine + Qik * cosine);
+ }
+ }
+ }
+
+ /* Sweep completed. */
+ sweep++;
+ }
+
+ /*
+ * Orthogonalization complete. Compute singular values.
+ */
+
+ {
+ double prev_norm = -1.0;
+
+ for (j = 0; j < N; j++)
+ {
+ gsl_vector_view column = gsl_matrix_column (A, j);
+ double norm = gsl_blas_dnrm2 (&column.vector);
+
+ /* Determine if singular value is zero, according to the
+ criteria used in the main loop above (i.e. comparison
+ with norm of previous column). */
+
+ if (norm == 0.0 || prev_norm == 0.0
+ || (j > 0 && norm <= tolerance * prev_norm))
+ {
+ gsl_vector_set (S, j, 0.0); /* singular */
+ gsl_vector_set_zero (&column.vector); /* annihilate column */
+
+ prev_norm = 0.0;
+ }
+ else
+ {
+ gsl_vector_set (S, j, norm); /* non-singular */
+ gsl_vector_scale (&column.vector, 1.0 / norm); /* normalize column */
+
+ prev_norm = norm;
+ }
+ }
+ }
+
+ if (count > 0)
+ {
+ /* reached sweep limit */
+ GSL_ERROR ("Jacobi iterations did not reach desired tolerance",
+ GSL_ETOL);
+ }
+
+ return GSL_SUCCESS;
+ }
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-30 21:10:10 +0000