计算科学-第一次课程设计

第一次课程设计：

1. 把矩阵向量相乘示例中的for_each()，for_each_n()，generate()，transform()等stl库函数调用改写成普通for循环。测试算例。

2. 保留上面的stl库函数调用，但添加execution的执行策略，测试并行的效果。

3. 调用utilities_sc.h中Timer计时器，给程序添加计时功能。

4. 改写给矩阵赋初值的方式，允许非方阵的矩阵向量相乘。

按照发布的实验报告模板，提交课程设计。谢谢。

一、代码理解或课程设计原理

1. for_each()，for_each_n()，generate()，transform()函数

四个函数的结构都如下：

for_each( _first, _last, func() )
其中：
_first和_last都为正向迭代器，[first, last) 是用于赋值的范围；
func()是函数，可以是lambda函数或者一般函数函数，或者函数对象。

这些函数都可以改为如下所示的for循环体：

for( i=_first;i<_last;i++)
{
   func();
}

2.execution——标准库的并行算法

能够标定了程序应该串行、并行，还是与向量化并行。

std::execution::seq: 串行执行
std::execution::par: 多线程并行执行
std::execution::par_unseq: 多个线程上并行，可以循环交叉，也能使用SIMD(单指令多数据)
std::execution::par或std::execution::par_unseq允许算法并行或向量化并行。。

二、课程设计过程及结果
1.把矩阵向量相乘示例中的for_each()，for_each_n()，generate()，transform()等stl库函数调用改写成普通for循环
（1） for_each函数与for循环的等价：

for_each(begin(vec4), end(vec4), [&norm1](ZPLX &x){ norm1 = norm1< AbsT(x) ? AbsT(x):norm1; } );   

for(size_t i=0; i<m; i++)
{
       ZPLX x=vec4[i];
       norm1 = norm1< AbsT(x) ? AbsT(x):norm1;
}

（2） for_each_n函数与for循环的等价：

for_each_n(std::execution::seq, arr.begin()+i*cols, i+1, [&urd, &mre](auto& e){ e = urd(mre); } );  

for (size_t j = 0; j < i+1; j++) 
{
      arr[i*cols + j] = urd(mre);
}

（3）transform函数与for循环的等价：

transform(begin(out), end(out), begin(vr), begin(out), minus<>());

for(size_t i=0; i<out.size(); ++i) out[i] -= vr[i];

（4） generate函数与for循环的等价：

generate(begin(vec1), end(vec1), [&n, i=-1]() mutable {i++; return ZPLX(0.8*(i-0.5*n), i); } );

for (size_t i=0; i< m; i++) 
{
    vec1[i]= ZPLX(0.8*(i-0.5*n), i);
}

2.保留上面的stl库函数调用，但添加execution的执行策略，测试并行的效果。
如下，分别加入execution的不同执行策略。

std::execution::seq: 串行执行
std::execution::par: 多线程并行执行
std::execution::par_unseq: 多个线程上并行，可以循环交叉，也能使用SIMD(单指令多数

1.transform(std::execution::seq,begin(out), end(out), 
2.      begin(vr), begin(out), minus<>());
3.
4.generate(std::execution::seq,begin(vec1), end(vec1), 
5.      [&n, i=-1]() mutable {i++; return ZPLX(0.8*(i-0.5*n), i); } );
6.
7.for_each(std::execution::seq,begin(vec4), end(vec4),
8.      [&norm1](ZPLX &x){ norm1 = norm1< AbsT(x) ? AbsT(x):norm1; } );
9.
10.for_each_n(std::execution::seq, arr.begin()+i*cols, i+1, 
11.      [&urd, &mre](auto& e){ e = urd(mre); } );  

3.调用utilities_sc.h中Timer计时器，给程序添加计时功能。

1.#include "utilities_sc.h"     //加入头文件
2.    ......
3.
4.    Timer t;
5.    vec4 -= vec2;
6.
7.    norm1 = 0.0;
8.    for_each(std::execution::seq,begin(vec4), end(vec4), [&norm1](ZPLX &x){ norm1 = norm1< AbsT(x) ? AbsT(x):norm1; } );
9. 
10.    norm2 = norm_L2(vec4);
11.    t.printDiff("for_each: ");
12.    ......

如上，通过调用Timer类中的printDiff函数计时。

4.改写给矩阵赋初值的方式，允许非方阵的矩阵向量相乘。
矩阵初始化头文件中只能初始化方阵。示例中sparse_mat/danse_mat是nm维的矩阵，vec1长度为m，即是m1维的矩阵，两者相乘结果是n1维的。我将其更改为nm的矩阵与m*q的矩阵相乘。
更改如下：
（1）将矩阵赋初值的方式由方阵改为非方阵：
原init_matrix.h中通过上下对角矩阵赋值的方式只能对方阵赋值，我将其更改为对全矩阵赋值。

1.    for (size_t i = 0; i < rows; i++) 
2.        for_each_n(std::execution::seq, arr.begin()+i*cols, cols, 
3.            [&urd, &mre](auto& e){ e = urd(mre); } ); 

（2）除nm维矩阵外，再引入一个mq维的矩阵与其相乘：
在两个.cc文件中，除了n，m外，引入整数q，用vec1[mq]存储这个mq维的矩阵的元素，并为其赋值。

1.size_t n{20}, m{25}, q{20};
2.vector<double> vec1(m*q);
3.generate(begin(vec1), end(vec1), [&n, i=-1]() mutable {i++; return (0.8*(i-0.5*n));});

（3）将nm维的矩阵与mq维的矩阵相乘：
各个vector的作用如下：

1.vector<double> vec1(m*q);   //存储初始的m*q维的矩阵
2.vector<double> vec1_1(m);   //每次将vec1中的一列共m个数据传给vec1_1用于乘法运算
3.vector<double> vec2(n*q);   //存储sparse_mat*vec1计算结果，共n*q个数据
4.vector<double> vec4(n*q);   //存储dense_mat*vec1计算结果，共n*q个数据
5.vector<double> vec2_1;  //存储每次让sparse_mat与vec1_1相乘得到的结果，后续要传递到vec2中
6.vector<double> vec4_1; //存储每次让dense_mat与vec_1相乘得到的结果，后续要传递到vec4中

用vec1_1[m]每次从vec1中取出一列的数据与sparse_mat/danse_mat相乘，并将计算结果存储在vec2_1[m]/vec4_1[m]中，然后将vec2_1[m]/vec4_1[m]中的数据映射到vec2[nq]/vec4[nq]相应列中，最后打印输出vec2/vec4。

1.   
2.    size_t n{20}, m{25},q{10};
3.
4.    SparseMat<double> sparse_mat(n, m, n*m);
5.    DenseMat<double> dense_mat(n, m);
6.
7.    vector<double> vec1(m*q);    //存储一个包含m*q个数据的矩阵
8.    generate(begin(vec1), end(vec1), [&n, i=-1]() mutable {i++; return (0.8*(i-0.5*n)); } );
9.    vector<double> vec1_1(m);   //每次将vec1中的一列m个数据传给vec1_1
10.    vector<double> vec2(n*q);   //存储sparse_mat*vec1的计算结果，共n*q个数据
11.    vector<double> vec4(n*q);   //存储dense_mat*vec1的计算结果，共n*q个数据
12.    
13.    for(size_t x = 0; x < q;x++)
14.    {
15.        for(size_t y = 0;y < m;y++)
16.        {
17.            vec1_1[y] = vec1[x+y*q];//每次将vec1中的一列传给vec1_1
18.        }
19.
20.        vector<double> vec2_1 = sparse_mat*vec1_1;
21.        vector<double> vec4_1 = dense_mat*vec1_1;
22.
23.        for(size_t z = 0;z < n;z++)
24.        {
25.            vec2[z*q+x] = vec2_1[z];     //每次将结果的一列传入vec2
26.            vec4[z*q+x] = vec4_1[z];    //每次将结果的一列传入vec4
27.        }
28.    }

三、结果分析

1.把矩阵向量相乘示例中的for_each()，for_each_n()，generate()，transform()等stl库函数调用改写成普通for循环。
见上，略。

2.保留上面的stl库函数调用，但添加execution的执行策略，测试并行的效果。
对一些程序，std::execution::seq 串行执行结果正确，而std::execution::par多线程并行执行和std::execution::par_unseq由于并行使数据传入不同而导致结果不同。

3.调用utilities_sc.h中Timer计时器，给程序添加计时功能。
见上，略。

4.改写给矩阵赋初值的方式，允许非方阵的矩阵向量相乘。
详细结果可以运行代码，略。

详细报告及修改后的代码包见：
https://download.csdn.net/download/weixin_50836014/85007999