最近导师给了个Maple实现并行计算的任务，具体就是求一个马尔可夫链中的状态转移概率矩阵，由于每个转移概率的公式比较复杂（不只是因为矩阵本身的尺寸比较大，矩阵大小范围从100x100到2000x2000）。在使用串行计算的情况下，计算一个200x200的矩阵，所需时间约为330s，2000x2000时耗时5天不止。

Maple Programming Guide中给出了两种并行模型，一种称为Task Programming Model（任务编程模型），它的并行是线程级的（在一个进程下运行多个线程实现并行），线程间的通信通过共享内存实现。另一种称为Grid Programming Model（不知怎么翻译），实现进程级的并行（也可以实现在一个网络上的分布计算，参看Grid:-Setup帮助），进程间的通信通过Send/Receive函数实现。

就我们的问题而言，使用Task Programming Model是最自然的想法，各个线程只需共享一个用于存储矩阵的全局变量，然后分别计算矩阵的不同区域并写入，线程之间不需要交流什么信息，也不需要考虑各个线程的不同运行顺序（因为他们计算不同的区域）。我认为这个想法没有错，但在实际操作上遇到了一些问题。调试之后我认为是计算概率时使用了numbcomb函数来求组合数，而这个函数并非是Thread-safe的，也就是说多个线程同时调用这个函数会产生预料之外的问题。

这里可能可以通过重写一个numbcomb实现，但是我还是决定改用Grid来做。用Grid时一定要明白每个节点有一个单独的运行环境，这里的运行环境包括执行任务的函数，环境变量（包括像Digits这样的变量）。Launch任务时需要将这些环境发送给每一个节点。节点运行结束后，只会接受第一个节点（Node 0）返回的结果。所以一般来说Node 0都是负责任务的调度的调度，以及接受处理其他节点的运算结果。

具体代码：

CalculateTransitionMatrix := proc()    uses Grid;    global m,r,CV,ROU,P;    local this_node,nodes_num,current_job,j,k,result,count;        this_node := MyNode();    nodes_num := NumNodes();        if this_node <> 0 then         while true do            current_job := Receive(0);            if current_job < 0 then break; end if;            result := NULL;            for j from 1 to m+r+2 do                result := result, CalculateTransitionProbability(current_job-1,j-1,m,r,CV,ROU);            end do;            result := this_node, current_job, result;            Send(0, result);         end do;    else        P := Matrix(m+r+2,m+r+2);        for current_job from 1 to nodes_num-1 do            Send(current_job,current_job);        end do;        count := 0;        while count < nodes_num - 1 do            result := Receive();                        j := result[1];            Send(j,current_job);                        if current_job = -1 then                 count := count + 1;            end if;                        if current_job <> -1 then                 current_job := current_job + 1;            end if;                        if current_job > m+r+2 then                current_job := -1;            end if;            j := result[2];            for k from 1 to m+r+2 do                P[j,k] := result[k+2];            end do;        end do;    end if;    Barrier();end proc;Grid:-Launch(CalculateTransitionMatrix,imports=['CalculateTransitionProbability','m','r','CV','ROU','Px','Py','Digits'],exports=['P']);

其中CalculateTransitionProbability计算矩阵中的一个元素。整个模型很简单除了Node 0之外的其他节点每次计算矩阵中的一行。Node 0负责接受其他节点计算出的结果，然后要计算的下一行的行号发给该节点，没有多余的计算任务时发送-1。接收到-1的节点终止任务，当0号节点给其他每个节点都发送了-1之后结束。

Node 0中的几个if逻辑稍微有些混乱，应该按照Code Complete中的几个原则改写一下的。之前有个改好的版本，但没在这台机子上。

另外值得注意的一点事，Grid:-Launch()在不指定numnodes的情况下，默认启动的节点数为当前CPU的逻辑核心数，比如在E3-1240 V2的机器上会启动8个节点进行运算。