并行矩阵求逆代码说明
1. vector
在 C++ 中,vector 是一个十分有用的容器。它能够像容器一样存放各种类型的对象,简单地说,vector是一个能够存放任意类型的动态数组,能够增加和压缩数据。
vector 是同一种类型的对象的集合,每个对象都有一个对应的整数索引值。和 string 对象一样,标准库将负责管理与存储元素相关的内存。我们把 vector 称为容器,是因为它可以包含其他对象。一个容器中的所有对象都必须是同一种类型的。
向量(vector)是一个封装了动态大小数组的顺序容器(Sequence Container)。跟任意其它类型容器一样,它能够存放各种类型的对象。可以简单的认为,向量是一个能够存放任意类型的动态数组。
2. generate_random_double_array
生成随机double型数组(一维向量)
vector<double> generate_random_double_array(int size)
{
vector<double> random_array;
for (int i = 0; i < size; i++)
{
// 生成随机数 10~108
double random_number = (rand() % 99) + 10;
// 向量 random_array 尾部增加一个元素 random_number
random_array.push_back(random_number);
}
return random_array;
}
3. generate_random_double_matrix
生成随机double型矩阵(二维向量)
vector<vector<double>> generate_random_double_matrix(int size)
{
// declare 2D vector 声明二维向量
vector<vector<double>> matrix;
for (int i = 0; i < size; i++)
{
// 调用函数 generate_random_double_array
vector<double> row = generate_random_double_array(size);
// 向量 matrix 尾部增加一个元素 row
matrix.push_back(row);
}
return matrix;
}
4. generate_identity
打一个flag
vector<vector<double>> generate_identity(int size)
{
vector<vector<double>> I;
for (int i = 0; i < size; i++)
{
vector<double> row;
for (int j = 0; j < size; j++)
{
// 行号和列号相等时,向量 row 尾部增加一个元素 1
if (i == j)
{
row.push_back(1);
continue;
}
// 向量 row 尾部增加一个元素 0
row.push_back(0);
}
// 向量 I 尾部增加一个元素 row
I.push_back(row);
}
return I;
}
5. generate_inverse_parallel
并行求逆
vector<vector<double>> generate_inverse_parallel(vector<vector<double>> input_matrix, int threads)
{
int size = input_matrix.size();
// 设置主对角线的标记
vector<vector<double>> I = generate_identity(size);
for (int i = 0; i < size; i++)
{
if (input_matrix[i][i] == 0)
{
// swap nearest subsequent row s.t input_matrix[i][i] != 0 after swapping
// 如果交换后 input_matrix[i][i] != 0 ,则交换最近的后一行
for (int j = i + 1; j < size; j++)
{
if (input_matrix[j][i] != 0.0)
{
swap(input_matrix[i], input_matrix[j]);
break;
}
if (j == size - 1)
{
// 输出 “该矩阵的逆矩阵不存在”
cout << "Inverse does not exist for this matrix";
exit(0);
}
}
}
double scale = input_matrix[i][i];
// 设置即将发生的并行区域中的线程数为 threads
omp_set_num_threads(threads);
#pragma omp parallel for
for (int col = 0; col < size; col++)
{
input_matrix[i][col] = input_matrix[i][col] / scale;
I[i][col] = I[i][col] / scale;
}
if (i < size - 1)
{
#pragma omp parallel for
for (int row = i + 1; row < size; row++)
{
double factor = input_matrix[row][i];
for (int col = 0; col < size; col++)
{
input_matrix[row][col] -= factor * input_matrix[i][col];
I[row][col] -= factor * I[i][col];
}
}
}
}
for (int zeroing_col = size - 1; zeroing_col >= 1; zeroing_col--)
{
#pragma omp parallel for
for (int row = zeroing_col - 1; row >= 0; row--)
{
double factor = input_matrix[row][zeroing_col];
for (int col = 0; col < size; col++)
{
input_matrix[row][col] -= factor * input_matrix[zeroing_col][col];
I[row][col] -= factor * I[zeroing_col][col];
}
}
}
return I;
}
6. generate_inverse_serial
串行求逆
vector<vector<double>> generate_inverse_serial(vector<vector<double>> input_matrix)
{
signed int size = input_matrix.size();
signed int i = 0;
vector<vector<double>> I = generate_identity(size);
for (i = 0; i < size; i++)
{
if (input_matrix[i][i] == 0)
{
// swap nearest subsequent row s.t input_matrix[i][i] != 0 after swapping
// 如果交换后 input_matrix[i][i] != 0 ,则交换最近的后一行
for (int j = i + 1; j < size; j++)
{
if (input_matrix[j][i] != 0.0)
{
swap(input_matrix[i], input_matrix[j]);
break;
}
if (j == size - 1)
{
cout << "Inverse does not exist for this matrix";
exit(0);
}
}
}
double scale = input_matrix[i][i];
for (int col = 0; col < size; col++)
{
input_matrix[i][col] = input_matrix[i][col] / scale;
I[i][col] = I[i][col] / scale;
}
if (i < size - 1)
{
for (int row = i + 1; row < size; row++)
{
double factor = input_matrix[row][i];
for (int col = 0; col < size; col++)
{
input_matrix[row][col] -= factor * input_matrix[i][col];
I[row][col] -= factor * I[i][col];
}
}
}
}
for (int zeroing_col = size - 1; zeroing_col >= 1; zeroing_col--)
{
for (int row = zeroing_col - 1; row >= 0; row--)
{
double factor = input_matrix[row][zeroing_col];
for (int col = 0; col < size; col++)
{
input_matrix[row][col] -= factor * input_matrix[zeroing_col][col];
I[row][col] -= factor * I[zeroing_col][col];
}
}
}
return I;
}
7. print_matrix
打印矩阵
void print_matrix(vector<vector<double>> matrix)
{
for (size_t i = 0; i < matrix.size(); i++)
{
cout << endl;
for (size_t j = 0; j < matrix.size(); j++)
{
cout << matrix[i][j] << " ";
}
}
cout << endl
<< endl;
}
8. main
主函数
int main()
{
// 初始化随机数发生器
srand(time(NULL));
double dtime;
vector<vector<double>> matrix;
vector<vector<double>> matrix_serial; //串行
vector<vector<double>> matrix_parallel; //并行
int SIZE;
// checking for accuracy
// 检查准确性
SIZE = 10;
matrix = generate_random_double_matrix(SIZE);
cout << "Input Matrix is:";
print_matrix(matrix);
matrix_parallel = generate_inverse_parallel(matrix, 2);
cout << "Parallel Matrix output is:";
print_matrix(matrix_parallel);
matrix_serial = generate_inverse_serial(matrix);
cout << "Serial Matrix output is:";
print_matrix(matrix_serial);
// changing the size for higher value
// 改变矩阵大小,测试运行时间( 100阶矩阵)
SIZE = 100;
matrix = generate_random_double_matrix(SIZE);
cout << "New Size of Matrix is:" << SIZE << endl;
cout << "Number of threads: 2" << endl;
dtime = omp_get_wtime();
matrix_parallel = generate_inverse_parallel(matrix, 2);
cout << "Time taken in parallel:" << omp_get_wtime() - dtime << endl;
cout << "Number of threads: 4" << endl;
dtime = omp_get_wtime();
matrix_parallel = generate_inverse_parallel(matrix, 4);
cout << "Time taken in parallel:" << omp_get_wtime() - dtime << endl;
cout << "Number of threads: 5" << endl;
dtime = omp_get_wtime();
matrix_parallel = generate_inverse_parallel(matrix, 5);
cout << "Time taken in parallel:" << omp_get_wtime() - dtime << endl;
cout << "Number of threads: 6" << endl;
dtime = omp_get_wtime();
matrix_parallel = generate_inverse_parallel(matrix, 6);
cout << "Time taken in parallel:" << omp_get_wtime() - dtime << endl;
cout << endl;
dtime = omp_get_wtime();
matrix_serial = generate_inverse_serial(matrix);
cout << "Time taken in serial:" << omp_get_wtime() - dtime << endl;
cout << endl;
// 改变矩阵大小,测试运行时间( 500阶矩阵)
SIZE = 500;
matrix = generate_random_double_matrix(SIZE);
cout << "New Size of Matrix is:" << SIZE << endl;
cout << "Number of threads: 2" << endl;
dtime = omp_get_wtime();
matrix_parallel = generate_inverse_parallel(matrix, 2);
cout << "Time taken in parallel:" << omp_get_wtime() - dtime << endl;
cout << "Number of threads: 4" << endl;
dtime = omp_get_wtime();
matrix_parallel = generate_inverse_parallel(matrix, 4);
cout << "Time taken in parallel:" << omp_get_wtime() - dtime << endl;
cout << "Number of threads: 5" << endl;
dtime = omp_get_wtime();
matrix_parallel = generate_inverse_parallel(matrix, 5);
cout << "Time taken in parallel:" << omp_get_wtime() - dtime << endl;
cout << "Number of threads: 6" << endl;
dtime = omp_get_wtime();
matrix_parallel = generate_inverse_parallel(matrix, 6);
cout << "Time taken in parallel:" << omp_get_wtime() - dtime << endl;
cout << endl;
dtime = omp_get_wtime();
matrix_serial = generate_inverse_serial(matrix);
cout << "Time taken in serial:" << omp_get_wtime() - dtime << endl;
cout << endl;
return 0;
}
9. 实现功能
- 求出分别使用 2~6 个线程进行并行时,实现矩阵求逆的时间。
- 求出串行矩阵求逆的时间。
- 使用10阶矩阵测试串行程序与并行程序正确性。
- 分别用100阶和500阶测试并行程序的加速比与效率。
- 数据自己写程序生成后写出独立测试文件,不需要手动输入。