eigen_internal.c

浏览该文件的文档.

 
#include <math.h>
#include <string.h>
#include "memory_bridge.h"
#include "internal.h"
#include "eigen_internal.h"
#include "lmmc/config.h"
#include "lmmc/dense.h"
#include "lmmc/eigen.h"
#include "lmmc/linear_algebra.h"
 
lmmc_status_t lmmc_eigen_matrix_scale(const lmmc_mat_t* a,
                                      lmmc_real_t* scale)
{
    *scale = 0.0;
    for (size_t i = 0; i < a->rows; ++i) {
        for (size_t j = 0; j < a->cols; ++j) {
            lmmc_real_t magnitude = fabs(MAT_ELEM(a, i, j));
            if (!isfinite(magnitude)) {
                return LMMC_STATUS_NUMERICAL_FAILURE;
            }
            if (magnitude > *scale) {
                *scale = magnitude;
            }
        }
    }
    if (*scale == 0.0) {
        *scale = 1.0;
    }
    return LMMC_STATUS_OK;
}
 
void householder_make(lmmc_real_t *x, size_t len, lmmc_real_t *tau_out,
                      lmmc_real_t *beta_out) {
    lmmc_scaled_sumsq_t acc;
    lmmc_real_t alpha;
    lmmc_real_t tail_norm;
    lmmc_real_t norm;
    lmmc_real_t beta;
    lmmc_real_t denominator_scaled;
 
    if (len == 0) { *tau_out = 0.0; *beta_out = 0.0; return; }
    if (len == 1) { *tau_out = 0.0; *beta_out = x[0]; return; }
 
    lmmc_scaled_sumsq_init(&acc);
    for (size_t i = 1; i < len; ++i) {
        lmmc_scaled_sumsq_add(&acc, x[i]);
    }
    tail_norm = lmmc_scaled_sumsq_norm(&acc);
    alpha = x[0];
    if (tail_norm == 0.0) {
        *tau_out = 0.0;
        *beta_out = alpha;
        return;
    }
 
    norm = hypot(alpha, tail_norm);
    beta = -copysign(norm, alpha);
    *beta_out = beta;
    *tau_out = 1.0 - alpha / beta;
    denominator_scaled = alpha / norm - beta / norm;
    for (size_t i = 1; i < len; ++i) {
        x[i] = (x[i] / norm) / denominator_scaled;
    }
    x[0] = 1.0;
}
 
void householder_apply_left(lmmc_mat_t *M, size_t i0, size_t len,
                                   size_t j0, size_t j1,
                                   const lmmc_real_t *v, lmmc_real_t tau) {
    size_t i, j;
    if (tau == 0.0) return;
    for (j = j0; j < j1; j++) {
        lmmc_real_t s = MAT_ELEM(M, i0, j);
        for (i = 1; i < len; i++) s += v[i] * MAT_ELEM(M, i0 + i, j);
        s *= tau;
        MAT_ELEM(M, i0, j) -= s;
        for (i = 1; i < len; i++) MAT_ELEM(M, i0 + i, j) -= s * v[i];
    }
}
 
void householder_apply_right(lmmc_mat_t *M, size_t i0, size_t i1,
                                    size_t j0, size_t len,
                                    const lmmc_real_t *v, lmmc_real_t tau) {
    size_t i, j;
    if (tau == 0.0) return;
    for (i = i0; i < i1; i++) {
        lmmc_real_t s = MAT_ELEM(M, i, j0);
        for (j = 1; j < len; j++) s += v[j] * MAT_ELEM(M, i, j0 + j);
        s *= tau;
        MAT_ELEM(M, i, j0) -= s;
        for (j = 1; j < len; j++) MAT_ELEM(M, i, j0 + j) -= s * v[j];
    }
}
 
void givens_right(lmmc_mat_t *M, size_t rows, size_t p, size_t q,
                         lmmc_real_t c, lmmc_real_t s) {
    size_t i;
    for (i = 0; i < rows; i++) {
        lmmc_real_t mp = MAT_ELEM(M, i, p);
        lmmc_real_t mq = MAT_ELEM(M, i, q);
        MAT_ELEM(M, i, p) =  c * mp + s * mq;
        MAT_ELEM(M, i, q) = -s * mp + c * mq;
    }
}
 
void givens_compute(lmmc_real_t a, lmmc_real_t b,
                           lmmc_real_t *c, lmmc_real_t *s, lmmc_real_t *r) {
    if (b == 0.0) { *c = (a >= 0.0) ? 1.0 : -1.0; *s = 0.0; *r = lmmc_abs(a); }
    else if (a == 0.0) { *c = 0.0; *s = (b >= 0.0) ? 1.0 : -1.0; *r = lmmc_abs(b); }
    else {
        lmmc_real_t h = hypot(a, b);
        *r = h;
        *c = a / h;
        *s = b / h;
    }
}