关于fortran子例行程序用matlab实现出现的问题，如何解决？

大家好，我有一点问题想要请教一下，若是各位看到了且有能力有时间的话可否为我指点一二。
我现在在将一个远古的fortran程序用matlab复现，但是现在在实现一个子例行程序时发现原程序运行后实部与matlab复现程序相同，但是虚部却差了很多，整个数组中只有第一行第一个列的虚部是正确的，想问问大家有没有什么办法。
下面是fortran程序:


```c
        SUBROUTINE CMUL(A,B)                
        DIMENSION AR(4,4),AI(4,4),BR(4,4),BI(4,4),
     +  CR(4,4),CI(4,4)
        COMPLEX A(4,4),B(4,4)
        N=4
        DO 50 I=1,N
        DO 50 J=1,N
        AR(I,J)=REAL(A(I,J))
        AI(I,J)=AIMAG(A(I,J))
        BR(I,J)=REAL(B(I,J))
        BI(I,J)=AIMAG(B(I,J))
50      CONTINUE
        DO 100 I=1,N
        DO 100 J=1,N
        CR(I,J)=0.0
        CI(I,J)=0.0
        DO 60 L=1,N
        P=AR(I,L)*BR(L,J)
        Q=AI(I,L)*BI(L,J)
        S=(AR(I,L)+AI(I,L))*(BR(L,J)+BI(L,J))
        CR(I,J)=CR(I,J)+P-Q
        CI(I,J)=CI(I,J)+S-P-Q
60      CONTINUE
100     CONTINUE
        DO 110 I=1,N
        DO 110 J=1,N
        A(I,J)=CMPLX(CR(I,J),CI(I,J))
110     CONTINUE
        RETURN
        END

下面的是matlab程序


```c
function [A,B]=CMUL(A,B)
AR=zeros(4,4);
AI=zeros(4,4);
BR=zeros(4,4);
BI=zeros(4,4);
CR=zeros(4,4);
CI=zeros(4,4);
N = 4;
for I = 1:N
    for J = 1:N
        AR(I,J)=real(A(I,J));
        AI(I,J)=imag(A(I,J));
        BR(I,J)=real(B(I,J));
        BI(I,J)=imag(B(I,J));
    end
end
for I = 1:N
    for J = 1:N
        CR(I,J)=0.0;
        CI(I,J)=0.0;
        for L = 1:N
            P=AR(I,L)*BR(L,J);
            Q=AI(I,L)*BI(L,J);
            S=(AR(I,L)+AI(I,L))*(BR(L,J)+BI(L,J));
            CR(I,J)=CR(I,J)+P-Q;
            CI(I,J)=CI(I,J)+S-P-Q;
        end
    end
end
for I = 1:N
    for J = 1:N
        A(I,J)=complex(CR(I,J),CI(I,J));
    end
end
end

这里是用到的两个变量

A = [0.00000000000000 + 0.00000000000000i,0.00000000000000 - 8965350.00000000i,0.00000000000000 + 0.00000000000000i,-23809.5238095238 + 0.00000000000000i;
    0.00000000000000 + 0.00000000000000i,0.00000000000000 + 8965350.00000000i,0.00000000000000 + 0.00000000000000i,-23809.5238095238 + 0.00000000000000i;
    0.00000000000000 - 1799950.00000000i,0.00000000000000 + 0.00000000000000i,-23809.5238095238 + 0.00000000000000i,0.00000000000000 + 0.00000000000000i;
    0.00000000000000 + 1799950.00000000i,0.00000000000000 + 0.00000000000000i,-23809.5238095238 + 0.00000000000000i,0.00000000000000 + 0.00000000000000i];
B = [0.00000000000000 + 0.00000000000000i,0.00000000000000 + 0.00000000000000i,-4.12735218388592e-12 + 4.54277018106903e-07i,-4.12735218388592e-12 - 4.54277018106903e-07i;
    -7.29789917994251e-14 + 6.04065359871497e-08i,-7.29789917994251e-14 - 6.04065359871497e-08i,0.00000000000000 + 0.00000000000000i,0.00000000000000 + 0.00000000000000i;
    0.00000000000000 + 0.00000000000000i,0.00000000000000 + 0.00000000000000i,-1.99999999991745e-05 - 1.81710807247761e-10i,-1.99999999991745e-05 + 1.81710807247761e-10i;
    -1.99999999999854e-05 - 2.41626143948717e-11i,-1.99999999999854e-05 + 2.41626143948717e-11i,0.00000000000000 + 0.00000000000000i,0.00000000000000 + 0.00000000000000i];

fortran数组是列优先，你检查一下输入。另外这个就是复矩阵乘法，为啥要写得这么复杂。直接用复数类型就好。

不知道你这个问题是否已经解决, 如果还没有解决的话:

帮你找了个相似的问题, 你可以看下: https://ask.csdn.net/questions/7659230
除此之外, 这篇博客: 从数学的角度来理解傅里叶变换，离散傅里叶变换，傅里叶矩阵以及在matlab里的实现（避开电信号，频率，正弦波等）中的 一维离散傅里叶变换 部分也许能够解决你的问题, 你可以仔细阅读以下内容或者直接跳转源博客中阅读:
其实函数可以被理解为无限维的向量（这里考虑有限向量），只是需要把它们的定义域理解成自然数集合{0,1,2,3… }\{0,1 ,2 ,3\dots\}{0,1,2,3…}。所以上上面的连续傅里叶变换公式实际上也可以被离散称为向量的形式，关键在于怎么把连续的自变量离散成有限的自然数点集。
现在我们用向量 x={x(k)}k=0N−1x=\{x(k)\}_{k=0}^{N-1}x={x(k)}k=0N−1 代表被离散的函数 fff ， y={y(n)}n=0N−1y=\{y(n)\}_{n=0}^{N-1}y={y(n)}n=0N−1 代表被离散的函数 FFF。
再次将傅里叶变换公式抄下来：
F(ω)=F[f(t)]=∫−∞+∞f(t)e−iωtdtF(\omega)=\mathcal{F}[f(t)]=\int_{-\infty}^{+\infty}f(t)e^{-i\omega t}dt F(ω)=F[f(t)]=∫−∞+∞f(t)e−iωtdt
- 注意到需要被离散的自变量有 ttt 和 ω\omegaω :
  t:{0,1,2,…,N−1}t:\{0,1,2,\dots,N-1\}t:{0,1,2,…,N−1};
  ω:{0,2πN,22πN˙,…,(N−1)2πN}\omega:\{0,\dfrac{2 \pi}{N},2 \dot \dfrac{2 \pi}{N},\dots,(N-1)\dfrac{2 \pi}{N}\}ω:{0,N2π,2N2π˙,…,(N−1)N2π}
所以：
1. 一维离散傅里叶公式：
  y(k)=∑n=oN−1x(n)(e−2πNi)kn,k=0,1,…N−1(3) y(k)=\sum_{n=o}^{N-1}x(n)(e^{-\frac{2 \pi}{N}i})^{kn},\quad k=0,1,\dots N-1\quad\tag{3}y(k)=n=o∑N−1x(n)(e−N2πi)kn,k=0,1,…N−1(3)
2. 一维离散逆傅里叶公式：
  x(n)=1n∑k=oN−1y(k)(e2πNi)kn,n=0,1,…N−1(4)x(n)=\dfrac{1}{n}\sum_{k=o}^{N-1}y(k)(e^{\frac{2 \pi}{N}i})^{kn},\quad n=0,1,\dots N-1 \quad\tag{4}x(n)=n1k=o∑N−1y(k)(eN2πi)kn,n=0,1,…N−1(4)
将一维离散傅里叶公式展开成矩阵向量相乘的形式：
(y(0)y(1)⋮y(N−1))=(w−0⋅0w−0⋅1⋯w−0⋅(N−1)w−1⋅1w−1⋅1⋯w−1⋅(N−1)⋮⋮⋱⋮w−(N−1)⋅1w−(N−1)⋅1⋯w−(N−1)⋅(N−1))(x(0)x(1)⋮x(N−1))(5) \left( \begin{matrix} y(0)\\ y(1)\\ \vdots \\y(N-1) \end{matrix} \right)= \left( \begin{matrix} w^{-0\cdot0}& w^{-0\cdot1} &\cdots & w^{-0\cdot(N-1)}\\ w^{-1\cdot1}& w^{-1\cdot1} &\cdots & w^{-1\cdot(N-1)}\\ \vdots&\vdots&\ddots & \vdots\\ w^{-(N-1)\cdot1}& w^{-(N-1)\cdot1} &\cdots & w^{-(N-1)\cdot(N-1)} \end{matrix} \right) \left( \begin{matrix} x(0)\\ x(1)\\ \vdots \\x(N-1) \end{matrix} \right) \quad(5) ⎝⎜⎜⎜⎛y(0)y(1)⋮y(N−1)⎠⎟⎟⎟⎞=⎝⎜⎜⎜⎛w−0⋅0w−1⋅1⋮w−(N−1)⋅1w−0⋅1w−1⋅1⋮w−(N−1)⋅1⋯⋯⋱⋯w−0⋅(N−1)w−1⋅(N−1)⋮w−(N−1)⋅(N−1)⎠⎟⎟⎟⎞⎝⎜⎜⎜⎛x(0)x(1)⋮x(N−1)⎠⎟⎟⎟⎞(5)
- 其中 w=e2πNiw=e^{\frac{2 \pi}{N}i}w=eN2πi.

如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^