含有Rayleigh商加速的幂法求矩阵的特征值及特征向量代码c语言

以下是含有Rayleigh商加速的幂法求矩阵的特征值及特征向量的C语言代码：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAX_ITER 1000     // 最大迭代次数
#define EPSILON 1e-10     // 精度

void matrix_vector_multiply(double *matrix, double *vector, double *result, int n) {
    int i, j;
    for (i = 0; i < n; i++) {
        *(result + i) = 0;
        for (j = 0; j < n; j++) {
            *(result + i) += *(matrix + i * n + j) * *(vector + j);
        }
    }
}

double vector_norm(double *vector, int n) {
    int i;
    double norm = 0;
    for (i = 0; i < n; i++) {
        norm += *(vector + i) * *(vector + i);
    }
    return sqrt(norm);
}

void matrix_eigenvalue_eigenvector_power(double *matrix, double *eigenvalue, double *eigenvector, int n) {
    int i, j, k, iter = 0;
    double norm_v, lambda, mu, delta;
    double *v, *u, *temp;

    // 初始化特征向量
    for (i = 0; i < n; i++) {
        *(eigenvector + i) = 1;
    }

    while (iter < MAX_ITER) {
        norm_v = vector_norm(eigenvector, n);
        for (i = 0; i < n; i++) {
            *(eigenvector + i) /= norm_v;
        }
        v = (double *) malloc(n * sizeof(double));
        u = (double *) malloc(n * sizeof(double));
        temp = (double *) malloc(n * sizeof(double));
        matrix_vector_multiply(matrix, eigenvector, v, n);
        mu = vector_norm(v, n);
        lambda = 0;
        for (i = 0; i < n; i++) {
            lambda += *(eigenvector + i) * *(v + i);
        }
        lambda /= norm_v * norm_v;
        for (k = 0; k < 10; k++) {
            for (i = 0; i < n; i++) {
                *(u + i) = *(v + i) - lambda * *(eigenvector + i);
            }
            norm_v = vector_norm(u, n);
            for (i = 0; i < n; i++) {
                *(u + i) /= norm_v;
            }
            matrix_vector_multiply(matrix, u, temp, n);
            mu = vector_norm(temp, n);
            delta = 0;
            for (i = 0; i < n; i++) {
                delta += *(u + i) * *(temp + i);
            }
            delta /= norm_v * norm_v;
            if (fabs(lambda - delta) < EPSILON) {
                lambda = delta;
                break;
            } else {
                lambda = delta;
            }
        }
        for (i = 0; i < n; i++) {
            *(eigenvector + i) = *(u + i);
        }
        *(eigenvalue + iter) = lambda;
        iter++;
        free(v);
        free(u);
        free(temp);
    }
}

int main() {
    int n, i, j;
    double *matrix, *eigenvalue, *eigenvector;

    printf("请输入矩阵的维数：");
    scanf("%d", &n);

    // 申请内存
    matrix = (double *) malloc(n * n * sizeof(double));
    eigenvalue = (double *) malloc(MAX_ITER * sizeof(double));
    eigenvector = (double *) malloc(n * sizeof(double));

    printf("请输入矩阵的元素：\n");
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            scanf("%lf", matrix + i * n + j);
        }
    }

    // 求解特征值和特征向量
    matrix_eigenvalue_eigenvector_power(matrix, eigenvalue, eigenvector, n);

    printf("特征值为：\n");
    for (i = 0; i < MAX_ITER; i++) {
        printf("%lf ", *(eigenvalue + i));
    }
    printf("\n");

    printf("特征向量为：\n");
    for (i = 0; i < n; i++) {
        for (j = 0; j < MAX_ITER; j++) {
            printf("%lf ", *(eigenvector + i) * pow(*(eigenvalue + j), j));
        }
        printf("\n");
    }

    // 释放内存
    free(matrix);
    free(eigenvalue);
    free(eigenvector);

    return 0;
}