1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
|
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_linalg.h>
#include "gsl_functions_pinv.h"
int
gsl_elinalg_complex_QR_decomp (gsl_matrix * A, gsl_vector * tau)
{
const size_t M = A->size1;
const size_t N = A->size2;
if (tau->size != GSL_MIN (M, N))
{
GSL_ERROR ("size of tau must be MIN(M,N)", GSL_EBADLEN);
}
else
{
size_t i;
for (i = 0; i < GSL_MIN (M, N); i++)
{
/* Compute the Householder transformation to reduce the j-th
column of the matrix to a multiple of the j-th unit vector */
gsl_vector_view c_full = gsl_matrix_column (A, i);
gsl_vector_view c = gsl_vector_subvector (&(c_full.vector), i, M-i);
double tau_i = gsl_linalg_householder_transform (&(c.vector));
gsl_vector_set (tau, i, tau_i);
/* Apply the transformation to the remaining columns and
update the norms */
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, &(c.vector), &(m.matrix));
}
}
return GSL_SUCCESS;
}
}
[/quote] |
Partager