论文中经常要在图中再加插图，这样可以省下论文的版图，在Letters这类文章中很重要，因为一般限制你最多3-4版。
现在要在y = sin(x)这幅图中插入y = cos(x)和y = sin(x) * cos(x)这两个小图，分别位于原图的右上方和左下方：
>> x = 0 : 0.01 : 2 * pi;
>> y1 = sin(x); y2 = cos(x); y3 = sin(x) .* cos(x);
>> plot(x, y1)
>> axes('position', [0.55 0.65 0.3 0.2])
>> plot(x, y2)
>> axes('position', [0.15 0.15 0.3 0.2])
>> plot(x, y3)
其中axes函数中用'position'参数就是为了在原图中加多一个坐标轴，后的参数说明该坐标轴在原图的位置和其大小。

有时要在MATLAB中调用Maple的函数参与数值计算，像一些MATLAB本身没有的特殊函数，如Laguerre多项式和Hermite多项式等。此时要注意，Maple是会尽可能地返回解析的表达式，而不是浮点数。比如在MATLAB中输入：
>> maple('LaguerreL', 3, 1, 2)
ans =
LaguerreL(3,1,2)
返回的是解析的结果，这种细果参与不了MATLAB的数值计算，我们需要的是浮点数的结果。根据Maple的基础知识可知，你输入的参数是浮点数，那么Maple就给你返回浮点数结果。比如在Maple中输入：
> LaguerreL(3, 1., 2);
                                -1.333333333
即在输入参数1后面加多一点，代表1.0，这便是个浮点数，于是LaguerreL函数给你返回浮点数结果。据此，在MATLAB中这般输入：
>> a = maple('LaguerreL', 3, '1.', 2)
a =
-1.3333333333333333333333333333333
再将a用str2double输入双精度浮点数即可使用了。结总起来，MATLAB中调用Maple的LaguerreL函数的具体函数为：
% LaguerreL.m
function f = LaguerreL(n, k, x)
xn = length(x);
y = zeros(xn, 1);
for i = 1 : xn
    y(i) = str2double(maple('LaguerreL', [num2str(n) '.'], k, x(i)));
end
f = y;
end

 

然而，MATLAB中已提供了函数做这种事情，这就是mfun函数。例如（下面测试不计启动Maple内核时间）：
>> x = 0 : 0.01 : 3;
>> tic, y = LaguerreL(3, 1, x); toc
Elapsed time is 0.560188 seconds.
>> tic, y2 = mfun('LaguerreL', 3, 1, x); toc
Elapsed time is 0.053402 seconds.
即提供矢量输入时，mfun要比自己所编写的函数要快得多。但有时却不能直接用矢量形式输入，要反覆去调用该函数（LaguerreL）时，由于mfun要做许多格式判断和转换，反而令计算速度下降，如：
>> tic, for i = 1 : 300, y1 = LaguerreL(3, 1, 2); end, toc
Elapsed time is 0.408626 seconds.
>> tic, for i = 1 : 300, y2 = mfun('LaguerreL', 3, 1, 2); end, toc
Elapsed time is 0.920374 seconds.

如有一An x An矩阵A和Bn x Bn矩阵B，求其Kronecker tensor product kron(A, B)与(An x Bn) x 1矢量phi的相乘，即计算

y = kron(A, B) * phi.

这个计算当An和（或）Bn很大时，会变得难以计算。例如An = Bn = 100，那么kron(A, B)就是一个10000 x 10000的矩阵。如果不使用稀疏矩阵（假设A和B够稀疏）的储存方法，在MATLAB里（32位）根本表示不出来。在这种情况下，y = kron(A, B) * phi就无法计算了。

在很多情况下，kron(A, B)会变得更大，例如在DMRG的计算中。为了解决这个问题，是有办法在不算出kron(A, B)的情况下算出y = kron(A, B) * phi的。具体的做法参考了R. M. Noack, S. R. Manmana, arxiv:cond-mat/0510321，可把y写成

y = A * (B * psi')'.

其中psi就是phi，只不过是对phi做了个reshape，使其变成psi = reshape(phi, Bn, An)'。注意到 ' 是转置的意思。利用这个技巧，可以轻易在MATLAB中编写程序，不需要算出kron(A, B)，就完成了y = kron(A, B) * phi这个式子的计算。不用这个技巧，就难以提高DMRG中保留的状态数。

处理这个问题的另一个原因是大型矩阵的对角化问题。要知道，如果矩阵太大，我们是无法用一般的对角化手段去处理的（如MATLAB里那eig函数）。只有使用迭代算法，如Lanczos算法，求出该矩阵最大和最小的本征值本征矢。而这类迭代算法中，最繁重的一步计算就是对该大型矩阵H乘上一矢量phi。由于这是一个迭代的过程，H * phi要计算很多很多次。如果H能写成kron(A, B)，利用上述的技巧，能大大减少H * x的计算量，具体分析可看U. Schollwöck, Rev. Mod. Phys., 77, 259 (2005)中的第二节I.2。

昨天又用不太正大光明的手段，搞了个Matlab R2007a，决定“试用”一下。我也不会“试用”太久，就半年左右吧……（接着“试用R2007b”……）

个人认为Matlab R2007a比起上一个版本来说改进不大，其中最重要的改进是增加了多线程的启动开关，可以让Matlab自动决定用多少个线程。Intel MKL那些能多线程化的函数，应该都能多线程运算。另外，Matlab那些element-wise，即逐个元素操作的函数，如sin，log等都可以多线程运算。另外，现在可以用键盘箭头键来移动图中的物件，这功能对我来说有用……

Matlab语句如果太长，就可以用Dot-Dot-Dot，即...来换行，...后面的代码都会忽略，但昨天在解一个方程组时遇到一个奇怪的问题。该方程组很长很长，每一个方程写出来都要三四百个字符，列出方程时免不了要用...换行。因为是用fsolve做数值解，那就要它该方程组写成一个函数，需要返回代入变量值后每个方程的值，列成一个矢量的形式。例如：
function F = myfun(x)
F = [2*x(1) - x(2);
      -x(1) + 1/2 * x(2)];
如果方程很长，那么换行时，比如写成：
function F = myfun(x)
F = [2*x(1)...
    - x(2);
      -x(1) + 1/2 * x(2)];
即第一个方程做了换行，但这种换行格式是错误的，导致方程老是解不出来。正确应写成：
function F = myfun(x)
F = [2*x(1) -...
    x(2);
      -x(1) + 1/2 * x(2)];
即换行符...要加在操作符“-”后面，否则会报错。但如果是这种语句：
function F = myfun(x)
F = 2*x ...
    + 1;
这又是对的，那问题出在那呢？问题在于矩阵中的写法是不同的，...一定要加操作符后，否则就是错误的。

之前讲过，实验室来了五台双核Pentium D 925计算机，这正适合用来做分布式或并行式计算。我打算只调用那些计算机中的一个核参与计算，留下一个核可以让其他人正常地使用该计算机。我们在这里将会介绍Matlab中Distributed Computing Toolbox的基本使用方法，目标是实现简单的分布式计算。
Distributed Computing Toolbox就是分布式计算工具箱，简称DCT，其可以在多台计算机组成的Cluster中实现分布式或并行式计算。简单来说，我们是把一个很繁重的工作，分解成许多小任务，然后分给不同的计算机去处理，最后把计算结果汇总，以达到提高计算效率的目的。
Matlab的做法是这样的：在每台参与计算的计算机中启动一个叫Matlab Distributed Computing Engine的服务，该服务能启动参与计算的worker的Matlab session和管理各台计算机workers的job manager。Job manager对workers进行管理，给workers分配计算任务，接收workers计算后的结果。而你本人就是client，你要把你的工作分解为多个任务，然后把任务给job manager。job manager就会根据workers的多少和空闲情况，适当地把任务分配给workers去做。workers完成任务后，会把结果返回给job manager。当所有workers都完成任务后，你，即是client，便可以从job manager里取回结果。
具体的概念可以参考Matlab的帮助，我们也不能说得很准确。我们在这里只想给出使用Matlab实现分布式计算的简单步骤，以便初学者快速入门。
１、首先第一步要做的，就是令每台要参与计算的计算机组成局域网。比如我有三台计算机，其IP地址分别为192.168.1.101－192.168.1.103，以下简称计算机名为101，102和103。
２、在三台计算机中安装Matlab Distributed Computing Engine（mdce）服务。安装方法为：如Matlab的安装地址为C:\Program Files\MATLAB\R2006b，则Start->Run->cmd到命令行窗口，进入C:\Program Files\MATLAB\R2006b\toolbox\distcomp\bin目录，运行mdce install命令安装mdce服务。接着去控制台->管理工具->服务，查看Matlab Distributed Computing Engine的属性。进入登录页，选择“此帐户”，输入NT AUTHORITY\NetworkService，删除下面的密码，让该服务以NetworkService的形式登入，以便该服务存取共享的映射网络驱动器中的原程序文件。接着便可以启动该服务了。注意以后重新开机，该服务都会启动，当然你可以设置让它手动启动。
３、启动job manager。任一台计算机都可以启动job manager，只要mdce服务启动了即可。比如使用计算机101，在C:\Program Files\MATLAB\R2006b\toolbox\distcomp\bin目录下，运行以下命令：
startjobmanager -name frenseljobm
该命令启动jobmanager，其名字叫frenseljobm，启动地点为计算机101。
４、启动workers。任一台计算机都可以启动workers，只要mdce服务启动了即可。比如使用计算机101，在C:\Program Files\MATLAB\R2006b\toolbox\distcomp\bin目录下，运行以下命令：
startworker -jobmanagerhost 192.168.1.101 -jobmanager frenseljobm -name worker1
此命令指明在计算机192.168.1.101中，启动名为worker1的worker，而该worker受名为frenseljobm的jobmanager管理。就是说来自乡下101的可怜工人worker1，成为万恶的监工frenseljobm的“马仔”了。接着，监工frenseljobm要在不同村102和103中雇用更多的工人worker2、worker3。运行如下的命令：
startworker -jobmanagerhost 192.168.1.101 -jobmanager frenseljobm -name worker2 -remotehost 192.168.1.102
即可在102计算机中启动一个新的，名为worker2的worker，如此类推启动103计算机的worker3。使用nodestatus命令可以查看节点的状态，加上-remotehost可以查看其他节点的状态。
５、如令计算机101为client，即我们的程序在这里编写的。设程序文件位于D:\Matlab_code\testDCT中。共享出文件夹Matlab_code，在文件夹中按工具->映射网络驱动器->令盘符为Z:->文件夹里填\\192.168.1.101\Matlab_code。于是Z:\testDCT便成为放置你程序的地方了。以同样的方法，让计算机102和103都建立映射网络驱动器，令盘符为Z:，文件夹里填\\192.168.1.101\Matlab_code。这时三台机都可以通过Z:\testDCT访问原程序文件。
６、现在便可以进行计算了。这里给出测试的代码。首先写一个函数，模拟我们实际的工作。
% hp.m
function f = hp(m, n)
H1 = zeros(n);
H2 = zeros(n);
for i = 1 : m
    H = H1 + H2;
end
f = H;
end
将此程序hp.m放在D:\Matlab_code\testDCT中。此函数计算n维随机矩阵的加法m次。接着建立另一个m文件，做具体的分布式计算。
% runDCT.m
tic
% 寻找资源，比如jobmanager在什么地方，叫什么名字。
jm = findResource('scheduler', 'type', 'jobmanager', 'name',...
    'frenseljobm', 'LookupURL', '192.168.1.101');
% 使用刚才找到的资源建立一个工作
job = createJob(jm);
% 设置该工作的文件关联，让所有workers都可以找到原程序文件。
set(job, 'PathDependencies', {'Z:\testDCT'})
% 另一种方法，把用到的原程序文件传给所有workers。
% set(job, 'FileDependencies', {'hp.m'})
N = 100;
M = 1000000;
% 建立三个任务，每任务都是算hp(M, N)。
createTask(job, @hp, 1, {M, N});
createTask(job, @hp, 1, {M, N});
createTask(job, @hp, 1, {M, N});
% 提交工作给jobmanager。
submit(job)
% 等待所有workers都把任务做完。
waitForState(job, 'finished')
% 取出计算结果。
results = getAllOutputArguments(job);
toc
同样地，该程序runDCT.m也是放在D:\Matlab_code\testDCT中。该程序计算了三次100维矩阵的加法1000000次，即算了100维矩阵的加法3000000次。如果在单机上运行：
>> tic, a = hp(3000000, 100); toc
Elapsed time is 63.096369 seconds.
而使用三台机作分布式计算时：
>> runDCT
Elapsed time is 24.323556 seconds.
效率有明显的提升。但注意到，当第一次进行分布式计算时，其他几台机要从Z:\testDCT中读取原程序文件，会使得计算速度降低。
总结来说，Matlab的Distributed Computing Toolbox为我们提供了一种简便的分布式或并行式计算的实现方法。以上所写的是为了对DCT具体做法的整个过程做一次简单的介绍，可能很粗糙和存在许多谬误，敬请指正。

上次说到，Matlab R2006a开始使用Intel MKL，只要设好OMP_NUM_THREADS这个环境变量，即可以BLAS Level 3的运算在任意线程数下。然而，尽管设定的线程数和CPU核的数目相同，但这样也并不能保证能提升计算效率。主要的原因在于建立线程也是需要时间的。如果你的任务只要0.0001秒就能算完，但建立线程也要用0.0001秒，那就根本没有必要把该任务多线程化。
麻烦的是，当我们设定好环境变量OMP_NUM_THREADS后启动Matlab，那么这个进程运行的线程数就定下来了，不能中途改变。于是，我们没有办法根据具体问题随时改变使用的线程数，使得在我们的程序中，能提升某部份的效率，但另一部份的效率却可能降低。
我们可以对矩阵乘法做一点测试，在不同的线性数下，看看不同的矩阵大小，其乘法的效率如何。首先编写程序如下：
THREADS = 2;    % 线程数
N = 2000000;    % N / 矩阵大小 = 每种维度的矩阵要做乘法的次数
% 设定测试的矩阵大小
MN = 50;
step = 2;
x = 10 : step : MN;
n = max(size(x));
T = zeros(2, n);    % 用作存放结果
for i = 1 : n
    M1 = zeros(x(i));
    M2 = zeros(x(i));
    M1 = rand(x(i));
    M2 = rand(x(i));
    t = cputime;    % 准备计时
    for j = 1 : N / x(i)    % 令计算次数随矩阵增大而减少
        tmp = M1 * M2;
    end
    T(:, i) = [x(i) (cputime - t) / THREADS];   % 存放结果
end
可以在不同线程下计算，把T存起来。接着可以比较单线程下的计算结果T1和双线程下计算的结果T2，有：
>> T = T1(2, :) ./ T2(2, :);
>> plot(T1(1, :), T, 'k')
可得出下图。说明当矩阵大小约为20到30时，双线程反而令矩阵乘法效率降低。当矩阵更小时，是否多线程效率更好，我不能保证，应该要作更仔细的验算。

Matlab本身是不可以做多线程计算的，即不可像Fortran/C++使用OpenMP那样，可以指明程序的那个部份应该用多少个线程去执行，但却可以使用间接地调用已经多线程化的程序，如Intel MKL中的BLAS Level 3去实现多线程运算。于是我们想到，是否可以把用Fortran/C++写的OpenMP多线程程序，编译成Matlab可以调用的动态链接库MEX文件。通过调用拥有多线程运算能力的MEX文件，提高Matlab在某些应用中的效率。这种想法是可以实现的，这里简述具体做法。

我们这里使用Fortran做例子，C++中的实现方法和Fortran是类似的。我们使用的是Widnows XP、Microsoft Visual Studio .NET 2005、Intel Visual Fortran Compiler 9.1、Matlab R2006b、CPU是Intel Core 2 Duo E6400。首先设置好MEX要用到的一些路径，如include和lib的位置，具体可参考Matlab帮助中的MATLAB -> External Interfaces -> Calling C and Fortran programs from Matlab -> Custum Building on Windows中的Compiling MEX-Files with the Microsoft Visual C++ IDE。只要按着要求的步骤做即可，要注意编写def文件时，如用IVF编译器的，要写成EXPORTS MEXFUNCTION。

现在可以做个小测试。先从Matlab帮助中，找一个最简单的Fortran MEX文件的例子，我使用的是”A First Example — Passing a Scalar”（只要搜寻”Passing a Scalar”，即可找出此例子）。这个例子是通过接口子程序mexFunction，把输入参数从Matlab里调入，同时给定Matlab的输出变量，再把该输入参数通过另一个子程序timestwo乘以2，将结果放入参数y中，最后把y返回到Matlab的输出变量中。就是说，子程序timestwo是该MEX文件的功能部份，只要将该功能部份并行化，即可提高运算效率。但子程序timestwo太简单了，我们没必要做什么。但我们却可以人为地加上一个很“艰巨”的任务，以此测试OpenMP的使用。

我们考虑做两个1000 X 1000矩阵的加法，要执行2000次。这任务看来并不轻松，足够CPU忙一阵子。既然这事情要做2000次，不如让每个线程做1000次，这可以使用OpenMP的SECTIONS Directive实现。现把子程序timestwo改动如下：
    subroutine timestwo(y, x)
        implicit none
        real(8) :: x, y
        integer, parameter :: N = 1000
        integer :: i
        real(8), allocatable :: a(:, :), b(:, :), c(:, :), d(:, :)
        ! dynamically creates storage
        allocate(a(N, N), b(N, N), c(N, N), d(N, N))
        ! initialization
        a = 1.23456789d0
        b = 9.87654321d0
 
        !$OMP PARALLEL
        !$OMP SECTIONS
 
        !$OMP SECTION
        do i = 1, 1000
            c = a + b
        end do
 
        !$OMP SECTION
        do i = 1, 1000
            d = a + b
        end do
 
        !$OMP END SECTIONS
        !$OMP END PARALLEL
 
        y = 2.0 * x
        return
    end subroutine timestwo
这程序相当简单，它并没有影响原本timestwo的功能，而是在其中加上一些无关的计算，测试OpenMP。我们使用了OpenMP的SECTIONS directive，将矩阵加法分为两个部份，每个部份都是只有一个线程执行，这样便可以使CPU的两个核同时工作。

接着设定好OpenMP的编译命令，即可生成MEX文件。把该MEX文件放在Matlab的相关目录下，假设文件名为timestwo.mexw32，在Matlab里运行：
>> tic, a = timestwo(2), toc
a =
     4
Elapsed time is 22.742711 seconds.
>> clear timestwo
可见，做2000次1000 X 1000矩阵的加法，在Intel Core 2 Duo E6400上，要约22.74秒。在Release模式下，需时约12.50秒（在Pentium D 925中，要29.57秒）。要注意使用timestwo后，最好先把它clear了，否则该MEX文件不能更删。

这里还有一点是必须提一下的。我们在子程序timestwo中，矩阵使用了动态分配内存，而不是静态地预先分配好。这是由于动态链接库是不会分配stack空间的，而是调用该dll（MEX）的程序Matlab分配。由于要分配1000 X 1000的double矩阵，而且每个线程都要分配一份，使得Matlab调用该MEX文件时，很可能出现stack overflow错误。要改动Matlab所分配的stack的大小看来不是个好主意，我也不清楚该如何做，那就只好动态地分配内存，使其数据放在heap中，那就可以避免stack overflow。

Matlab从R14开始使用Intel MKL的BLAS，在Windows平台下，设定环境变量BLAS_VERSION=mkl.dll即可使用（默认还是使用ATLAS BLAS，在R2006a版本这个环境变量不用设了，都用MKL）。而且通过MKL，可实现BLAS Level 3的多线程运行。设定环境变量OMP_NUM_THREADS=2，重新启动Matlab，运行昨天的测试程序：
    function test(l)
    N = 1000;
    a = zeros(N, N);
    b = zeros(N, N);
    a(:) = 1.23456789;
    b(:) = 9.87654321;
    tic
    for i = 1 : l
        c = a * b;
    end
    toc
    end
在Matlab的Command Window运行：
>> test(30)
Elapsed time is 6.300522 seconds.
结果不比在Fortran直接调用MKL的要差。用双核的朋友们，别忘了设定环境变量OMP_NUM_THREADS啊！

近来在写一份note，要经常写上Matlab算出来的矩阵，由本我是用Latex来编写，所以需要把Matlab的数值矩阵转成Latex语句。Matlab中看来没有把数值矩阵转成Latex语句的程序，而只有关于符号矩阵的程序（其实是Maple的），所以我就自己写了一个。我的程序可能写得很烂，但还能凑合用一用。

 

% 数值矩阵输出成Latex代码
% A为矩阵，precision为输出精度
% 如输入：
% >> a = [1 2 3; 4 5 6; 7 8 9];
% >> latex2(a)
% \left(
%  \begin{array}{ccc}
%   1 & 2 & 3 \\
%   4 & 5 & 6 \\
%   7 & 8 & 9 \\
%  \end{array}
% \right)
function f = latex2(A, precision)
% 没输入精度参数时，默认精度为小数后4位
if nargin == 1
    precision = '4';
else
    precision = int2str(precision);
end
% 定义单一元素输出格式
out_num = [' %0.' precision 'f &'];
% 用作整数输出判断
z = zeros(1, str2num(precision) + 1);
z(1) = '.';
z(2 : end) = '0';
z = char(z);
% 求矩阵大小
[r c] = size(A);
nc = zeros(1, c);
nc(:) = 99;  % 存放character c
% 生成第一句Latex语句
out = sprintf('\\left(\n\t\\begin{array}{%s}', char(nc));
% 二重循环，用作生成整个矩阵的Latex语句
for i = 1 : r
    out = [out sprintf('\n\t')]; % 换行
    for j = 1 : c
        temp = sprintf(out_num, A(i, j));
        % 小数位皆为零时，把数取整。如1.0001取为1
        dot_position = find(temp == '.');
        if temp(dot_position : end - 2) == z
            temp = temp(1 : dot_position - 1);
            temp = [temp ' &'];
            % 要取整时，如有负号，则必须丢掉
            if temp(2) == '-'
               temp = [temp(1) temp(3 : end)];
            end
        end
        out = [out temp];
    end
    % 丢掉最后的'&'号
    out = out(1 : end - 1);
    % 行末加上'\\'号
    out = [out '\\'];
end
% 加上最后一句结束代码
out = [out sprintf('\n\t\\end{array}\n\\right)')];
f = out;

如果我有一个用C语言写的函数，实现了一个功能，如一个简单的函数：
double add(double x, double y) {
　return x + y;
}
现在我想要在Matlab中使用它，比如输入：
>> a = add(1.1, 2.2)
                                                     3.3000
要得出以上的结果，那应该怎样做呢？
解决方法之一是要通过使用MEX文件，MEX文件使得调用C函数和调用Matlab的内置函数一样方便。MEX文件是由原C代码加上MEX文件专用的接口函数后编译而成的。
可以这样理解，MEX文件实现了一种接口，它把在Matlab中调用函数时输入的自变量通过特定的接口调入了C函数，得出的结果再通过该接口调回Matlab。该特定接口的操作，包含在mexFunction这个函数中，由使用者具体设定。
所以现在我们要写一个包含add和mexFunction的C文件，Matlab调用函数，把函数中的自变量（如上例中的1.1和2.2）传给mexFunction的一个参数，mexFunction把该值传给add，把得出的结果传回给mexFunction的另一个参数，Matlab通过该参数来给出在Matlab语句中调用函数时的输出值（如上例中的a）。
比如该C文件已写好，名为add.c。那么在Matlab中，输入：
>> mex add.c
就能把add.c编译为MEX文件（编译器的设置使用指令mex -setup），在Windows中，MEX文件类型为mexw32，即现在我们得出add.mexw32文件。现在，我们就可以像调用M函数那样调用MEX文件，如上面说到的例子。所以，通过MEX文件，使用C函数就和使用M函数是一样的了。
我们现在来说mexFunction怎样写。
mexFunction的定义为：
void mexFunction(
　　　　　int nlhs,
　　　　　mxArray *plhs[],
　　　　　int nrhs,
　　　　　const mxArray *prhs[]) {
/*....................................*/
}
可以看到，mexFunction是没返回值的，它不是通过返回值把结果传回Matlab的，而是通过对参数plhs的赋值。mexFunction的四个参数皆是说明Matlab调用MEX文件时的具体信息，如这样调用函数时：
>> b = 1.1; c = 2.2;
>> a = add(b, c)
mexFunction四个参数的意思为：
nlhs = 1，说明调用语句左手面（lhs－left hand side）有一个变量，即a。
nrhs = 2，说明调用语句右手面（rhs－right hand side）有两个自变量，即b和c。
plhs是一个数组，其内容为指针，该指针指向数据类型mxArray。因为现在左手面只有一个变量，即该数组只有一个指针，plhs[0]指向的结果会赋值给a。
prhs和plhs类似，因为右手面有两个自变量，即该数组有两个指针，prhs[0]指向了b，prhs[1]指向了c。要注意prhs是const的指针数组，即不能改变其指向内容。
因为Matlab最基本的单元为array，无论是什么类型也好，如有double array、 cell array、 struct array……所以a,b,c都是array，b = 1.1便是一个1x1的double array。而在C语言中，Matlab的array使用mxArray类型来表示。所以就不难明白为什么plhs和prhs都是指向mxArray类型的指针数组。
完整的add.c如下：
// add.c
#include "mex.h" // 使用MEX文件必须包含的头文件

// 执行具体工作的C函数
double add(double x, double y) {
　return x + y;
}

// MEX文件接口函数
void mexFunction(
　int nlhs,
　mxArray *plhs[],
　int nrhs,
　const mxArray *prhs[]) {
　double *a;
　double b, c;
　plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
　a = mxGetPr(plhs[0]);
　b = *(mxGetPr(prhs[0]));
　c = *(mxGetPr(prhs[1]));
　*a = add(b, c);
}

mexFunction的内容是什么意思呢？我们知道，如果这样调用函数时：
>> output = add(1.1, 2.2);
在未涉及具体的计算时，output的值是未知的，是未赋值的。所以在具体的程序中，我们建立一个1x1的实double矩阵（使用mxCreateDoubleMatrix函数，其返回指向刚建立的mxArray的指针），然后令plhs[0]指向它。接着令指针a指向plhs[0]所指向的mxArray的第一个元素（使用mxGetPr函数，返回指向mxArray的首元素的指针）。同样地，我们把prhs[0]和prhs[1]所指向的元素（即1.1和2.2）取出来赋给b和c。于是我们可以把b和c作自变量传给函数add，得出给果赋给指针a所指向的mxArray中的元素。因为a是指向plhs[0]所指向的mxArray的元素，所以最后作输出时，plhs[0]所指向的mxArray赋值给output，则output便是已计算好的结果了。
上面说的一大堆指向这指向那，什么mxArray，初学者肯定都会被弄到头晕眼花了。很抱歉，要搞清楚这些乱糟糟的关系，只有多看多练。
实际上mexFunction是没有这么简单的，我们要对用户的输入自变量的个数和类型进行测试，以确保
输入正确。如在add函数的例子中，用户输入char array便是一种错误了。
从上面的讲述中我们总结出，MEX文件实现了一种接口，把C语言中的计算结果适当地返回给Matlab罢了。当我们已经有用C编写的大型程序时，大可不必在Matlab里重写，只写个接口，做成MEX文件就成了。另外，在Matlab程序中的部份计算瓶颈（如循环），可通过MEX文件用C语言实现，以提高计算速度。

我们经常把语句或函数写在一个M文件里，比如：
% ----------------------
% test.m
a = 1;
b = 2;
% ----------------------
这个M语句文件，或：
% ----------------------
% test2.m
function a = test2(b)
a = b * 2;
% ----------------------
这个M函数。
我们在Command Window里输入：
>> test
即是等于输入了
>> a = 1;
>> b = 2;
而输入：
>> a = test2(2)
a =
     4
便是执行了一个函数。
以上的内容我想大部份人都知道是怎么一回事了，以下说一说P代码文件。
如果在Command Window里输入：
>> pcode test
便会在相应test.m的文件夹里产生了test.p。如果这时在Command Window里输入：
>> test
其实是运行了test.p，而不是test.m。
P文件是对应M文件的一种预解析版本（preparsed version）。因为当你第一次执行M文件时，Matlab需要将其解析（parse）一次（第一次执行后的已解析内容会放入内存作第二次执行时使用，即第二次执行时无需再解析），这无形中增加了执行时间。所以我们就预先作解释，那么以后再使用该M文件时，便会直接执行对应的已解析版本，即P文件。但又因为Matlab的解析速度非常快，一般不用自己作预解析。只有当一些程序要调用到非常多的M文件时，如GUI应用程序时，才会作预解析，以增加以后的调用速度。
如Matlab的当前目录（Current Directory）有test.m文件，作预解析后，又有test.p文件。因为P文件的调用优先级比M文件要高，所以当你调用test时，会作优先选择而调用了test.p。
可以修改test.m的代码为：
% ----------------------
% test.m
a = 3;
b = 4;
% ----------------------
再在Command Window里调用test：
>> test
Warning: P-file C:\Program Files\MATLAB71\work\test.p is older than M-file C:\Program Files\MATLAB71\work\test.m.
C:\Program Files\MATLAB71\work\test.p may be obsolete and may need to be regenerated.
Type "help pcode" for information about generating P-files from M-files.
会出现一个Warning，警告你所调用的P文件比同名的M文件要旧，即表示M文件已被修改了。尽管如此，调用的还是旧的P文件，即得出a = 1和 b = 2。
P文件可以用来作保密代码之用，如果你给别人一个M文件，别人可以打开来看到你所有的代码和算法。如果你的代码不想被别人看到，那可以给他P文件。
pcode函数也可以应用在M函数文件。

Matlab 7.x的编译器版本为4.x，和以前的3.x版本大不一样。主要体现在不会产生完整的C/C++代码以实现独立的可执行程序，而是通过另外一种间接的途径。

为了解决3.x版本编译器不够兼容Matlab程序的问题，在4.x版本的编译器中，编译器实际上只会产生一些“接口函数”（Wrapper），实际上M文件里的内容被加密压缩存放在一个CTF（Component Technology File）的文档里。具体计算工作由一个名为MCR（MATLAB Component Runtime）的共享库集合工具提供。

MCR就好像一个Matlab一样，只不过是集成了Matlab所有的计算方法，所以通过MCR，4.x版本的编译器几乎能编译所有的Matlab程序。尽管能编译几乎所有的Matlab程序，但这带来了一个麻烦的问题。

当你要把编译好的二进制执行文件（如Windows中的exe）发给其他机器上使用时，该机器上必须安装好MCR。然而，MCR的安装文件大小为99.9MB……所以尽管你的程序可能只完成一次加法，但你要发给其他人的机器中使用你的程序时，请一并把99.9MB的MCR附上。

而且，通过MCR的计算，实际上和在Matlab上计算没差别，速度是一样的。所以，4.x编译器的存在只是为了使你的算法保密。当你的算法程序不想被其他使用者看到时，就要用4.x编译器产生exe。

上次讲到用Maple可以把式子显示得比较漂亮,Matlab中也是可以的,要用pretty函数.

pretty函数可以把一个"符号类型"式子变成漂亮的形式,如:

这个式子,只要输入:

就得到上面的效果了!