【并行程序设计实践小组作业】并行矩阵求逆代码说明

在 C++ 中，vector 是一个十分有用的容器。它能够像容器一样存放各种类型的对象，简单地说，vector是一个能够存放任意类型的动态数组，能够增加和压缩数据。

vector 是同一种类型的对象的集合，每个对象都有一个对应的整数索引值。和 string 对象一样，标准库将负责管理与存储元素相关的内存。我们把 vector 称为容器，是因为它可以包含其他对象。一个容器中的所有对象都必须是同一种类型的。

向量（vector）是一个封装了动态大小数组的顺序容器（Sequence Container）。跟任意其它类型容器一样，它能够存放各种类型的对象。可以简单的认为，向量是一个能够存放任意类型的动态数组。

生成随机double型数组（一维向量）

vector generate_random_double_array(int size)
{
    vector 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;
}

生成随机double型矩阵（二维向量）

vector> generate_random_double_matrix(int size)
{
    // declare 2D vector 声明二维向量
    vector> matrix;
    for (int i = 0; i < size; i++)
    {
        // 调用函数 generate_random_double_array
        vector row = generate_random_double_array(size);
        // 向量 matrix 尾部增加一个元素 row
        matrix.push_back(row);
    }
    return matrix;
}

打一个flag

vector> generate_identity(int size)
{
    vector> I;
    for (int i = 0; i < size; i++)
    {
        vector 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;
}

并行求逆

vector> generate_inverse_parallel(vector> input_matrix, int threads)
{
    int size = input_matrix.size();
    // 设置主对角线的标记
    vector> 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;
}

串行求逆

vector> generate_inverse_serial(vector> input_matrix)
{
    signed int size = input_matrix.size();
    signed int i = 0;
    vector> 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;
}

打印矩阵

void print_matrix(vector> 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;
}

主函数

int main()
{
    // 初始化随机数发生器
    srand(time(NULL));
    double dtime;
    vector> matrix;
    vector> matrix_serial; //串行
    vector> 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;
}

1. 求出分别使用 2～6 个线程进行并行时，实现矩阵求逆的时间。
2. 求出串行矩阵求逆的时间。
3. 使用10阶矩阵测试串行程序与并行程序正确性。
4. 分别用100阶和500阶测试并行程序的加速比与效率。
5. 数据自己写程序生成后写出独立测试文件，不需要手动输入。