我有一个调用 Fortran DLL 进行数值计算的 C# 应用程序。已确定 DLL 调用至少比从 Fortran 控制台应用程序调用相同计算慢 10 倍,因此我已开始将代码移植到 C#。移植过程是逐行复制代码并修复语法。所以 C# 和 Fortran 看起来非常相似。Fortran 中的公共块成为类中的字段。在移植了一些核心例程后,我开始测试,发现双精度 C# 代码比双精度 Fortran 慢 30 倍(比单精度 Fortran 慢 50 倍)。我循环测试代码 100 倍,以最大限度地减少 C# JIT 编译器开销的影响。
该代码使用了复杂的算术和复杂的函数,例如 SQRT 和 LOG10。我提供了一个 C# 结构来处理该数学。我怀疑这就是问题所在,但是 VS2010 Pro 中似乎没有分析器,所以我不确定。
Fortran 编译器是 Intel 的最新版本。我没有对 C# 或 Fortran 代码进行任何特殊优化。我比较了使用发布版本的时间。我只有一个两核 CPU,所以并行化可能不会有太大帮助。
我可以使用一些关于如何加速此代码的建议。
这是我的一种方法。C# 看起来就像我没有编写的 Fortran。
public void Forwx(double rshldr, double rbed, double[][] resdep, double[]toolparm)
{
var zex = new Complex[70];
var zey = new Complex[70];
var px2 = new Complex[70];
var px4 = new Complex[70];
var rlt = new double[2];
var trsmt = new double[2];
var fr = new double[2];
var dnc = -0.02;
var factr = 26332.65;
var rh2 = Math.Max(0.1, rbed);
var rh1 = Math.Max(0.1, rshldr);
const double e1 = 1.0;
const double e2 = 1.0;
const double er = 0.818730753077982;
const double re = 1.0 / er;
var ii = Complex.I;
const double pi = Math.PI;
const double eps0 = 8.854e-12;
const double amu0 = 4.0e-7 * pi;
for (var ktool = 3; ktool <= 6; ktool++)
{
if (ktool == 3) // Integrated 2MHz
{
dnc = -0.02;
rlt[0] = 0.2794;
rlt[1] = -0.2794;
trsmt[0] = 0.904875;
trsmt[1] = -0.904875;
fr[0] = 2000.0;
factr = 26332.65;
}
if (ktool == 4) // Integrated 400kHz
{
dnc = -0.02;
rlt[0] = 0.2794;
rlt[1] = -0.2794;
trsmt[0] = 0.904875;
trsmt[1] = -0.904875;
fr[0] = 400.0;
factr = 26811.866;
}
if (ktool == 5) // Option 5 20kHz
{
dnc = -0.1;
rlt[0] = 0.0;
rlt[1] = 0.0;
trsmt[0] = 5.75;
trsmt[1] = 5.75;
fr[0] = 20.0;
factr = 26811.866 * 2.516 * toolparm[1] / 0.28e8;
}
if (ktool == 6) // Option 6 50kHz
{
dnc = -0.1;
rlt[0] = 0.0;
rlt[1] = 0.0;
trsmt[0] = 5.75;
trsmt[1] = 5.75;
fr[0] = 50.0;
factr = 26811.866 * 6.291 * toolparm[2] / 0.7e8;
}
var r1 = trsmt[0] - rlt[0];
var r2 = trsmt[0] - rlt[1];
var omega = 2000.0 * pi * fr[0];
var k12 = omega*amu0*(omega*e1*eps0 + ii/rh1);
var k22 = omega*amu0*(omega*e2*eps0 + ii/rh2);
var krat = (k22 - k12)/k12;
for (var iz = 0; iz < 601; iz++)
{
var recx1 = new Complex(0.0, 0.0);
var rx1 = new Complex(0.0, 0.0);
var recy1 = new Complex(0.0, 0.0);
var ry1 = new Complex(0.0, 0.0);
var lam = new Complex(3.01517934056e-04 / (Math.Pow(er, 5) * r1));
Complex c1;
Complex c2;
for (var i = 0; i < 70; i++)
{
if (iz == 0)
{
lam = lam * re;
var lam2 = lam * lam;
var p11 = lam2 - k12;
var p1 = Complex.Sqrt(p11);
var p22 = lam2 - k22;
var p2 = Complex.Sqrt(p22);
zex[i] = Complex.Exp(dnc * p2);
zey[i] = Complex.Exp(dnc * p1);
c1 = p2 * k12;
c2 = p1 * k22;
var t3 = lam / p2;
var t2 = t3 * (c1 - c2) / (c1 + c2);
var q2 = lam * krat * (t2 + t3) / (p1 + p2);
px2[i] = (lam2 * q2 + lam * p2 * t2);
px4[i] = px2[i];
}
else
{
px2[i] = px2[i] * zex[i];
px4[i] = px4[i] * zey[i];
}
recx1 = recx1 + a1[i] * px2[i];
recy1 = recy1 + a1[i] * px4[i];
rx1 = rx1 + px2[i] * as1i[i];
ry1 = ry1 + px4[i] * as1i[i];
}
if (ktool <= 4)
{
c1 = recx1*r1;
c2 = rx1*r2;
c2 = c2 - Math.Pow(r1/r2,3)*c1;
resdep[12 - ktool][iz + 600] = c2.Re*factr;
c1 = recy1*r1;
c2 = ry1*r2;
c2 = c2 - Math.Pow(r1 / r2,3) * c1;
resdep[12 - ktool][600 - iz] = c2.Re*factr;
}
else
{
c1 = recx1*r1;
//c2 = rx1*r2;
//c2 = c2 - Math.Pow(r1 / r2,3) * c1;
resdep[ktool + 5][iz + 600] = c1.Re * factr;
c1 = recy1*r1;
//c2 = ry1*r2;
//c2 = c2 - Math.Pow(r1 / r2,3) * c1;
resdep[ktool + 5][600 - iz] = c1.Re * factr;
}
}
}
}
以下是复杂结构中的一些方法。
public static Complex Sqrt(double x)
{
return x >= 0 ? new Complex(Math.Sqrt(x)) : new Complex(0, Math.Sqrt(-x));
}
public static Complex Exp(Complex z)
{
return new Complex(Math.Exp(z.Re) * Math.Cos(z.Im), Math.Exp(z.Re) * Math.Sin(z.Im));
}
public static Complex Log(Complex z)
{
return new Complex(Math.Log(Abs(z)), Arg(z));
}
这是复杂结构的一部分。
public struct Complex
{
private readonly double _re;
private readonly double _im;
#region Properties
public double Re
{
get { return _re; }
//set { re = value; }
}
public double Im
{
get { return _im; }
//set { im = value; }
}
public static Complex I
{
get { return new Complex(0.0, 1.0); }
}
public static Complex Zero
{
get { return new Complex(0.0, 0.0); }
}
public static Complex One
{
get { return new Complex(1.0, 0.0); }
}
#endregion
#region constructors
public Complex(double x)
{
_re = x;
_im = 0.0;
}
public Complex(Complex z)
{
_re = z.Re;
_im = z.Im;
}
public Complex(double x, double y) //constructor
{
_re = x;
_im = y;
}
#endregion
}