简介这是一份C语言项目源码将复数矩阵特征值计算与黑白棋搜索决策整合在一个.c文件中面向数值算法和C语言实战学习者。代码实现了自定义复数结构体、矩阵结构体和复数矩阵乘法、加法等基本运算并通过幂迭代法近似求解特征值同时包含黑白棋规则判断、深度搜索及边角权重评估逻辑清晰展现了从数学理论到编程实现的转化过程。压缩包内仅1个c文件大小约7KB结构紧凑便于逐个函数阅读。已有139人学习浏览在矩阵计算与游戏AI结合的小型项目中有一定参考价值适合课程设计或算法实验场景。通过阅读源码可以同时掌握复数矩阵运算的实现细节和黑白棋评估函数的权重设定思路是一个较完整的数值计算加博弈算法样例。1. 复数矩阵特征值C语言项目真正要攻克的三个难点复数矩阵特征值在数值线性代数里是一个可以被一句话定义、却很难在三十分钟内写完的问题。对C语言来说尤其如此没有numpy.linalg.eig可用没有MATLAB的eig函数连复数类型也要先在标准库和自建结构之间做取舍。以“final”命名的这类源码常见于课程设计与算法原型需求通常是把一个n阶复矩阵的所有特征值算出来精度要求到1e-8还要交付一个完整可编译的C语言项目。这个题目真正考察的不是背公式而是三件事复数矩阵怎么存、QR迭代怎么收敛、结果怎么验证。把这三件事分开处理整个源码的模块边界就清晰了后面的编码和排错才有抓手。2. 复数与矩阵运算基础复数矩阵特征值c语言源码的地基2.1 用C99 complex.h还是自建复数结构体C语言项目处理复数的第一道选择题是类型方案。C99标准提供了double complex内建类型和complex.h头文件支持 - * /四则运算以及conj、cabs、creal、cimag等函数。如果编译环境是GCC或Clang直接用double complex最省事代码里每个复数运算都和数学表达式一一对应不容易写错。自建结构体的典型写法是typedef struct { double re; double im; } MyComplex;好处是不依赖C99老版本MSVC也能编译坏处是所有四则运算都要手写代码量会成倍增加。这个项目选C99路线但会把复数类型统一改成别名方便日后整体替换。#include stdio.h #include stdlib.h #include math.h #include complex.h typedef double complex cx; #define IDX(i, j, n) ((i) * (n) (j))IDX宏是行主序矩阵的下标缩写后面所有源码都靠它降低噪音。复数类型别名cx让函数签名短一截改动复数表示方式时也只动这一行。2.2 行主序存储与矩阵内存管理复数矩阵按一维数组存行主序A[i * n j]是第i行第j列。n阶矩阵需要连续n * n个复数单元每个double complex通常占16字节。分配时用calloc而不是malloc因为calloc会把内存清零避免把未初始化的数据带进计算malloc不会清零第一次打印矩阵时看到随机虚部多半就是这个原因。C语言项目里内存管理是评分重灾区。复数矩阵特征值计算的每一次分配都要有对应的free否则Valgrind一跑final代码直接扣印象分。常见做法是把分配和释放封装在同一层函数里保证成对出现。static cx *mat_new(int n) { return (cx *)calloc((size_t)n * n, sizeof(cx)); } static void mat_free(cx *A) { free(A); }注意calloc参数顺序是数量、大小别写反n为0时返回什么由实现决定调用侧最好在n较小的时候直接判空退出。矩阵尺寸超过几万时n * n * sizeof(cx)可能溢出size_t这在课设规模下遇不到但做通用库时要在入口处检查。2.3 乘法、共轭转置与F范数三个必写的矩阵工具特征值计算绕不开三个基础操作复矩阵乘法、共轭转置、Frobenius范数。复矩阵乘法和实数写法在结构上完全一样但每个乘加都是复数乘法编译器会链接复数运算支持速度比实矩阵慢这是正常的。共轭转置在厄米矩阵判定和Householder变换里都会用到公式是(A^H)[j][i] conj(A[i][j])少了conj就是普通转置后续算法会全部走偏。static void mat_mul(cx *C, const cx *A, const cx *B, int n) { for (int i 0; i n; i) for (int j 0; j n; j) { C[IDX(i, j, n)] 0; for (int k 0; k n; k) C[IDX(i, j, n)] A[IDX(i, k, n)] * B[IDX(k, j, n)]; } } static void mat_conj_transpose(cx *At, const cx *A, int n) { for (int i 0; i n; i) for (int j 0; j n; j) At[IDX(j, i, n)] conj(A[IDX(i, j, n)]); } static double mat_fnorm(const cx *A, int n) { double s 0.0; for (int i 0; i n; i) for (int j 0; j n; j) s creal(A[IDX(i, j, n)] * conj(A[IDX(i, j, n)])); return sqrt(s); }mat_mul要求C、A、B三块内存互不重叠否则循环中会读到已经被覆盖的旧值QR迭代里我会用独立矩阵暂存来规避。mat_fnorm返回矩阵的F范数它在收敛判断中用来衡量下三角残余作用比名字看起来重要得多。这三个函数建议单独放到cmat.c里导出头文件后文所有算法只需要包含接口。3. 特征值算法选型为什么这个C语言项目绕不开QR迭代3.1 特征多项式不是正路数值稳定性的坏消息一看到特征值就想到det(A - λI) 0这是线性代数课堂路线但不是数值计算路线。Abel-Ruffini定理决定了五次以上多项式没有通用根式解这是理论层面的死路即使对低阶矩阵强行展开特征多项式得到的系数也是矩阵元素的高次组合微小浮点误差会被放大求根结果常出现完全虚假的复数共轭对。如果final源码走特征多项式路线n4以上开始翻车n8以上基本不可用。原因很直接展开行列式本身就是大量乘加运算浮点误差在过程中累积多项式求根又依赖系数系数误差被放大后根的位置完全失控。QR迭代不构造特征多项式而是通过正交相似变换逐步把矩阵化成上三角对角线直接就是特征值这条路数值上稳定得多也是LAPACK、GSL等库的实际做法。3.2 厄米矩阵可走Jacobi一般复矩阵走QR迭代写代码前先判断矩阵类型。若A^H A即共轭转置等于自身称为厄米矩阵其特征值全为实数自研代码可以用复数Jacobi旋转每次消一个非对角元实现简单。但复数矩阵特征值的一般情况没有厄米约束Jacobi对非厄米矩阵不收敛必须上QR迭代。LAPACK的zgeev就是对一般复矩阵做原位QR迭代并叠加了Hessenberg化、移位和平衡化。先看矩阵有没有结构再看是否需要自己造轮子矩阵类型推荐方法特征值类型实现成本厄米矩阵复数Jacobi旋转实数低一般复矩阵移位QR迭代复数中高任意矩阵工程快速交付LAPACKE zgeev复数极低为什么厄米矩阵特征值一定是实数对特征对Ax λx两边左乘x^H得到x^H A x λ x^H x因为A A^Hx^H A x的共轭等于自身所以λ等于自己的共轭虚部只能为0。这个结论能用来做验证如果输入是厄米矩阵输出却出现明显虚部说明QR迭代代码有bug。3.3 移位QR迭代的收敛逻辑与停止条件QR迭代的基本流程是三步对当前矩阵A_k做QR分解得到酉矩阵Q_k和上三角R_k计算A_{k1} R_k Q_k重复直到下三角元素小到可忽略。因为A_{k1} Q_k^H A_k Q_k每轮都是相似变换特征值集合同原来完全一样。当k足够大时A_k的下三角元素趋向0对角线趋向特征值。收敛速度取决于相邻特征值的模之比模相近时很慢所以实用代码都要加移位每轮从A_k右下角减一个数μ做QR分解得到R Q后再加回μ。复矩阵的μ也可以是复数后面调参章节再展开。停止条件用一个绝对阈值tol。计算下三角元素的模方和开方后小于tol就停止double off 0.0; for (int i 0; i n; i) for (int j 0; j i; j) off creal(A[IDX(i, j, n)] * conj(A[IDX(i, j, n)])); off sqrt(off);阈值取太大会得到错误特征值取太小小则迭代次数爆炸。我一般从1e-10起步高精度需求下放宽到1e-12同时配合最大迭代次数双重保护。4. 源码实现跑通复数矩阵特征值的最终代码4.1 工程最快方案用LAPACKE封装zgeev如果final项目允许链接外部库直接调用LAPACKE_zgeev最省事。这是LAPACK对一般复矩阵zgeev的C接口内部完成了平衡化、Hessenberg化、移位QR以及特征向量的可选计算数值稳定性远好于课程作业级自研代码。调用代码很短#include stdio.h #include complex.h #include lapacke.h void eig_lapack(int n, double complex *A, double complex *w) { int info LAPACKE_zgeev(LAPACK_ROW_MAJOR, N, N, n, A, n, w, NULL, 1, NULL, 1); if (info ! 0) fprintf(stderr, zgeev failed, info%d\n, info); }LAPACK_ROW_MAJOR表示调用方按行主序传入矩阵和前面的IDX布局一致两个N表示不计算左右特征向量A会被内部工作空间覆盖需要保留原矩阵就先拷贝一份w是长度为n的复数数组接收特征值。编译时链接-llapacke -llapack -lblas -lm系统需要装liblapacke-dev。缺点是部分环境没有预装LAPACK交叉编译到嵌入式设备时更麻烦这时自研QR迭代就成了实际可行的备选。4.2 自主实现复数Householder QR分解与QR迭代主循环自己写QR迭代复杂度集中在QR分解。这里用Householder变换实现复数版QR分解对第k列取出子列x构造反射向量v使Hx除第一个元素外都变成0。复数版与实数版的差别在两点一是反射向量要考虑x[0]的相位alpha -(x[0] / |x[0]|) * ||x||这样算出的v模最大数值最稳二是更新矩阵时用共轭转置v^H而不是普通转置。static void qr_decomp(cx *A, cx *Q, int n) { for (int i 0; i n; i) for (int j 0; j n; j) Q[IDX(i, j, n)] (i j) ? 1.0 : 0.0; for (int k 0; k n - 1; k) { int m n - k; cx *x (cx *)malloc(sizeof(cx) * m); for (int i 0; i m; i) x[i] A[IDX(k i, k, n)]; double nrm 0.0; for (int i 0; i m; i) nrm creal(x[i] * conj(x[i])); nrm sqrt(nrm); if (nrm 1e-300) { free(x); continue; } cx alpha (cabs(x[0]) 1e-300) ? -nrm : -(x[0] / cabs(x[0])) * nrm; cx *v (cx *)malloc(sizeof(cx) * m); v[0] x[0] - alpha; for (int i 1; i m; i) v[i] x[i]; double vh_v 0.0; for (int i 0; i m; i) vh_v creal(v[i] * conj(v[i])); if (vh_v 1e-300) { free(x); free(v); continue; } double scale 2.0 / vh_v; // A H A A - scale * v * (v^H A) for (int j k; j n; j) { cx dot 0.0; for (int l 0; l m; l) dot conj(v[l]) * A[IDX(k l, j, n)]; dot * scale; for (int i 0; i m; i) A[IDX(k i, j, n)] - v[i] * dot; } // Q Q H Q - scale * (Q v) * v^H cx *u (cx *)malloc(sizeof(cx) * n); for (int i 0; i n; i) { cx s 0.0; for (int l 0; l m; l) s Q[IDX(i, k l, n)] * v[l]; u[i] s * scale; } for (int i 0; i n; i) for (int j 0; j m; j) Q[IDX(i, k j, n)] - u[i] * conj(v[j]); free(x); free(v); free(u); } } void eig_qr(cx *A, cx *w, int n, int max_iter, double tol) { cx *Q mat_new(n); cx *R mat_new(n); int it; for (it 0; it max_iter; it) { qr_decomp(A, Q, n); // A_new R * Q用独立矩阵 R 暂存再拷回 for (int i 0; i n; i) for (int j 0; j n; j) { cx s 0.0; for (int k 0; k n; k) s A[IDX(i, k, n)] * Q[IDX(k, j, n)]; R[IDX(i, j, n)] s; } for (int i 0; i n; i) for (int j 0; j n; j) A[IDX(i, j, n)] R[IDX(i, j, n)]; double off 0.0; for (int i 0; i n; i) for (int j 0; j i; j) off creal(A[IDX(i, j, n)] * conj(A[IDX(i, j, n)])); if (sqrt(off) tol) break; } for (int i 0; i n; i) w[i] A[IDX(i, i, n)]; mat_free(Q); mat_free(R); }qr_decomp把A原地覆盖成R同时生成Q。更新A与Q的顺序可以互换因为二者都只依赖旧A和旧Q不存在相互覆盖。eig_qr中R * Q用独立矩阵暂存是必须的A此时就是R直接在A上做A A * Q会覆盖尚未读取的元素。收敛后A接近上三角对角线就是特征值。这个实现是无移位版本矩阵条件数好时能收敛遇到收敛慢的矩阵把max_iter调大或者按下一章的加移位方式改造。4.3 编译命令与链接参数把上述函数和一个main放在eigen.c里编译命令是gcc -stdc99 -O2 -Wall -o eigen eigen.c -lm-lm不能省复数运算中的conj、cabs、sqrt都在libm里。-O2对数值代码影响很大不开优化时浮点运算顺序可能不同特征值结果会在最后几位抖动提交final项目前务必用优化编译重测一遍。如果加了LAPACKE分支编译命令变成gcc -stdc99 -O2 -o eigen eigen.c -llapacke -llapack -lblas -lm链接顺序要按依赖关系排-llapacke在最前因为它依赖-llapack和-lblas。报错出现undefined reference to zgeev_多半是lapack没装或链接顺序写反了。5. 验证与调参final版复数矩阵特征值C语言源码的收尾检查5.1 用已知复特征值的矩阵验证输出验证是final项目最该做却最容易被省掉的一步。用一个已知特征值的实矩阵[[2, 1], [-1, 2]]其特征多项式是λ^2 - 4λ 5特征值应为2 i与2 - i。测试代码cx A[4] {2.0, 1.0, -1.0, 2.0}; cx w[2]; eig_qr(A, w, 2, 1000, 1e-10); for (int i 0; i 2; i) printf(%.6f %.6f i\n, creal(w[i]), cimag(w[i]));理想输出是2.000000 1.000000i和2.000000 -1.000000i。如果虚部符号相反只是特征值顺序不同这是正常的特征值本身没有固定顺序比较前先排序。如果虚部完全消失大概率是QR分解更新A时把conj(v[l])写成了v[l]这是复数Householder变换最经典的错误。5.2 收敛公差、最大迭代次数与移位策略参数建议初值后果与调试方向tol1e-10偏大时特征值误差明显偏小则迭代次数多余max_iter500 * n达到上限未break先检查是否忘记加移位shiftA[n-1][n-1]加速收敛的关键复矩阵用复μ矩阵规模n 2..50n大于50应先用Householder化为Hessenberg形无移位QR对模相近的特征值收敛很慢。常见做法是每轮取μ A[n-1][n-1]作为移位先算A A - μI对A做QR更新A_new R * Q μI。右下角元素在迭代中逐步接近某个特征值这个移位能显著加快收敛复矩阵就取复μ。收敛后对角线就是特征值不需要再减μ。5.3 三个容易被忽略的边界问题特征值顺序不稳定。同一矩阵在不同n或不同轮数下输出顺序可能变化做对比测试时按实部排序实部相近再按虚部排序。矩阵尺度差异过大。元素跨10个数量级时直接QR可能因上溢得到NaN先做平衡化行列同时缩放LAPACK的zgebal专门做这件事。子列范数恰好为0。qr_decomp遇到nrm 1e-300直接跳过该列这对奇异或稀疏结构是合法的但主循环收敛判据仍会统计该列下三角残余迭代次数可能变多不应误判为死循环。把平衡化加进final项目通常能让特征值精度回到1e-10量级。本文还有配套的精品资源点击获取
复数矩阵特征值C语言实现：QR迭代与LAPACK实战指南
简介这是一份C语言项目源码将复数矩阵特征值计算与黑白棋搜索决策整合在一个.c文件中面向数值算法和C语言实战学习者。代码实现了自定义复数结构体、矩阵结构体和复数矩阵乘法、加法等基本运算并通过幂迭代法近似求解特征值同时包含黑白棋规则判断、深度搜索及边角权重评估逻辑清晰展现了从数学理论到编程实现的转化过程。压缩包内仅1个c文件大小约7KB结构紧凑便于逐个函数阅读。已有139人学习浏览在矩阵计算与游戏AI结合的小型项目中有一定参考价值适合课程设计或算法实验场景。通过阅读源码可以同时掌握复数矩阵运算的实现细节和黑白棋评估函数的权重设定思路是一个较完整的数值计算加博弈算法样例。1. 复数矩阵特征值C语言项目真正要攻克的三个难点复数矩阵特征值在数值线性代数里是一个可以被一句话定义、却很难在三十分钟内写完的问题。对C语言来说尤其如此没有numpy.linalg.eig可用没有MATLAB的eig函数连复数类型也要先在标准库和自建结构之间做取舍。以“final”命名的这类源码常见于课程设计与算法原型需求通常是把一个n阶复矩阵的所有特征值算出来精度要求到1e-8还要交付一个完整可编译的C语言项目。这个题目真正考察的不是背公式而是三件事复数矩阵怎么存、QR迭代怎么收敛、结果怎么验证。把这三件事分开处理整个源码的模块边界就清晰了后面的编码和排错才有抓手。2. 复数与矩阵运算基础复数矩阵特征值c语言源码的地基2.1 用C99 complex.h还是自建复数结构体C语言项目处理复数的第一道选择题是类型方案。C99标准提供了double complex内建类型和complex.h头文件支持 - * /四则运算以及conj、cabs、creal、cimag等函数。如果编译环境是GCC或Clang直接用double complex最省事代码里每个复数运算都和数学表达式一一对应不容易写错。自建结构体的典型写法是typedef struct { double re; double im; } MyComplex;好处是不依赖C99老版本MSVC也能编译坏处是所有四则运算都要手写代码量会成倍增加。这个项目选C99路线但会把复数类型统一改成别名方便日后整体替换。#include stdio.h #include stdlib.h #include math.h #include complex.h typedef double complex cx; #define IDX(i, j, n) ((i) * (n) (j))IDX宏是行主序矩阵的下标缩写后面所有源码都靠它降低噪音。复数类型别名cx让函数签名短一截改动复数表示方式时也只动这一行。2.2 行主序存储与矩阵内存管理复数矩阵按一维数组存行主序A[i * n j]是第i行第j列。n阶矩阵需要连续n * n个复数单元每个double complex通常占16字节。分配时用calloc而不是malloc因为calloc会把内存清零避免把未初始化的数据带进计算malloc不会清零第一次打印矩阵时看到随机虚部多半就是这个原因。C语言项目里内存管理是评分重灾区。复数矩阵特征值计算的每一次分配都要有对应的free否则Valgrind一跑final代码直接扣印象分。常见做法是把分配和释放封装在同一层函数里保证成对出现。static cx *mat_new(int n) { return (cx *)calloc((size_t)n * n, sizeof(cx)); } static void mat_free(cx *A) { free(A); }注意calloc参数顺序是数量、大小别写反n为0时返回什么由实现决定调用侧最好在n较小的时候直接判空退出。矩阵尺寸超过几万时n * n * sizeof(cx)可能溢出size_t这在课设规模下遇不到但做通用库时要在入口处检查。2.3 乘法、共轭转置与F范数三个必写的矩阵工具特征值计算绕不开三个基础操作复矩阵乘法、共轭转置、Frobenius范数。复矩阵乘法和实数写法在结构上完全一样但每个乘加都是复数乘法编译器会链接复数运算支持速度比实矩阵慢这是正常的。共轭转置在厄米矩阵判定和Householder变换里都会用到公式是(A^H)[j][i] conj(A[i][j])少了conj就是普通转置后续算法会全部走偏。static void mat_mul(cx *C, const cx *A, const cx *B, int n) { for (int i 0; i n; i) for (int j 0; j n; j) { C[IDX(i, j, n)] 0; for (int k 0; k n; k) C[IDX(i, j, n)] A[IDX(i, k, n)] * B[IDX(k, j, n)]; } } static void mat_conj_transpose(cx *At, const cx *A, int n) { for (int i 0; i n; i) for (int j 0; j n; j) At[IDX(j, i, n)] conj(A[IDX(i, j, n)]); } static double mat_fnorm(const cx *A, int n) { double s 0.0; for (int i 0; i n; i) for (int j 0; j n; j) s creal(A[IDX(i, j, n)] * conj(A[IDX(i, j, n)])); return sqrt(s); }mat_mul要求C、A、B三块内存互不重叠否则循环中会读到已经被覆盖的旧值QR迭代里我会用独立矩阵暂存来规避。mat_fnorm返回矩阵的F范数它在收敛判断中用来衡量下三角残余作用比名字看起来重要得多。这三个函数建议单独放到cmat.c里导出头文件后文所有算法只需要包含接口。3. 特征值算法选型为什么这个C语言项目绕不开QR迭代3.1 特征多项式不是正路数值稳定性的坏消息一看到特征值就想到det(A - λI) 0这是线性代数课堂路线但不是数值计算路线。Abel-Ruffini定理决定了五次以上多项式没有通用根式解这是理论层面的死路即使对低阶矩阵强行展开特征多项式得到的系数也是矩阵元素的高次组合微小浮点误差会被放大求根结果常出现完全虚假的复数共轭对。如果final源码走特征多项式路线n4以上开始翻车n8以上基本不可用。原因很直接展开行列式本身就是大量乘加运算浮点误差在过程中累积多项式求根又依赖系数系数误差被放大后根的位置完全失控。QR迭代不构造特征多项式而是通过正交相似变换逐步把矩阵化成上三角对角线直接就是特征值这条路数值上稳定得多也是LAPACK、GSL等库的实际做法。3.2 厄米矩阵可走Jacobi一般复矩阵走QR迭代写代码前先判断矩阵类型。若A^H A即共轭转置等于自身称为厄米矩阵其特征值全为实数自研代码可以用复数Jacobi旋转每次消一个非对角元实现简单。但复数矩阵特征值的一般情况没有厄米约束Jacobi对非厄米矩阵不收敛必须上QR迭代。LAPACK的zgeev就是对一般复矩阵做原位QR迭代并叠加了Hessenberg化、移位和平衡化。先看矩阵有没有结构再看是否需要自己造轮子矩阵类型推荐方法特征值类型实现成本厄米矩阵复数Jacobi旋转实数低一般复矩阵移位QR迭代复数中高任意矩阵工程快速交付LAPACKE zgeev复数极低为什么厄米矩阵特征值一定是实数对特征对Ax λx两边左乘x^H得到x^H A x λ x^H x因为A A^Hx^H A x的共轭等于自身所以λ等于自己的共轭虚部只能为0。这个结论能用来做验证如果输入是厄米矩阵输出却出现明显虚部说明QR迭代代码有bug。3.3 移位QR迭代的收敛逻辑与停止条件QR迭代的基本流程是三步对当前矩阵A_k做QR分解得到酉矩阵Q_k和上三角R_k计算A_{k1} R_k Q_k重复直到下三角元素小到可忽略。因为A_{k1} Q_k^H A_k Q_k每轮都是相似变换特征值集合同原来完全一样。当k足够大时A_k的下三角元素趋向0对角线趋向特征值。收敛速度取决于相邻特征值的模之比模相近时很慢所以实用代码都要加移位每轮从A_k右下角减一个数μ做QR分解得到R Q后再加回μ。复矩阵的μ也可以是复数后面调参章节再展开。停止条件用一个绝对阈值tol。计算下三角元素的模方和开方后小于tol就停止double off 0.0; for (int i 0; i n; i) for (int j 0; j i; j) off creal(A[IDX(i, j, n)] * conj(A[IDX(i, j, n)])); off sqrt(off);阈值取太大会得到错误特征值取太小小则迭代次数爆炸。我一般从1e-10起步高精度需求下放宽到1e-12同时配合最大迭代次数双重保护。4. 源码实现跑通复数矩阵特征值的最终代码4.1 工程最快方案用LAPACKE封装zgeev如果final项目允许链接外部库直接调用LAPACKE_zgeev最省事。这是LAPACK对一般复矩阵zgeev的C接口内部完成了平衡化、Hessenberg化、移位QR以及特征向量的可选计算数值稳定性远好于课程作业级自研代码。调用代码很短#include stdio.h #include complex.h #include lapacke.h void eig_lapack(int n, double complex *A, double complex *w) { int info LAPACKE_zgeev(LAPACK_ROW_MAJOR, N, N, n, A, n, w, NULL, 1, NULL, 1); if (info ! 0) fprintf(stderr, zgeev failed, info%d\n, info); }LAPACK_ROW_MAJOR表示调用方按行主序传入矩阵和前面的IDX布局一致两个N表示不计算左右特征向量A会被内部工作空间覆盖需要保留原矩阵就先拷贝一份w是长度为n的复数数组接收特征值。编译时链接-llapacke -llapack -lblas -lm系统需要装liblapacke-dev。缺点是部分环境没有预装LAPACK交叉编译到嵌入式设备时更麻烦这时自研QR迭代就成了实际可行的备选。4.2 自主实现复数Householder QR分解与QR迭代主循环自己写QR迭代复杂度集中在QR分解。这里用Householder变换实现复数版QR分解对第k列取出子列x构造反射向量v使Hx除第一个元素外都变成0。复数版与实数版的差别在两点一是反射向量要考虑x[0]的相位alpha -(x[0] / |x[0]|) * ||x||这样算出的v模最大数值最稳二是更新矩阵时用共轭转置v^H而不是普通转置。static void qr_decomp(cx *A, cx *Q, int n) { for (int i 0; i n; i) for (int j 0; j n; j) Q[IDX(i, j, n)] (i j) ? 1.0 : 0.0; for (int k 0; k n - 1; k) { int m n - k; cx *x (cx *)malloc(sizeof(cx) * m); for (int i 0; i m; i) x[i] A[IDX(k i, k, n)]; double nrm 0.0; for (int i 0; i m; i) nrm creal(x[i] * conj(x[i])); nrm sqrt(nrm); if (nrm 1e-300) { free(x); continue; } cx alpha (cabs(x[0]) 1e-300) ? -nrm : -(x[0] / cabs(x[0])) * nrm; cx *v (cx *)malloc(sizeof(cx) * m); v[0] x[0] - alpha; for (int i 1; i m; i) v[i] x[i]; double vh_v 0.0; for (int i 0; i m; i) vh_v creal(v[i] * conj(v[i])); if (vh_v 1e-300) { free(x); free(v); continue; } double scale 2.0 / vh_v; // A H A A - scale * v * (v^H A) for (int j k; j n; j) { cx dot 0.0; for (int l 0; l m; l) dot conj(v[l]) * A[IDX(k l, j, n)]; dot * scale; for (int i 0; i m; i) A[IDX(k i, j, n)] - v[i] * dot; } // Q Q H Q - scale * (Q v) * v^H cx *u (cx *)malloc(sizeof(cx) * n); for (int i 0; i n; i) { cx s 0.0; for (int l 0; l m; l) s Q[IDX(i, k l, n)] * v[l]; u[i] s * scale; } for (int i 0; i n; i) for (int j 0; j m; j) Q[IDX(i, k j, n)] - u[i] * conj(v[j]); free(x); free(v); free(u); } } void eig_qr(cx *A, cx *w, int n, int max_iter, double tol) { cx *Q mat_new(n); cx *R mat_new(n); int it; for (it 0; it max_iter; it) { qr_decomp(A, Q, n); // A_new R * Q用独立矩阵 R 暂存再拷回 for (int i 0; i n; i) for (int j 0; j n; j) { cx s 0.0; for (int k 0; k n; k) s A[IDX(i, k, n)] * Q[IDX(k, j, n)]; R[IDX(i, j, n)] s; } for (int i 0; i n; i) for (int j 0; j n; j) A[IDX(i, j, n)] R[IDX(i, j, n)]; double off 0.0; for (int i 0; i n; i) for (int j 0; j i; j) off creal(A[IDX(i, j, n)] * conj(A[IDX(i, j, n)])); if (sqrt(off) tol) break; } for (int i 0; i n; i) w[i] A[IDX(i, i, n)]; mat_free(Q); mat_free(R); }qr_decomp把A原地覆盖成R同时生成Q。更新A与Q的顺序可以互换因为二者都只依赖旧A和旧Q不存在相互覆盖。eig_qr中R * Q用独立矩阵暂存是必须的A此时就是R直接在A上做A A * Q会覆盖尚未读取的元素。收敛后A接近上三角对角线就是特征值。这个实现是无移位版本矩阵条件数好时能收敛遇到收敛慢的矩阵把max_iter调大或者按下一章的加移位方式改造。4.3 编译命令与链接参数把上述函数和一个main放在eigen.c里编译命令是gcc -stdc99 -O2 -Wall -o eigen eigen.c -lm-lm不能省复数运算中的conj、cabs、sqrt都在libm里。-O2对数值代码影响很大不开优化时浮点运算顺序可能不同特征值结果会在最后几位抖动提交final项目前务必用优化编译重测一遍。如果加了LAPACKE分支编译命令变成gcc -stdc99 -O2 -o eigen eigen.c -llapacke -llapack -lblas -lm链接顺序要按依赖关系排-llapacke在最前因为它依赖-llapack和-lblas。报错出现undefined reference to zgeev_多半是lapack没装或链接顺序写反了。5. 验证与调参final版复数矩阵特征值C语言源码的收尾检查5.1 用已知复特征值的矩阵验证输出验证是final项目最该做却最容易被省掉的一步。用一个已知特征值的实矩阵[[2, 1], [-1, 2]]其特征多项式是λ^2 - 4λ 5特征值应为2 i与2 - i。测试代码cx A[4] {2.0, 1.0, -1.0, 2.0}; cx w[2]; eig_qr(A, w, 2, 1000, 1e-10); for (int i 0; i 2; i) printf(%.6f %.6f i\n, creal(w[i]), cimag(w[i]));理想输出是2.000000 1.000000i和2.000000 -1.000000i。如果虚部符号相反只是特征值顺序不同这是正常的特征值本身没有固定顺序比较前先排序。如果虚部完全消失大概率是QR分解更新A时把conj(v[l])写成了v[l]这是复数Householder变换最经典的错误。5.2 收敛公差、最大迭代次数与移位策略参数建议初值后果与调试方向tol1e-10偏大时特征值误差明显偏小则迭代次数多余max_iter500 * n达到上限未break先检查是否忘记加移位shiftA[n-1][n-1]加速收敛的关键复矩阵用复μ矩阵规模n 2..50n大于50应先用Householder化为Hessenberg形无移位QR对模相近的特征值收敛很慢。常见做法是每轮取μ A[n-1][n-1]作为移位先算A A - μI对A做QR更新A_new R * Q μI。右下角元素在迭代中逐步接近某个特征值这个移位能显著加快收敛复矩阵就取复μ。收敛后对角线就是特征值不需要再减μ。5.3 三个容易被忽略的边界问题特征值顺序不稳定。同一矩阵在不同n或不同轮数下输出顺序可能变化做对比测试时按实部排序实部相近再按虚部排序。矩阵尺度差异过大。元素跨10个数量级时直接QR可能因上溢得到NaN先做平衡化行列同时缩放LAPACK的zgebal专门做这件事。子列范数恰好为0。qr_decomp遇到nrm 1e-300直接跳过该列这对奇异或稀疏结构是合法的但主循环收敛判据仍会统计该列下三角残余迭代次数可能变多不应误判为死循环。把平衡化加进final项目通常能让特征值精度回到1e-10量级。本文还有配套的精品资源点击获取
