对 Jacobi 程序的循环优化

最新推荐文章于 2024-09-16 17:46:36 发布
最新推荐文章于 2024-09-16 17:46:36 发布
阅读量75 收藏
点赞数
文章标签： c++ 算法 开发语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Fallen250/article/details/131681872
版权

Jacobi 程序是矩阵计算中常见的算法，程序中存在大量的循环语句，通过使用循环优化技术，可以降低程序的运行时间，以下进行详细介绍。

原始程序如下所示：

#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <stdlib.h>
#include <array>
#include <cmath>
#include <stdint.h>
#include <vector>
#include <algorithm>
#include <iterator>
#include <iostream>


#define N 202
#define K 202

struct RectanglePoints {
	double x_0, x_l;
	double y_0, y_m;
};

namespace InitialFunctions {
	inline double leftSideRectangle(double y) { return 15 * sqrt(y); }
	inline double rightSideRectangle(double y) { return 30 * y * (1 - y); }
	inline double upSideRectangle(double x) { return 15 * (1 - x); }
	inline double downSideRectangle(double x) { return 0; }

}

using namespace InitialFunctions;

class Grid {
private:

	RectanglePoints rectangle;

	double length_x;
	double length_y;
	double h_xyStep;
	
public:
	Grid(); //default constructor
	Grid(double l, double m, double h, std::array<double, N>& x, std::array<double, N>& y);

    void setPoints();
	void createGrid(std::array<double, N>& x, std::array<double, K>& y);
	void getGrid(std::array<double, N>& x, std::array<double, K>& y);
};

Grid::Grid()
{
	length_x = 1;
	length_y = 1;
	h_xyStep = 0.2;

	setPoints();
}

Grid::Grid(double l, double m, double h, std::array<double, N>& x, std::array<double, K>& y) {

	length_x = l;
	length_y = m;
	h_xyStep = h;

	setPoints();
	createGrid(x, y);
}

void Grid::setPoints()
{

	rectangle.x_0 = 0;
	rectangle.y_0 = 0;
	rectangle.x_l = rectangle.x_0 + length_x;
	rectangle.y_m = rectangle.y_0 + length_y;
}

void Grid::createGrid(std::array<double, N>& x, std::array<double, K>& y)
{
	x[0] = rectangle.x_0;
	y[0] = rectangle.y_0;
	
	for (uint8_t i = 0; i < N; ++i)
	{
		x[i] = i * h_xyStep;
	}

	for (uint8_t j = 0; j < K; ++j)
	{
		y[j] = j * h_xyStep;
	}
}

void Grid::getGrid(std::array<double, N>& x, std::array<double, K>& y)
{
	for (uint8_t i = 0; i < N; i++)
		printf("%lf ", x[i]);
	printf("\n\n");
	for (uint8_t j = 0; j < K; j++)
		printf("%lf ", y[j]);
	printf("\n\n");

}


class LaplasEquationSolver {
public:
	double w;
	LaplasEquationSolver(double(*u1)[N], std::array<double, N>& x, std::array<double, K>& y, double w);
};

LaplasEquationSolver::LaplasEquationSolver(double(*u1)[N], std::array<double, N>& x, std::array<double, K>& y, double w)
{
	double u2[N][K] = { 0 }, u_new[N][K] = { 0 };
	std::vector<double> errors;
	double max_error;
	double eps = 0.01;

	for (uint8_t i = 0; i < N; ++i)
		for (uint8_t j = 0; j < K; ++j) {
			u1[i][j] = 0;
			u2[i][j] = 0;
			u_new[i][j] = 0;
		}


	for (uint8_t i = 0; i < N; ++i)
	{
		u1[i][0] = downSideRectangle(x[i]);
		u1[i][K - 1] = upSideRectangle(x[i]);
	}

	for (uint8_t j = 0; j < K; ++j)
	{
		u1[0][j] = leftSideRectangle(y[j]);
		u1[N - 1][j] = rightSideRectangle(y[j]);
	}

	for (uint8_t i = 0; i < N; ++i)
	{
		u2[i][0] = downSideRectangle(x[i]);
		u2[i][K - 1] = upSideRectangle(x[i]);
	}

	for (uint8_t j = 0; j < K; ++j)
	{
		u2[0][j] = leftSideRectangle(y[j]);
		u2[N - 1][j] = rightSideRectangle(y[j]);
	}

	int k = 0;
	clock_t start_s=clock();
	do {
		for (uint8_t j = 1; j < K - 1; ++j) {
			for (uint8_t i = 1; i < N - 1; ++i) {
				u1[i][j] = (1.0 / 4.0) * (u1[i - 1][j] + u1[i + 1][j] + u1[i][j - 1] + u1[i][j + 1]);
			}
		}

		for (uint8_t j = 1; j < K - 1; ++j) {
			for (uint8_t i = 1; i < N - 1; ++i) {
				u2[i][j] = (1.0 / 4.0) * (u1[i - 1][j] + u1[i + 1][j] + u1[i][j - 1] + u1[i][j + 1]);
			}
		}

		for (uint8_t j = 0; j < K; ++j) {
			for (uint8_t i = 0; i < N; ++i) {
				u_new[i][j] = u1[i][j] + w * (u2[i][j] - u1[i][j]);
			}
		}

		for (uint8_t i = 0; i < N; ++i) {
			for (uint8_t j = 0; j < K; ++j) {
				errors.push_back(fabs(u_new[i][j] - u1[i][j]));;
			}
		}
		max_error = *std::max_element(errors.begin(), errors.end());
		//if(k%100==0)std::cout<<k<<" "<<max_error<<"\n";
		errors.clear();
		for (uint8_t i = 0; i < N; ++i) {
			for (uint8_t j = 0; j < K; ++j) {
				u1[i][j] = u_new[i][j];
			}
		}
		k++;
	} while (max_error > eps);
	clock_t end_s=clock();

	std::cout << "Input = "<<N<<" "<<K <<std::endl;
	std::cout << "w= " << w << ", iteration= " << k << ", max error = " << max_error << std::endl;
	std::cout << "cost time = "<< end_s - start_s << " clocks"<< std::endl;
}




int main(void) 
{
	std::array<double, N> x = { 0 };
	std::array<double, K> y = { 0 };
	double u[N][K] = { 0 };

	/*calling 2nd constructor*/
	Grid grid(1, 1, 0.2, x, y); 

	/*print grid to console*/
	//grid.getGrid(x, y);

	double w = 1.9;
	LaplasEquationSolver mySolver(u, x, y, w);



	return 0;
}

 1、对原始程序进行编译，并查看运行时间

 2、使用循环交换进行优化

do {
		
		for (uint8_t i = 1; i < N - 1; ++i) {
			for (uint8_t j = 1; j < K - 1; ++j) {
				u1[i][j] = (1.0 / 4.0) * (u1[i - 1][j] + u1[i + 1][j] + u1[i][j - 1] + u1[i][j + 1]);
			}
		}

		for (uint8_t i = 1; i < N - 1; ++i) {
			for (uint8_t j = 1; j < K - 1; ++j) {
				u2[i][j] = (1.0 / 4.0) * (u1[i - 1][j] + u1[i + 1][j] + u1[i][j - 1] + u1[i][j + 1]);
			}
		}

		
		for (uint8_t i = 0; i < N; ++i) {
			for (uint8_t j = 0; j < K; ++j) {
				u_new[i][j] = u1[i][j] + w * (u2[i][j] - u1[i][j]);
			}
		}

		
		for (uint8_t i = 0; i < N; ++i) {
			for (uint8_t j = 0; j < K; ++j) {
				errors.push_back(fabs(u_new[i][j] - u1[i][j]));;
			}
		}
		max_error = *std::max_element(errors.begin(), errors.end());
		//if(k%100==0)std::cout<<k<<" "<<max_error<<"\n";
		errors.clear();
		for (uint8_t i = 0; i < N; ++i) {
			for (uint8_t j = 0; j < K; ++j) {
				u1[i][j] = u_new[i][j];
			}
		}
		k++;
	} while (max_error > eps);
	clock_t end_s=clock();

编译后查看运行时间

3、 使用循环不变量外提进行优化

do {
		double coff = 1.0 / 4.0;
		for (uint8_t j = 1; j < K - 1; ++j) {
			for (uint8_t i = 1; i < N - 1; ++i) {
				u1[i][j] = coff * (u1[i - 1][j] + u1[i + 1][j] + u1[i][j - 1] + u1[i][j + 1]);
			}
		}

		for (uint8_t j = 1; j < K - 1; ++j) {
			for (uint8_t i = 1; i < N - 1; ++i) {
				u2[i][j] = coff * (u1[i - 1][j] + u1[i + 1][j] + u1[i][j - 1] + u1[i][j + 1]);
			}
		}

		for (uint8_t j = 0; j < K; ++j) {
			for (uint8_t i = 0; i < N; ++i) {
				u_new[i][j] = u1[i][j] + w * (u2[i][j] - u1[i][j]);
			}
		}

		for (uint8_t i = 0; i < N; ++i) {
			for (uint8_t j = 0; j < K; ++j) {
				errors.push_back(fabs(u_new[i][j] - u1[i][j]));;
			}
		}
		max_error = *std::max_element(errors.begin(), errors.end());
		//if(k%100==0)std::cout<<k<<" "<<max_error<<"\n";
		errors.clear();
		for (uint8_t i = 0; i < N; ++i) {
			for (uint8_t j = 0; j < K; ++j) {
				u1[i][j] = u_new[i][j];
			}
		}
		k++;
	} while (max_error > eps);
	clock_t end_s=clock();

编译后查看运行时间

 4、使用循环展开进行优化

5、使用循环交换、循环不变量外提、循环展开混合技术进行优化

do {
		double coff = 1.0 / 4.0;
		for (uint8_t i = 1; i < N - 1; ++i) {
			for (uint8_t j = 1; j < K - 1; j += 4) {
				u1[i][j] = coff * (u1[i - 1][j] + u1[i + 1][j] + u1[i][j - 1] + u1[i][j + 1]);
				u1[i][j + 1] = coff * (u1[i - 1][j + 1] + u1[i + 1][j + 1] + u1[i][j] + u1[i][j + 2]);
				u1[i][j + 2] = coff * (u1[i - 1][j + 2] + u1[i + 1][j + 2] + u1[i][j + 1] + u1[i][j + 3]);
				u1[i][j + 3] = coff * (u1[i - 1][j + 3] + u1[i + 1][j + 3] + u1[i][j + 2] + u1[i][j + 4]);
			}
		}

		for (uint8_t i = 1; i < N - 1; ++i) {
			for (uint8_t j = 1; j < K - 1; j += 4) {
				u2[i][j] = coff * (u1[i - 1][j] + u1[i + 1][j] + u1[i][j - 1] + u1[i][j + 1]);
				u2[i][j + 1] = coff * (u1[i - 1][j + 1] + u1[i + 1][j + 1] + u1[i][j] + u1[i][j + 2]);
				u2[i][j + 2] = coff * (u1[i - 1][j + 2] + u1[i + 1][j + 2] + u1[i][j + 1] + u1[i][j + 3]);
				u2[i][j + 3] = coff * (u1[i - 1][j + 3] + u1[i + 1][j + 3] + u1[i][j + 2] + u1[i][j + 4]);
			}
		}

		for (uint8_t i = 0; i < N; ++i) {
			for (uint8_t j = 0; j < K - 2; j += 4) {
				u_new[i][j] = u1[i][j] + w * (u2[i][j] - u1[i][j]);
				u_new[i][j + 1] = u1[i][j + 1] + w * (u2[i][j + 1] - u1[i][j + 1]);
				u_new[i][j + 2] = u1[i][j + 2] + w * (u2[i][j + 2] - u1[i][j + 2]);
				u_new[i][j + 3] = u1[i][j + 3] + w * (u2[i][j + 3] - u1[i][j + 3]);
			}
			for (uint8_t j = 200; j < 202; j++)
				u_new[i][j] = u1[i][j] + w * (u2[i][j] - u1[i][j]);
		}

		for (uint8_t i = 0; i < N; ++i) {
			for (uint8_t j = 0; j < K - 2; j += 4) {
				errors.push_back(fabs(u_new[i][j] - u1[i][j]));
				errors.push_back(fabs(u_new[i][j + 1] - u1[i][j + 1]));
				errors.push_back(fabs(u_new[i][j + 2] - u1[i][j + 2]));
				errors.push_back(fabs(u_new[i][j + 3] - u1[i][j + 3]));
			}
			for (uint8_t j = 200; j < 202; j++)
				errors.push_back(fabs(u_new[i][j] - u1[i][j]));
		}
		max_error = *std::max_element(errors.begin(), errors.end());
		//if(k%100==0)std::cout<<k<<" "<<max_error<<"\n";
		errors.clear();
		for (uint8_t i = 0; i < N; ++i) {
			for (uint8_t j = 0; j < K - 2; j += 4) {
				u1[i][j] = u_new[i][j];
				u1[i][j + 1] = u_new[i][j + 1];
				u1[i][j + 2] = u_new[i][j + 2];
				u1[i][j + 3] = u_new[i][j + 3];
			}
			for (uint8_t j = 200; j < 202; j++)
				u1[i][j] = u_new[i][j];
		}
		k++;
	} while (max_error > eps);
	clock_t end_s=clock();

 编译后查看运行时间

总结：

在迭代次数和误差不变的情况下，通过使用循环交换、循环不变量外提、循环展开三种技术对程序进行了优化，每一种都可以提高效率，降低运行时间，混合后，更能提升程序的效率，值得推广。

板牙奇奇
关注
[Math] 常见的几种最优化方法
weixin_34162695的博客
08-23 6261
　　我们每个人都会在我们的生活或者工作中遇到各种各样的最优化问题，比如每个企业和个人都要考虑的一个问题“在一定成本下，如何使利润最大化”等。最优化方法是一种数学方法，它是研究在给定约束之下如何寻求某些因素(的量)，以使某一(或某些)指标达到最优的一些学科的总称。随着学习的深入，博主越来越发现最优化方法的重要性，学习和工作中遇到的大多问题都可以建模成一种最优化模型进行求解，比如我们现在学习的机器学习...
并行程序设计基础——非阻塞通信实现Jacobi迭代
欢迎各位大佬光临“寒舍”！请多多指点，共同进步！
08-20 1060
主要介绍了MPI的非阻塞通信模式和重复非阻塞通信，并实现Jacobi迭代
jacobi迭代的并行优化
06-22
并行计算课程作业，实现jacobi迭代的串行优化，主要是一级和二级cache优化。
机器学习中导数最优化方法(基础篇)
weixin_30266885的博客
10-19 665
1. 前言 熟悉机器学习的童鞋都知道，优化方法是其中一个非常重要的话题，最常见的情形就是利用目标函数的导数通过多次迭代来求解无约束最优化问题。实现简单，coding 方便，是训练模型的必备利器之一。这篇博客主要总结一下使用导数的最优化方法的几个基本方法，梳理梳理相关的数学知识，本人也是一边写一边学，如有问题，欢迎指正，共同学习，一起进步。 2. 几个数学概念 1) 梯度（一阶导数...
约束问题的最优化
qq_54809548的博客
04-15 1808
目录 ACO（蚁群算法） 旅行商问题(TSP)优化 ACATSP.m DrawRoute.m run.m citys_data.mat: 我的资源 最小二乘法 L-M法 lmm.m FK.m JFK.m run.m 高斯-牛顿法 GaussNewton.m 外罚函数+粒子群 PSO.m penalty.m initpop.m getH.m getHeq.m fun.m run.m 线性规划问题-单纯形法 linp.m main.m BFGS+乘子法 ..
Jacobi的GPU并行迭代 OpenACC
admin11111111的博客
06-06 841
GPU上并行，用kernels构件加速二重循环，不计算误差
基于QR分解与Jacobi方法的SVD分解
热门推荐
chenaiyanmie的博客
04-24 1万+
基于QR分解方法的SVD分解：矩阵的 SVD 分解并不唯一。主要的并行算法子程序都是基于经典求解矩阵奇异值的串行方法而实现的。基于 QR 迭代求解矩阵奇异值的方法，是求解双对角矩阵所有奇异值最快速的算法，求解出的奇异值可以达到较高的相对精度；分而治之方法，它是求解全部矩阵奇异值和奇异向量速度很快的算法，但是当求解微小的奇异值时，往往不能保证很高的相对精度，对于一般应用来说，分而治之算法求解出的奇异值
最优化方法及其Matlab程序设计
chen030209的博客
09-14 875
最优化方法及MATLAB程序设计学习笔记
非线性优化-GaussNewton-g2o应用详解
sunxp1985的博客
04-04 912
从案例代码，g2o架构，实际算法多维度剖析《视觉slam十四讲》书中g2o拟合曲线例子
jacobi_jacobi_jacobi迭代matlab_
10-01
然后，通过循环调用`jacobi`函数，更新每个变量的值，直到满足收敛条件。 **应用和局限** Jacobi迭代法适用于大型稀疏矩阵，因为它的计算量相对较小。然而，它可能在非对角占优或病态问题上不收敛，或者收敛速度较...
jacobi迭代求解线性方程的程序
11-12
编写 Jacobi 迭代法的程序需要对 C 语言的基本语法、数据类型、控制结构（如循环和条件语句）以及数组有深入理解。C 语言允许直接操作内存，这在处理大量数据的数值计算中是至关重要的。 ** Jacobi 迭代算法 ** 1....
mj.c.tar.gz_Jacobi c
09-22
mj.c 是一个C语言编写的程序，用于实现 Jacobi 方法，这是一种在数值分析中常见的迭代方法，主要用于求解线性方程组。...同时，这个程序也为我们提供了一个研究如何优化迭代算法和数值稳定性的实例。
机器学习中最优化问题总结
anshuai_aw1的博客
02-14 1269
记录一下自己学习机器学习中最优化问题的过程。 一、基础知识 1.1 梯度向量、Jacobian矩阵和Hessian矩阵 这里讨论的三个概念：梯度向量、Jacobian矩阵和Hessian矩阵 它的自变量：x=(x1,x2,⋯&amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;amp;ThinSpace;,xn)T\mathbf{x} = (x_1,x_2,\cdots,x_n)^Tx=(x1​,x2​,⋯,xn​)T 因变量有两种情况： 一维f
C语言 12 函数
天河有尽后为涯，星海无边前作岸。
09-13 1113
程序的入口点就是main函数，只需要将程序代码编写到主函数中就可以运行了，不过这个函数只是由我们来定义，而不是我们来调用。当然，除了主函数之外，一直在使用的printf那么，函数的具体定义是什么呢？函数是完成特定任务的独立程序代码单元。每次要做这个任务时，都要完完整整地将任务的每一行代码都写下来，如果程序中多处都需要执行这个任务，每个地方都完整地写一遍，实在是太臃肿了，有没有一种更好的办法能优化代码呢？
CF C. Nikita and LCM
qq_74190237的博客
09-13 278
思维题.首先先对数组排序,考虑一种特殊情况,那就是如果数组最大值不能被其他数整除,这代表什么?代表其他数里面有某个质因子,最大值里面没有或者数量不管,所有最小公倍数都会大于最大值.那么其他情况的话,就可以确定所有的数都是最大数的因子,那么形成的最小公倍数也一定是最大数的因子,那么就可以枚举一下最大数的因子,就可以了.给长度为n的数组,求出最长的特殊子数组,如果子数组的最小公倍数不在原数组中,那么就是特殊的.
C++与QML的数据交互
qq_50679120的博客
09-11 531
这里没什么好说的，在C++中组合好json后，通过信号发送到qml中，qml中是可以解析json的。通过Q_PROPERTY属性，可以轻松地在QML中读取和修改数据，并将其更新到C++中，先创建一个类。在QML接收，直接通过data.can_id这样的形式去访问结构体中的元素。这两个也没什么好说的，在C++中直接发送信号就好了，在C++中声明并注册一个结构体。然后发送信号到QML中。
C++_map_set详解
最新发布
2302_77644537的博客
09-16 658
set是按照一定次序存储元素的容器在set中，元素的value也标识它(value就是key，类型为T)，并且每个value必须是唯一的。set中的元素不能在容器中修改(元素总是const)，但是可以从容器中插入或删除它们。在内部，set中的元素总是按照其内部比较对象(类型比较)所指示的特定严格的排序准则（默认是小于）进行排序。set在底层是用二叉搜索树(红黑树)实现的。map是关联容器，它按照特定的次序(按照key来比较)存储由键值key和值value组合而成的元素。
C语言从头学59——学习头文件math.h（二）
m0_72128260的博客
09-11 571
继续学习头文件的数学函数。
jacobi迭代法python三次循环
05-29
以下是使用Jacobi迭代法求解线性方程组的Python代码，其中包含三层循环： ```python import numpy as np # 定义线性方程组的系数矩阵A和常数向量b A = np.array([[4, -1, 0, 0], [-1, 4, -1, 0], [0, -1, 4, -1], [0, 0, -1, 3]]) b = np.array([15, 10, 10, 10]) # 定义初始解x和迭代次数n x = np.zeros_like(b) n = 3 # Jacobi迭代法求解线性方程组 for k in range(n): x_new = np.zeros_like(x) for i in range(A.shape[0]): s = 0 for j in range(A.shape[1]): if j != i: s += A[i, j] * x[j] x_new[i] = (b[i] - s) / A[i, i] x = x_new # 输出最终的解 print("The solution of Ax=b is:", x) ``` 在上面的代码中，首先定义了线性方程组的系数矩阵A和常数向量b，然后定义了初始解x和迭代次数n。接下来，使用三层循环进行了Jacobi迭代法的计算，最后输出了求解得到的解。
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值