diff --git a/CustomMath_lib/IterationMethod/PowerIteration.cpp b/CustomMath_lib/IterationMethod/PowerIteration.cpp index 4e05b39..b51bbb4 100644 --- a/CustomMath_lib/IterationMethod/PowerIteration.cpp +++ b/CustomMath_lib/IterationMethod/PowerIteration.cpp @@ -4,17 +4,42 @@ #include "PowerIteration.h" -lld PowerIteration(const PowerIterationInput &input) { +#define ITERATION_METHOD_MAX_ITERATION 100000 + +#define ENABLE_POWER_ITERATION_METHOD_TIMING + +#ifdef ENABLE_POWER_ITERATION_METHOD_TIMING +#define ITERATION_METHOD_TIMING_START auto start = std::chrono::high_resolution_clock::now(); +#define ITERATION_METHOD_TIMING_END auto end = std::chrono::high_resolution_clock::now(); \ + auto duration = std::chrono::duration_cast(end - start); +#define ITERATION_METHOD_RETURN_DURATION duration +#else +#define ITERATION_METHOD_TIMING_START +#define ITERATION_METHOD_TIMING_END +#define ITERATION_METHOD_RETURN_DURATION std::chrono::microseconds(0) +#endif + +PowerIterationOutput PowerIteration(const PowerIterationInput &input) { auto A = input.A; auto x = input.x_default; + + ITERATION_METHOD_TIMING_START + for (int i = 0; i < input.iteration_times; i++) { x = x / VectorNorm_Infinity(x); x = A * x; } - return VectorNorm_Infinity(x); + + ITERATION_METHOD_TIMING_END + + return { // + VectorNorm_Infinity(x), // + input.iteration_times, // + ITERATION_METHOD_RETURN_DURATION // + }; } -lld MaxRootForPolynomial(const Vector &coefficients) { +PowerIterationOutput MaxRootForPolynomial(const Vector &coefficients) { auto n = coefficients.size; auto A = Matrix(n, n); @@ -29,6 +54,5 @@ lld MaxRootForPolynomial(const Vector &coefficients) { auto x = Vector(n, 0); x.array[0] = 1; - auto result = PowerIteration({A, x, 100000}); - return result; + return PowerIteration({A, x, ITERATION_METHOD_MAX_ITERATION}); } \ No newline at end of file diff --git a/CustomMath_lib/IterationMethod/PowerIteration.h b/CustomMath_lib/IterationMethod/PowerIteration.h index 62979bf..ae35193 100644 --- a/CustomMath_lib/IterationMethod/PowerIteration.h +++ b/CustomMath_lib/IterationMethod/PowerIteration.h @@ -7,6 +7,7 @@ #include "Matrix/Matrix.h" #include "Vector/Vector.h" +#include "chrono" #include "InfinityNorm/InfinityNorm.h" typedef struct { @@ -15,7 +16,13 @@ typedef struct { int iteration_times; } PowerIterationInput; -lld PowerIteration(const PowerIterationInput &input); +typedef struct { + lld x; + int iteration_times; + std::chrono::microseconds time_cost; +} PowerIterationOutput; + +PowerIterationOutput PowerIteration(const PowerIterationInput &input); /* * @brief: Find the max root for a polynomial @@ -23,6 +30,6 @@ lld PowerIteration(const PowerIterationInput &input); * @param: coefficients: The coefficients of the polynomial, x^n + a_{1}x^{n-1} + ... + a_{n-1}x + a_{n} is assumed * @return: The max root for the polynomial */ -lld MaxRootForPolynomial(const Vector &coefficients); +PowerIterationOutput MaxRootForPolynomial(const Vector &coefficients); #endif //NUMERICAL_ALGEBRA_POWERITERATION_H diff --git a/CustomMath_lib/QRMethod/DoubleStepQRIteration.cpp b/CustomMath_lib/QRMethod/DoubleStepQRIteration.cpp index a9d2a71..794ea23 100644 --- a/CustomMath_lib/QRMethod/DoubleStepQRIteration.cpp +++ b/CustomMath_lib/QRMethod/DoubleStepQRIteration.cpp @@ -12,8 +12,8 @@ Matrix DoubleStepQRIteration(const Matrix &matrix) { auto n = matrix.rows; auto m = n - 1; - auto s = H.matrix[m-1][m-1] + H.matrix[n-1][n-1]; - auto t = H.matrix[m-1][m-1] * H.matrix[n-1][n-1] - H.matrix[m-1][n-1] * H.matrix[n-1][m-1]; + auto s = H.matrix[m - 1][m - 1] + H.matrix[n - 1][n - 1]; + auto t = H.matrix[m - 1][m - 1] * H.matrix[n - 1][n - 1] - H.matrix[m - 1][n - 1] * H.matrix[n - 1][m - 1]; auto x = H.matrix[0][0] * H.matrix[0][0] + H.matrix[0][1] * H.matrix[1][0] - s * H.matrix[0][0] + t; auto y = H.matrix[1][0] * (H.matrix[0][0] + H.matrix[1][1] - s); auto z = H.matrix[1][0] * H.matrix[2][1]; @@ -21,38 +21,38 @@ Matrix DoubleStepQRIteration(const Matrix &matrix) { Vector b; lld beta; - for(ull k = 0; k < n - 2; k++) { - auto tmpVector = Vector(std::vector {x, y, z}); + for (ull k = 0; k < n - 2; k++) { + auto tmpVector = Vector(std::vector{x, y, z}); HouseHolderMethod(tmpVector, b, beta); - auto w = product(b,b)* beta; + auto w = product(b, b) * beta; auto q = MAX(1, k); - auto H_sub = H.sub_matrix(k,k+3,q-1,n); + auto H_sub = H.sub_matrix(k, k + 3, q - 1, n); H_sub = H_sub - w * H_sub; - H.set(k,k+3,q-1,n, H_sub); + H.set(k, k + 3, q - 1, n, H_sub); - auto r = MIN(k + 3,n); - H_sub = H.sub_matrix(0,r,k,k+3); + auto r = MIN(k + 3, n); + H_sub = H.sub_matrix(0, r, k, k + 3); H_sub = H_sub - H_sub * w; - H.set(0,r,k,k+3, H_sub); + H.set(0, r, k, k + 3, H_sub); - x = H.matrix[k+1][k]; - y = H.matrix[k+2][k]; - if(k < n - 3) { - z = H.matrix[k+3][k]; + x = H.matrix[k + 1][k]; + y = H.matrix[k + 2][k]; + if (k < n - 3) { + z = H.matrix[k + 3][k]; } } - auto tmpVector = Vector(std::vector {x, y}); + auto tmpVector = Vector(std::vector{x, y}); HouseHolderMethod(tmpVector, b, beta); - auto w = product(b,b)* beta; + auto w = product(b, b) * beta; - auto H_sub = H.sub_matrix(n-2,n,n-3,n); + auto H_sub = H.sub_matrix(n - 2, n, n - 3, n); H_sub = H_sub - w * H_sub; - H.set(n-2,n,n-3,n, H_sub); + H.set(n - 2, n, n - 3, n, H_sub); - H_sub= H.sub_matrix(0,n,n-2,n); + H_sub = H.sub_matrix(0, n, n - 2, n); H_sub = H_sub - H_sub * w; - H.set(0,n,n-2,n, H_sub); + H.set(0, n, n - 2, n, H_sub); return H.clean(); } @@ -65,8 +65,8 @@ Matrix DoubleStepQRIteration(Matrix &matrix, Matrix &P) { auto n = matrix.rows; auto m = n - 1; - auto s = H.matrix[m-1][m-1] + H.matrix[n-1][n-1]; - auto t = H.matrix[m-1][m-1] * H.matrix[n-1][n-1] - H.matrix[m-1][n-1] * H.matrix[n-1][m-1]; + auto s = H.matrix[m - 1][m - 1] + H.matrix[n - 1][n - 1]; + auto t = H.matrix[m - 1][m - 1] * H.matrix[n - 1][n - 1] - H.matrix[m - 1][n - 1] * H.matrix[n - 1][m - 1]; auto x = H.matrix[0][0] * H.matrix[0][0] + H.matrix[0][1] * H.matrix[1][0] - s * H.matrix[0][0] + t; auto y = H.matrix[1][0] * (H.matrix[0][0] + H.matrix[1][1] - s); auto z = H.matrix[1][0] * H.matrix[2][1]; @@ -74,42 +74,42 @@ Matrix DoubleStepQRIteration(Matrix &matrix, Matrix &P) { Vector b; lld beta; - for(ull k = 0; k < n - 2; k++) { - auto tmpVector = Vector(std::vector {x, y, z}); + for (ull k = 0; k < n - 2; k++) { + auto tmpVector = Vector(std::vector{x, y, z}); HouseHolderMethod(tmpVector, b, beta); - auto w = product(b,b)* beta; + auto w = product(b, b) * beta; auto q = MAX(1, k); - auto H_sub = H.sub_matrix(k,k+3,q-1,n); + auto H_sub = H.sub_matrix(k, k + 3, q - 1, n); H_sub = H_sub - w * H_sub; - H.set(k,k+3,q-1,n, H_sub); + H.set(k, k + 3, q - 1, n, H_sub); - auto r = MIN(k + 3,n); - H_sub = H.sub_matrix(0,r,k,k+3); + auto r = MIN(k + 3, n); + H_sub = H.sub_matrix(0, r, k, k + 3); H_sub = H_sub - H_sub * w; - H.set(0,r,k,k+3, H_sub); + H.set(0, r, k, k + 3, H_sub); - auto P_sub = P.sub_matrix(0,n,k,k+3); + auto P_sub = P.sub_matrix(0, n, k, k + 3); P_sub = P_sub - P_sub * w; - P.set(0,n,k,k+3, P_sub); + P.set(0, n, k, k + 3, P_sub); - x = H.matrix[k+1][k]; - y = H.matrix[k+2][k]; - if(k < n - 3) { - z = H.matrix[k+3][k]; + x = H.matrix[k + 1][k]; + y = H.matrix[k + 2][k]; + if (k < n - 3) { + z = H.matrix[k + 3][k]; } } - auto tmpVector = Vector(std::vector {x, y}); + auto tmpVector = Vector(std::vector{x, y}); HouseHolderMethod(tmpVector, b, beta); - auto w = product(b,b)* beta; + auto w = product(b, b) * beta; - auto H_sub = H.sub_matrix(n-2,n,n-3,n); + auto H_sub = H.sub_matrix(n - 2, n, n - 3, n); H_sub = H_sub - w * H_sub; - H.set(n-2,n,n-3,n, H_sub); + H.set(n - 2, n, n - 3, n, H_sub); - H_sub= H.sub_matrix(0,n,n-2,n); + H_sub = H.sub_matrix(0, n, n - 2, n); H_sub = H_sub - H_sub * w; - H.set(0,n,n-2,n, H_sub); + H.set(0, n, n - 2, n, H_sub); return H.clean(); } \ No newline at end of file diff --git a/CustomMath_lib/QRMethod/QRMethod.cpp b/CustomMath_lib/QRMethod/QRMethod.cpp index 841f999..70ca6de 100644 --- a/CustomMath_lib/QRMethod/QRMethod.cpp +++ b/CustomMath_lib/QRMethod/QRMethod.cpp @@ -4,15 +4,31 @@ #include "QRMethod.h" -Matrix QRMethod(const Matrix &matrix) { +#define MAX_ITERATION 100000 +#define ENABLE_TIMING + +#ifdef ENABLE_TIMING +#define ITERATION_METHOD_TIMING_START auto start = std::chrono::high_resolution_clock::now(); +#define ITERATION_METHOD_TIMING_END auto end = std::chrono::high_resolution_clock::now(); \ + auto duration = std::chrono::duration_cast(end - start); +#define ITERATION_METHOD_RETURN_DURATION duration +#else +#define ITERATION_METHOD_TIMING_START +#define ITERATION_METHOD_TIMING_END +#define ITERATION_METHOD_RETURN_DURATION std::chrono::microseconds(0) +#endif + +QRMethodOutput QRMethod(const Matrix &matrix) { Matrix H = matrix; Matrix Q; auto P = Matrix::identity(H.rows); auto n = H.rows; + ITERATION_METHOD_TIMING_START + HessenbergMethod_Inplace(H, Q); - do { + for (int count = 0; count < MAX_ITERATION; count++) { // set all abs(h_{i,i-1}) < 1e-10 to 0 for (ull i = 0; i < H.rows - 1; i++) { if (std::abs(H.matrix[i + 1][i]) < 1e-6) { @@ -54,7 +70,13 @@ Matrix QRMethod(const Matrix &matrix) { l = n - m - k; if (m == n || n - m - l <= 2) { - break; + ITERATION_METHOD_TIMING_END + + return { // + H.clean(), // + count, // + ITERATION_METHOD_RETURN_DURATION // + }; } auto H22 = H.sub_matrix(l, n - m, l, n - m); @@ -72,12 +94,10 @@ Matrix QRMethod(const Matrix &matrix) { H.set(0, l, l, n - m, H12); H.set(l, n - m, n - m, n, H23); - } while (true); - - return H.clean(); + } } -Vector AllRootsForPolynomial(const Vector &coefficients) { +QRMethodOutput AllRootsForPolynomial(const Vector &coefficients) { auto n = coefficients.size; auto A = Matrix(n, n); @@ -93,18 +113,22 @@ Vector AllRootsForPolynomial(const Vector &coefficients) { auto result = Vector(n); for (ull i = 0; i < n; i++) { - result.array[i] = r.matrix[i][i]; + result.array[i] = r.result.matrix[i][i]; } - return result; + return { // + result, // + r.iteration_times, // + r.time_cost // + }; } std::vector AllEigenValues(const Matrix &R) { auto n = R.rows; auto result = std::vector(n); - for (ull i = 0; i +struct QRMethodOutput { + T result; + int iteration_times; + std::chrono::microseconds time_cost; +}; + + +QRMethodOutput QRMethod(const Matrix &matrix); + +QRMethodOutput AllRootsForPolynomial(const Vector &coefficients); + std::vector AllEigenValues(const Matrix &R); + void print_llc(const std::vector vec); #endif //NUMERICAL_ALGEBRA_QRMETHOD_H diff --git a/homeworks/homework_6.cpp b/homeworks/homework_6.cpp index dfcbe16..5a9c69a 100644 --- a/homeworks/homework_6.cpp +++ b/homeworks/homework_6.cpp @@ -9,7 +9,9 @@ void par_1_each(Vector &a) { std::cout << "Polynomial is "; a.print(); - std::cout << "Max root for polynomial is " << std::setprecision(10) << max_root << std::endl; + std::cout << "Max root for polynomial is " << std::setprecision(10) << max_root.x << std::endl; + std::cout<< "Iteration times is " << max_root.iteration_times << std::endl; + std::cout << "Time cost is " << max_root.time_cost.count() << " microseconds" << std::endl; } @@ -53,7 +55,9 @@ void par_2_2() { auto result = AllRootsForPolynomial(coefficients); std::cout << "Roots for polynomial are:" << std::endl; - result.print(); + result.result.print(); + std::cout << "Iteration times is " << result.iteration_times << std::endl; + std::cout << "Time cost is " << result.time_cost.count() << " microseconds" << std::endl; } void par_2_3() { @@ -81,9 +85,11 @@ void par_2_3() { {6.1, 4.9, 3.5, 6.2}}); auto h = QRMethod(m); - auto k = AllEigenValues(h); + auto k = AllEigenValues(h.result); std::cout << "Eigen values are:" << std::endl; print_llc(k); + std::cout << "Iteration times is " << h.iteration_times << std::endl; + std::cout << "Time cost is " << h.time_cost.count() << " microseconds" << std::endl; std::cout << std::endl; } diff --git a/reports/data/report_6_output.txt b/reports/data/report_6_output.txt new file mode 100644 index 0000000..b29aeb6 --- /dev/null +++ b/reports/data/report_6_output.txt @@ -0,0 +1,47 @@ +------ Q 6.1 ------ +Polynomial is [1,-5,3] +Max root for polynomial is 3 +Iteration times is 100000 +Time cost is 7580 microseconds +Polynomial is [0,-3,-1] +Max root for polynomial is 1.879385242 +Iteration times is 100000 +Time cost is 16660 microseconds +Polynomial is [101,208.01,10891.01,9802.08,79108.9,-99902,790,-1000] +Max root for polynomial is 100 +Iteration times is 100000 +Time cost is 39453 microseconds + +------ Q 6.2(2) ------ +Roots for polynomial are: +[0.9377142935,0.9968532977,0.8160073098,1.479041989,1.339234306,0.1165772655,0.4187154748,0.1967653274,0.4348249074,0.9172620712,0.8083910741,0.3900648089,0.4468288339,0.6955448095,0.6582632323,0.9368506979,0.2291885706,0.2278977162,0.0913795626,0.07848798633,-0.1041481049,-0.09830129001,0.1439366857,0.1649965994,-0.2819267664,-0.2842291233,-0.4558137458,-0.4457402621,-0.6045910143,-0.5946672966,-0.7371247322,-0.7165505332,-0.8389326301,-0.821359886,-0.9219128266,-0.8924451194,-0.940105485,-0.9696987081,-0.9365104319,-0.9983448119,-0.9516475702] +Iteration times is 4644 +Time cost is 1786975 microseconds + +------ Q 6.2(3) ------ +x = 0.9 +Eigen values are: +17.41702451 + 0i +2.955106057 + 0.8459755059i +2.955106057 + -0.8459755059i +6.672764737 + 0i +Iteration times is 10 +Time cost is 71 microseconds + +x = 1 +Eigen values are: +17.45296588 + 0i +6.520993181 + 0i +3.013020446 + 1.163763663i +3.013020446 + -1.163763663i +Iteration times is 8 +Time cost is 53 microseconds + +x = 1.1 +Eigen values are: +17.48872623 + 0i +6.506068455 + 0i +3.002575229 + 1.394322746i +3.002575229 + -1.394322746i +Iteration times is 7 +Time cost is 51 microseconds \ No newline at end of file diff --git a/reports/report_6.pdf b/reports/report_6.pdf new file mode 100644 index 0000000..e8c6b33 Binary files /dev/null and b/reports/report_6.pdf differ diff --git a/reports/report_6.tex b/reports/report_6.tex new file mode 100644 index 0000000..9364baa --- /dev/null +++ b/reports/report_6.tex @@ -0,0 +1,77 @@ +\documentclass{article} +\usepackage[UTF8]{ctex} +\usepackage{float,indentfirst,verbatim,fancyhdr,graphicx,listings,longtable,amsmath, amsfonts,amssymb} + +\textheight 23.5cm \textwidth 15.8cm +\topmargin -1.5cm \oddsidemargin 0.3cm \evensidemargin -0.3cm + +\title{数值代数实验报告 6} +\author{马天开} + +\begin{document} +\maketitle + +\section{问题描述} + +\subsection{求多项式方程的模最大根} + +\subsubsection{}用 C++ 编制利用幂法求多项式方程 +\[ + f(x) = x^n + \alpha_{n-1}x^{n-1} +\cdots+ \alpha_1 x + \alpha_0 = 0 +\] 的模最大根的通用子程序。 + +\subsubsection{}利用你所编制的子程序求下列各高次方程的模最大根 + +\begin{enumerate} + \item $x^3 + x^2 - 5x + 3 = 0;$ + \item $x^3 - 3x - 1 = 0;$ + \item $x^8 + 101x^7 + 208.01x^6 + 10891.01x^5 + 9802.08x^4 + 79108.9x^3-99902x^2 + 790x-1000 = 0.$ +\end{enumerate} + +\textbf{要求输出迭代次数,用时和最大根的值(注意正负)} + + +\subsection{求实矩阵的全部特征值} + +\subsubsection{}用 C++ 编制利用隐式 QR 算法 (课本算法 6.4.3) 求一个实矩阵的全部特征值的通用子程序。 + +\subsubsection{利用你所编制的子程序计算方程} +$$x^{41} + x^3 + 1 = 0 $$ + +的全部根。 + +\subsubsection{} +设 +\[ + A=\begin{bmatrix} + \;9.1 & 3.0 & 2.6 & 4.0 \;\\ + \;4.2 & 5.3 & 4.7 & 1.6 \;\\ + \;3.2 & 1.7 & 9.4 & x \;\\ + \;6.1 & 4.9 & 3.5 & 6.2\; + \end{bmatrix} +\] + +求当 $x = 0.9, 1.0, 1.1$ 时 A 的全部特征值,并观察并在报告中叙述分析特征值实部、虚部和模长的变化情况。 + +\textbf{要求输出迭代次数、用时和所有特征值。} + +\section{算法说明} + +必须实现的算法有: + +\begin{enumerate} + \item 幂法求模最大根 $\Rightarrow$ \verb|IterationMethod/PowerIteration| + \item Hessenberg 分解 $\Rightarrow$ \verb|HouseHolderMethod/HessenbergMethod| + \item 双重步位移的 QR 迭代 $\Rightarrow$ \verb|QRMethod/DoubleStepQRIteration| + \item 隐式 QR 算法 $\Rightarrow$ \verb|QRMethod/QRMethod| +\end{enumerate} + +有一些英文命名并不是很规范,有时间会调整的。 + +\section{运行结果} + +\verbatiminput{data/report_6_output.txt} + +\newpage + +\end{document} \ No newline at end of file