纯属娱乐，新年第一贴，举办一次编程语言的“性能大赛”
，本次比赛的参赛选手有：
Java C# C C++ Fortran Python IDL awk
共8位，We can see，既有像Java、C#这种高富帅选手（标准库很完备，饭来张口，粉丝无数），又包含了C、Fortran这样的纯屌丝选手（靠近底层，事事亲为），还有像Python、IDL这样的科研宠儿，而高大全、伟光正的C++自然不会甘落人后的，最后，awk作为“民科”类选手，也堂而皇之地参赛了，当然了，想也知道这货是成不了黑马的，
。言而总之，这场比赛着实令我期待。
前言：娱乐贴，不要太认真了，如果你想直观的看下这几种语言的实际“面目”，大体对这些语言的执行效率有些了解的话，本帖还是很有意思的，但是如果非要跟我较真“为什么不用所有编译器都编译一次，为什么不用编译优化，为什么不在所有操作系统上都试一遍，，，更有甚者，为什么要用官方库，为什么不自己写一个”的话，那就请绕行吧，且不说本人智商紧缺，我又不是富二代，没那么多闲功夫
。。。话虽如此，不过，I think 我使用的编译器/平台，抑或是操作系统，都是比较主流的，所以得到的结果应当是有一定代表性的。此外，我在编写本日志内的各种编程语言的代码时都尽量注意了“效率优先”（对于IDL、python、Java这一类的解释型语言，适当的编程方式会使性能得到显著改善，有时甚至是量级上的提高），我相信我所编写的代码应当属于“运行较快”的方式，以期充分发挥该语言的性能。
再扯几句，I know，对于很多“学院派”来说，这种对比很“轻浮”，很不“公平”，因为语言这东西说到底就是编译器实现，我对比各种语言实际上对比的是各个编译器实现的性能，所以单纯只用每个语言中的一个编译器来做会让这些同学无法接受。但是，我想对于“科研人员”来说，这样的对比“或许”已经足够了，很多“科研人员”都以“解决问题”为第一要务，无视或轻视性能，所以这篇日志应该会给一部分同学一些“启示”吧，I hope so.
比赛内容：
从一个特制的NetCDF格式的数据中，获得“目标”所在的时空位置（经纬度+日期）。这里的目标指的是，“data”值为1的时空位置。（我特制了一个NetCDF的数据，将它放置在了经度为29.875度，纬度为14.925度，时间为740小时的位置，其他位置的数据都是0）
比赛的内容就是遍历所有维度，找到1所在的位置，比得就是谁能更快的遍历所有数据，而这个目标位置仅仅是为了验证程序是否正确地执行了。总之，这个比赛就是做同样的事，我们就比谁更快！
赛前准备：
我简单地制做了一个NetCDF格式的数据，共有经度、纬度、时间三个维度，经度长600（0.025度起，间隔0.05度），纬度长300（0.025度起，间隔0.05度），时间长186（1月1日0时起，间隔4小时）。
至于使用的设备则是一台intel 4核笔记本，实际主频在满负荷运行时可达到2.6GHz
使用ncdump简单看一下我制作的这个数据的文件头：
netcdf test {
dimensions:
lon = 600 ;
lat = 300 ;
time = 186 ;
variables:
float lon(lon) ;
lon:units = “degrees_east” ;
lon:long_name = “longitude” ;
float lat(lat) ;
lat:units = “degrees_north” ;
lat:long_name = “latitude” ;
int time(time) ;
time:units = “hours from the beginning” ;
time:long_name = “time” ;
int data(time, lat, lon) ;
data:scale_factor = 1 ;
data:add_offset = 0 ;
data:description = “only one record is 1, the rest are all zero.” ;
}
此数据只在一个时空位置下的值为1，其余值都为0。
赛前预测：
按我之前对这几种语言的认知，我对他们的性能排序大体如下：Fortran >= C >= C++ >= IDL > C# >= Java > Python > Bash+awk
啰嗦几句，也许有些同学会问为什么不加入matlab呢，因为我个人没来由地非常讨厌matlab（其实我从来没有用过matlab， ）。
万事俱备，比赛马上开始：
（1）首先上场的是“Java”选手
（解说员：作为赛事的”组织者”，我个人是非常喜欢Java的，私以为Java是最“舒适”的编程语言，但是，大家都知道，Java是运行在虚拟机上的解释型语言，在这种比拼性能的场合里，，，可以说是没有任何优势的。）
编译环境：Java 1.7.0
操作系统：Windows 7
处理代码：
/*
* Author: zhb
* Organization: ITPCAS
* E-mail: zhanghongbo@itpcas.ac.cn
*/
import java.io.IOException;
import ucar.ma2.ArrayFloat;
import ucar.ma2.ArrayInt;
import ucar.nc2.dataset.*;
import ucar.nc2.*;
public class test {
public static void main(String[] args) {
long t0, t1;
t0 = System.currentTimeMillis();
for (int i=0; i<10; ++i)
test();
t1 = System.currentTimeMillis();
log(“Elapsed time: ” + (t1-t0)/1000.0 + ” seconds.”);
}
public static void test() {
NetcdfDataset nfd = null;
String fname = “G:/tout/test.nc”;
try {
nfd = NetcdfDataset.openDataset(fname);
} catch (IOException ioe) {
err(“Can’t open ” + fname + ” as a netcdf dataset.” + ioe.getMessage());
return;
}
Dimension dlon = nfd.findDimension(“lon”);
if (dlon == null) {
err(“can’t find dim lon”);
return;
}
//log(“lon: ” + dlon.getLength());
Variable vlon = nfd.findVariable(“lon”);
if (vlon == null) {
err(“can’t find var lon”);
return;
}
if (vlon.getDimensions().size() != 1 || vlon.getDimension(0) != dlon) {
err(“var lon doesn’t match dim lon.”);
return;
}
Dimension dlat = nfd.findDimension(“lat”);
if (dlat == null) {
err(“can’t find dim lat”);
return;
}
//log(“lat: ” + dlat.getLength());
Variable vlat = nfd.findVariable(“lat”);
if (vlat == null) {
err(“can’t find var lat”);
return;
}
if (vlat.getDimensions().size() != 1 || vlat.getDimension(0) != dlat) {
err(“var lat doesn’t match dim lat.”);
return;
}
Dimension dtime = nfd.findDimension(“time”);
if (dtime == null) {
err(“can’t find dim time”);
return;
}
//log(“time: ” + dtime.getLength());
Variable vtime = nfd.findVariable(“time”);
if (vtime == null) {
err(“can’t find var time”);
return;
}
if (vtime.getDimensions().size() != 1 || vtime.getDimension(0) != dtime) {
err(“var time doesn’t match dim time.”);
return;
}
Variable vdata = nfd.findVariable(“data”);
if (vdata == null) {
err(“can’t find var data”);
return;
}
if (vdata.getDimensions().size() != 3) {
err(“dims of data is not 3″);
return;
}
else if(vdata.getDimension(0) != dtime || vdata.getDimension(1) != dlat || vdata.getDimension(2) != dlon) {
err(“dims’ structure of data is unexpected.”);
return;
}
ArrayInt.D3 data = null;
try {
data = (ArrayInt.D3)vdata.read();
} catch (IOException ioe) {
err(“can’t read data.”);
ioe.printStackTrace();
return;
}
ArrayFloat.D1 lon = null;
try {
lon = (ArrayFloat.D1)vlon.read();
} catch (IOException ioe) {
err(“can’t read lon.”);
ioe.printStackTrace();
return;
}
ArrayFloat.D1 lat = null;
try {
lat = (ArrayFloat.D1)vlat.read();
} catch (IOException ioe) {
err(“can’t read lat.”);
ioe.printStackTrace();
return;
}
ArrayInt.D1 time = null;
try {
time = (ArrayInt.D1)vtime.read();
} catch (IOException ioe) {
err(“can’t read time.”);
ioe.printStackTrace();
return;
}
int nt = dtime.getLength();
int ny = dlat.getLength();
int nx = dlon.getLength();
int c = 0;
int[] data2 = (int[])data.get1DJavaArray(int.class);
for (int t = 0; t < nt; ++t) {
for (int y = 0; y < ny; ++y) {
for (int x = 0; x < nx; ++x) {
if (data2[c] == 1) {
log(“Find the target at lon-” + lon.get(x) + ” lat-” + lat.get(y) + ” time-” + time.get(t));
log(“Corresponding index: x-” + x + ” y-” + y + ” time-” + t);
}
++c;
}
}
}
}
public static void err(String msg) {
System.out.println(“Error: ” + msg);
}
public static void log(String msg) {
System.out.println(msg);
}
}
运行效果：
Find the target at lon-29.875 lat-14.925 time-740
Corresponding index: x-597 y-298 time-185
Elapsed time: 7.538 seconds.
为了更准确地获得实际运行时间，上面的结果是让此java程序运行十次得到的，可以看到，该程序正确地找到了“目标”，单次耗时为0.7538秒，在耗时的量级上，还是可以接受的。
（2） next——微软的“宠儿”: C#
编译环境： VS 2008 SP1
操作系统：Windows 7
处理代码：
private void dotest()
{
int ncidp;
int stat = NetCDF.nc_open(@”G:\tout\test.nc”, NetCDF.CreateMode.NC_NOWRITE, out ncidp);
check(stat);
//log(“open ok!”);
int lonID;
int latID;
int timeID;
int dataID;
int[] lon_dims = new int[1];
int[] lat_dims = new int[1];
int[] time_dims = new int[1];
int[] data_dims = new int[3];
int nx;
int ny;
int nt;
int did;
check(NetCDF.nc_inq_varid(ncidp, “lon”, out lonID));
//log(“lonid:” + lonID);
check(NetCDF.nc_inq_vardimid(ncidp, lonID, lon_dims));
check(NetCDF.nc_inq_dimid(ncidp, “lon”, out did));
if (did != lon_dims[0])
{
Console.WriteLine(“Error: var lon doesn’t match dim lon.”);
return;
}
check(NetCDF.nc_inq_dimlen(ncidp, did, out nx));
//log(“x counts:” + nx);
check(NetCDF.nc_inq_varid(ncidp, “lat”, out latID));
//log(“latid:” + latID);
check(NetCDF.nc_inq_vardimid(ncidp, latID, lat_dims));
check(NetCDF.nc_inq_dimid(ncidp, “lat”, out did));
if (did != lat_dims[0])
{
Console.WriteLine(“Error: var lat doesn’t match dim lat.”);
return;
}
check(NetCDF.nc_inq_dimlen(ncidp, did, out ny));
//log(“y counts:” + ny);
check(NetCDF.nc_inq_varid(ncidp, “time”, out timeID));
//log(“timeid:” + timeID);
check(NetCDF.nc_inq_vardimid(ncidp, timeID, time_dims));
check(NetCDF.nc_inq_dimid(ncidp, “time”, out did));
if (did != time_dims[0])
{
Console.WriteLine(“Error: var time doesn’t match dim time.”);
return;
}
check(NetCDF.nc_inq_dimlen(ncidp, did, out nt));
//log(“time counts:” + nt);
check(NetCDF.nc_inq_varid(ncidp, “data”, out dataID));
//log(“dataid:” + dataID);
check(NetCDF.nc_inq_vardimid(ncidp, dataID, data_dims));
if (data_dims[0] != time_dims[0] || data_dims[1] != lat_dims[0] || data_dims[2] != lon_dims[0])
{
Console.WriteLine(“Error: data’s dimids may be wrong.”);
return;
}
int[, ,] data = new int[nt, ny, nx];
float[] lon = new float[nx];
float[] lat = new float[ny];
int[] time = new int[nt];
stat = NetCDF.nc_get_var_float(ncidp, lonID, lon);
check(stat);
//log(“lon ok!”);
stat = NetCDF.nc_get_var_float(ncidp, latID, lat);
check(stat);
//log(“lat ok!”);
stat = NetCDF.nc_get_var_int(ncidp, timeID, time);
check(stat);
//log(“time ok!”);
stat = NetCDF.nc_get_vara_int(ncidp, dataID, new int[3] { 0, 0, 0 }, new int[3] { nt, ny, nx }, data);
check(stat);
//log(“data ok!”);
int[] scale = new int[1];
int[] add = new int[1];
stat = NetCDF.nc_get_att_int(ncidp, dataID, “scale_factor”, scale);
check(stat);
//log(“scale: ” + scale[0]);
stat = NetCDF.nc_get_att_int(ncidp, dataID, “add_offset”, add);
check(stat);
//log(“add: ” + add[0]);
int goal = (1 – add[0]) / scale[0];
for (int t = 0; t < nt; ++t)
{
for (int y = 0; y < ny; ++y)
{
for (int x = 0; x < nx; ++x)
{
if (data[t, y, x] == 1)
{
Console.WriteLine(“Find the target at lon-” + lon[x] + ” lat-” + lat[y] + ” time-” + time[t]);
Console.WriteLine(“Corresponding index: x-” + x + ” y-” + y + ” time-” + t);
}
}
}
}
NetCDF.nc_close(ncidp);
}
运行结果：
Find the target at lon-29.875 lat-14.925 time-740
Corresponding index: x-597 y-298 time-185
Executed time: 4.945628
同样也是运行了10次，结果正确，单次耗时为0.495秒
（3）这是一位值得期待的选手，它历来被认为是“速度最快”、“接近底层”，甚至被Raymond称为高级汇编器，它就是C
编译环境：GCC 4.7.2
操作系统：Fedora 18 32位
处理代码：
/*
* =====================================================================================
* Version: 1.0
* Created: 2014年01月28日 20时05分11秒
* Revision: none
* Compiler: gcc
* Author: zhanghongbo (zhb), zhanghongbo@itpcas.ac.cn
* Organization: ITPCAS
* =====================================================================================
*/
#include <stdio.h>
#include <stdlib.h>
#include <netcdf.h>
#include <time.h>
#include <sys/timeb.h>
#define ERR(e) \
{ \
printf(“Error: %s\n”, nc_strerror(e)); \
return -1; \
}
/* check one-dimension var */
int check(int ncid, const char *name, int *var_id, int *dims, size_t *n)
{
int status = 0;
if ((status = nc_inq_varid(ncid, name, var_id)))
ERR(status);
if ((status = nc_inq_vardimid(ncid, *var_id, dims)))
ERR(status);
int dim_id = -1;
if ((status = nc_inq_dimid(ncid, name, &dim_id)))
ERR(status);
if (dim_id != dims[0]) {
fprintf(stderr, “Error: var %s doesn’t dim %s\n”, name, name);
return -1;
}
if ((status = nc_inq_dimlen(ncid, dim_id, n)))
ERR(status);
return 1;
}
int test()
{
int ncid;
int status;
char fname[] = “../test.nc”;
if ((status = nc_open(fname, NC_NOWRITE, &ncid)))
ERR(status);
int lon_id;
int lat_id;
int time_id;
int data_id;
int lon_dims[1];
int lat_dims[1];
int time_dims[1];
int data_dims[3];
size_t nx, ny, nt;
if (check(ncid, “lon”, &lon_id, lon_dims, &nx) < 0)
return -1;
if (check(ncid, “lat”, &lat_id, lat_dims, &ny) < 0)
return -1;
if (check(ncid, “time”, &time_id, time_dims, &nt) < 0)
return -1;
if ((status = nc_inq_varid(ncid, “data”, &data_id)))
ERR(status);
if ((status = nc_inq_vardimid(ncid, data_id, data_dims)))
ERR(status);
if (data_dims[0] != time_dims[0] || data_dims[1] != lat_dims[0]
|| data_dims[2] != lon_dims[0]) {
fprintf(stderr, “Error: data’s dimids may be wrong.\n”);
return -1;
}
int *data = (int *)calloc(nx * ny * nt, sizeof(int));
if (data == NULL) {
fprintf(stderr, “Error: can’t calloc for data\n”);
return -1;
}
size_t data_start[] = {0, 0, 0};
size_t data_count[] = {nt, ny, nx};
if ((status = nc_get_vara_int(ncid, data_id, data_start, data_count,
data))) {
ERR(status);
}
float *lon = (float *)calloc(nx, sizeof(float));
float *lat = (float *)calloc(ny, sizeof(float));
int *time = (int *)calloc(nt, sizeof(int));
size_t start[] = { 0 };
size_t count[] = { nx };
if ((status = nc_get_vara_float(ncid, lon_id, start, count, lon)))
ERR(status);
count[0] = ny;
if ((status = nc_get_vara_float(ncid, lat_id, start, count, lat)))
ERR(status);
count[0] = nt;
if ((status = nc_get_vara_int(ncid, time_id, start, count, time)))
ERR(status);
int scale[1];
int add[1];
if ((status = nc_get_att_int(ncid, data_id, “scale_factor”, &scale[0])))
ERR(status);
if ((status = nc_get_att_int(ncid, data_id, “add_offset”, &add[0])))
ERR(status);
int t, y, x;
int c = 0;
int target = (1 – add[0]) / scale[0];
for (t = 0; t < nt; ++t)
for (y = 0; y < ny; ++y)
for (x = 0; x < nx; ++x) {
if (data[c] == target) {
printf(“Find the target at lon-%6.3f lat-%6.3f time-%d\n”,
lon[x], lat[y], time[t]);
printf(“Coresponding index: x-%d y-%d t-%d\n”,
x, y, t);
}
++c;
}
free(lon);
free(lat);
free(time);
free(data);
return 0;
}
int DiffMilliSecs(struct timeb tpstart, struct timeb tpend){
int ztime = (tpend.time – tpstart.time) * 1000 + tpend.millitm – tpstart.millitm;
return ztime;
}
int main(void)
{
int i = 0;
struct timeb t0, t1;
ftime(&t0);
for (; i<10; ++i)
test();
ftime(&t1);
printf(“Elapsed time: %4.2f seconds.\n”, DiffMilliSecs(t0, t1) / 1000.f);
return 0;
}
运行结果：
Find the target at lon-29.875 lat-14.925 time-740
Coresponding index: x-597 y-298 t-185
Elapsed time: 2.54 seconds.
果然不负老夫厚望， ，运行十次耗时2.54秒，单次只有0.254秒，此成绩足以傲视群雄了~
果然不负老夫厚望，
，运行十次耗时2.54秒，单次只有0.254秒，此成绩足以傲视群雄了~（4）全能战士C++要登场了~
编译环境：GCC 4.7.2
操作系统：Fedora 18 32位
处理代码：
/*
* =====================================================================================
* Filename: test.cpp
* Description: test nc lib for c++
* Version: 1.0
* Created: 2014年02月02日 14时42分02秒
* Revision: none
* Compiler: gcc
* Author: zhanghongbo (zhb), zhanghongbo@itpcas.ac.cn
* Organization: ITPCAS
* =====================================================================================
*/
#include <iostream>
#include <string>
#include <netcdf>
#include <time.h>
#include <sys/timeb.h>
using namespace std;
using namespace netCDF;
using namespace netCDF::exceptions;
void test();
int DiffMilliSecs(struct timeb tpstart, struct timeb tpend){
int ztime = (tpend.time – tpstart.time) * 1000 + tpend.millitm – tpstart.millitm;
return ztime;
}
int main()
{
struct timeb t0, t1;
ftime(&t0);
for (int i=0; i<10; ++i)
test();
ftime(&t1);
cout << “Elapsed time: ” << DiffMilliSecs(t0, t1) / 1000. << ” seconds” << endl;
return 0;
}
void test() {
char fname[] = “../test.nc”;
try {
NcFile nc = NcFile(fname, NcFile::read);
/* read all dims */
NcDim lon_dim = nc.getDim(“lon”);
if (lon_dim.isNull()) {
cout << “Error: can’t find dim lon” << endl;
return;
}
NcDim lat_dim = nc.getDim(“lat”);
if (lon_dim.isNull()) {
cout << “Error: can’t find dim lat” << endl;
return;
}
NcDim time_dim = nc.getDim(“time”);
if (lon_dim.isNull()) {
cout << “Error: can’t find dim time” << endl;
return;
}
/* check vars’ dims */
NcVar lon_var = nc.getVar(“lon”);
if (lon_var.isNull()) {
cout << “Error: can’t find var lon” << endl;
return;
}
vector<NcDim> lon_vdims = lon_var.getDims();
if (lon_vdims.size() != 1 || lon_vdims.at(0).getName() != “lon”) {
cout << “Error: var lon doesn’t match dim lon” << endl;
return;
}
NcVar lat_var = nc.getVar(“lat”);
if (lat_var.isNull()) {
cout << “Error: can’t find var lat” << endl;
return;
}
vector<NcDim> lat_vdims = lat_var.getDims();
if (lat_vdims.size() != 1 || lat_vdims.at(0).getName() != “lat”) {
cout << “Error: var lat doesn’t match dim lat” << endl;
return;
}
NcVar time_var = nc.getVar(“time”);
if (time_var.isNull()) {
cout << “Error: can’t find var time” << endl;
return;
}
vector<NcDim> time_vdims = time_var.getDims();
if (time_vdims.size() != 1 || time_vdims.at(0).getName() != “time”) {
cout << “Error: var time doesn’t match dim time” << endl;
return;
}
NcVar data_var = nc.getVar(“data”);
if (data_var.isNull()) {
cout << “Error: can’t find var data” << endl;
return;
}
vector<NcDim> data_vdims = data_var.getDims();
if (data_vdims.size() != 3 || data_vdims.at(0).getName() != “time”
|| data_vdims.at(1).getName() != “lat”
|| data_vdims.at(2).getName() != “lon” ) {
cout << “Error: var data doesn’t match dim data” << endl;
return;
}
/* read all vars */
int nx = lon_dim.getSize();
int ny = lat_dim.getSize();
int nt = time_dim.getSize();
float *lon = new float[nx];
float *lat = new float[ny];
int *time = new int[nt];
int *data = new int[nt * nx * ny];
lon_var.getVar(lon);
lat_var.getVar(lat);
time_var.getVar(time);
data_var.getVar(data);
/* read attrs */
NcVarAtt scale_att = data_var.getAtt(“scale_factor”);
if (scale_att.isNull()) {
cout << “Error: can’t find att scale_factor” << endl;
return;
}
int scale;
scale_att.getValues(&scale);
NcVarAtt add_att = data_var.getAtt(“add_offset”);
if (add_att.isNull()) {
cout << “Error: can’t find att add_offset” << endl;
return;
}
int add;
add_att.getValues(&add);
/* find the target */
int target = (1 – add) / scale;
int c = 0;
for (int t=0; t<nt; ++t)
for (int y=0; y<ny; ++y)
for (int x=0; x<nx; ++x) {
if (target == data[c]) {
cout << “Find the target at lon-” << lon[x] <<
” lat-” << lat[y] << ” time-” << time[t] << endl;
cout << “Corresponding index: x-” << x << ” y-“
<< y << ” t-” << t << endl;
}
++c;
}
delete[] lon;
delete[] lat;
delete[] time;
delete[] data;
}
catch (NcException& e) {
e.what();
}
}
运行结果：
Find the target at lon-29.875 lat-14.925 time-740
Find the target at lon-29.875 lat-14.925 time-740
Corresponding index: x-597 y-298 t-185
Elapsed time: 2.572 seconds
运行10次，共耗时2.572秒，单次0.257秒，与C的执行效率相当，而C++在面向对象和泛型编程方面又远远优于C，所以，跟C相比C++可谓是“占尽优势”了。
运行10次，共耗时2.572秒，单次0.257秒，与C的执行效率相当，而C++在面向对象和泛型编程方面又远远优于C，所以，跟C相比C++可谓是“占尽优势”了。
（5） 虽然我个人对Fortran比较“厌恶”，但是也被迫用了近三年了，不得不说，我觉得Fortran是本次夺冠的最大热门，
编译环境：Compaq Visual Fortran 6.6
操作系统：Windows 7
处理代码：
program test
implicit none
integer :: i
real :: t0, t1
call cpu_time( t0 )
do i = 1, 10
call dotest
enddo
call cpu_time( t1 )
write( *, ‘(“Elapsed time: “, f4.1, ” seconds”)’) t1-t0
end program test
subroutine dotest
use netcdf90
implicit none
include ‘netcdf.inc’
character*100, parameter :: FNAME = “G:/tout/test.nc”
integer :: ncidp
integer :: lon_id, lat_id, time_id, data_id
integer, dimension(1) :: lon_dims, lat_dims, time_dims
integer, dimension(3) :: data_dims
integer :: did
integer :: add
integer :: scale
integer :: nx, ny, nt
integer*4, allocatable :: vdata(:, :, :), vtime(:)
real, allocatable :: vlon(:), vlat(:)
integer :: tmp
integer :: goal
integer :: i, j, k
call check( nf90_open(FNAME, mode = nf90_nowrite, ncid = ncidp) )
call check( nf90_inq_varid(ncidp, “lon”, lon_id) )
call check( nf90_inquire_variable(ncidp, lon_id, dimids = lon_dims) )
call check( nf90_inq_dimid(ncidp, “lon”, dimid = did) )
if (did /= lon_dims(1)) then
write (*, *) “var lon doesn’t match dim lon.”
stop ‘Stopped’
endif
call check( nf90_inquire_dimension(ncidp, did, len = nx) )
call check( nf90_inq_varid(ncidp, “lat”, lat_id) )
call check( nf90_inquire_variable(ncidp, lat_id, dimids = lat_dims) )
call check( nf90_inq_dimid(ncidp, “lat”, dimid = did) )
if (did /= lat_dims(1)) then
write (*, *) “var lat doesn’t match dim lat.”
stop ‘Stopped’
endif
call check( nf90_inquire_dimension(ncidp, did, len = ny) )
call check( nf90_inq_varid(ncidp, “time”, time_id) )
call check( nf90_inquire_variable(ncidp, time_id, dimids = time_dims) )
call check( nf90_inq_dimid(ncidp, “time”, dimid = did) )
if (did /= time_dims(1)) then
write (*, *) “var time doesn’t match dim time.”
stop ‘Stopped’
endif
call check( nf90_inquire_dimension(ncidp, did, len = nt) )
call check( nf90_inq_varid(ncidp, “data”, data_id) )
call check( nf90_inquire_variable(ncidp, data_id, dimids = data_dims) )
if (data_dims(1) /= lon_dims(1) .OR. data_dims(2) /= lat_dims(1) &
.OR. data_dims(3) /= time_dims(1)) then
write (*, *) “var data doesn’t match dim data.”
stop ‘Stopped’
endif
allocate( vdata( 1:nx, 1:ny, 1:nt) )
call check( nf90_get_att(ncidp, data_id, “add_offset”, add) );
call check( nf90_get_att(ncidp, data_id, “scale_factor”, scale) );
call check( nf90_get_var(ncidp, data_id, vdata, (/ 1, 1, 1 /), (/ nx, ny, nt /)) );
allocate( vlon( 1:nx ) )
call check( nf90_get_var(ncidp, lon_id, vlon, (/ 1 /), (/ nx /)) );
allocate( vlat( 1:ny ) )
call check( nf90_get_var(ncidp, lat_id, vlat, (/ 1 /), (/ ny /)) );
allocate( vtime( 1:nt ) )
call check( nf90_get_var(ncidp, time_id, vtime, (/ 1 /), (/ nt /)) );
goal = (1 – add) / scale
do i = 1, nt
do j = 1, ny
do k = 1, nx
if (vdata(k, j, i) == 1) then
write( *, ‘(“Find the target at lon-“, f6.3, ” lat-“, f6.3, ” time-“, I3)’) &
vlon(k), vlat(j), vtime(i)
write( *, ‘(“Corresponding index: x-“, I3, ” y-“, I3, ” t-“, I3)’) &
k, j, i
endif
enddo
enddo
enddo
DEALLOCATE(vlon, vlat, vtime, vdata)
call check(nf90_close(ncidp))
contains
subroutine check (status)
integer, intent(in) :: status
if (status /= nf90_noerr) then
print *, trim( nf90_strerror(status) )
stop “Stopped”
endif
end subroutine check
end subroutine dotest
运行结果：
Find the target at lon-29.875 lat-14.925 time-740
Corresponding index: x-598 y-299 t-186
Elapsed time: 4.9 seconds
同样是运行10次，结果正确，单次耗时为0.49秒，结果显然没有我想象中的快，可能受我编写的代码影响较大（比如我没有选择直接读入整个变量的方式），如果是这样，我只能对Fortran的广大粉丝说I’m sorry了。
（6）终于要轮到近几年在科研圈子里非常火热的Python了，python是标准的脚本语言，代码清晰易读，易学易用，但是效率低下也是出了名的~
解释环境： Python 2.7
操作系统：Windows 7
处理代码：
# Author: zhb
# Contact: zhanghongbo@itpcas.ac.cn
import numpy as np
import netCDF4 as nC
import time
# test reading nc data
def test() :
# open specified nc file
fname = ‘G:/tout/test.nc’
try:
r = nC.Dataset(fname, ‘r’, format=’NETCDF’)
except:
print “Error: can’t open” + fname
return
# read all dims
try:
lon_dim = r.dimensions[‘lon’]
except:
print “Error: no dim called ‘lon’ is found!”
return
try:
lat_dim = r.dimensions[‘lat’]
except:
print “Error: no dim called ‘lat’ is found!”
return
try:
time_dim = r.dimensions[‘time’]
except:
print “Error: no dim called ‘time’ is found!”
return
# read all vars
try:
lon_var = r.variables[‘lon’]
except:
print “Error: no var called ‘lon’ is found!”
return
try:
lat_var = r.variables[‘lat’]
except:
print “Error: no var called ‘lat’ is found!”
return
try:
time_var = r.variables[‘time’]
except:
print “Error: no var called ‘time’ is found!”
return
# test whether var and dim is matching
if (lon_var.dimensions[0] != ‘lon’ or len(lon_var.dimensions) != 1):
print “Error: var lon and dim lon dose not match.”
return
if (lat_var.dimensions[0] != ‘lat’ or len(lat_var.dimensions) != 1):
print “Error: var lat and dim lat dose not match.”
return
if (time_var.dimensions[0] != ‘time’ or len(time_var.dimensions) != 1):
print “Error: var time and dim time dose not match.”
return
# read data
data_var = r.variables[‘data’]
dims = data_var.dimensions
if (len(dims) != 3 or dims[0] != ‘time’ or \
dims[1] != ‘lat’ or dims[2] != ‘lon’):
print “Error: dims of data are not as expected.”
return
# read atts
try:
scale = data_var.getncattr(‘scale_factor’)
except:
print “Error: can’t find the attribute called ‘scale_factor’.”
return
try:
add = data_var.getncattr(‘add_offset’)
except:
print “Error: can’t find the attribute called ‘add_offset’.”
return
# find the target
nx = len(lon_dim)
ny = len(lat_dim)
nt = len(time_dim)
data = data_var[:, :, :]
target = (1 – add) / scale
ti = 0;
yi = 0;
xi = 0;
for t in data:
yi = 0
for y in t:
xi = 0
for x in y:
if (x == target):
print “Find the target at lon-” + str(lon_var[xi]) \
+ ” lat-” + str(lat_var[yi]) + ” time-” + str(time_var[ti])
print “Corresponding index: x-” + str(xi) + ” y-” + str(yi) \
+ ” t-” + str(ti)
xi += 1
yi += 1
ti += 1
r.close();
return
t0 = time.time()
for i in range(0, 1):
test()
t1 = time.time()
print “Elapsed time: ” + str(t1-t0) + ” seconds”
运行结果：
Find the target at lon-29.875 lat-14.925 time-740
Corresponding index: x-597 y-298 t-185
Elapsed time: 10.0220000744 seconds
单次运行了10秒，
，其实也还在预料之中了，虽然有点让人失望。当然了，如果是因为我的代码“羞辱”了python同学的话，那烦请路过的大婶一定要告诉我，
，其实也还在预料之中了，虽然有点让人失望。当然了，如果是因为我的代码“羞辱”了python同学的话，那烦请路过的大婶一定要告诉我，
（7）神秘的IDL要登场了（当然了，神秘是对我而言的），我隐隐有一种要爆出黑马的感觉
编译环境： IDL 8.0.1
操作系统： Windows 7
处理代码：
PRO test_nc
t0 = SYSTIME(1)
for i=1, 10 do begin
test
endfor
t1 = SYSTIME(1)
print, FORMAT = ‘(“Elapsed time: “, F5.2, ” seconds”)’, t1-t0
END
PRO test
; open netcdf data
fname = ‘G:/tout/test.nc’
fid = NCDF_OPEN(fname, /NOWRITE)
; get all dims’ ids
lon_dimid = NCDF_DIMID(fid, ‘lon’)
if lon_dimid eq -1 then begin
print, “Error: can’t find dim lon”
return
endif
lat_dimid = NCDF_DIMID(fid, ‘lat’)
if lat_dimid eq -1 then begin
print, “Error: can’t find dim lat”
return
endif
time_dimid = NCDF_DIMID(fid, “time”)
if time_dimid eq -1 then begin
print, “Error: can’t find dim time”
return
endif
; get all vars’ ids
lon_varid = NCDF_VARID(fid, ‘lon’)
if lon_varid eq -1 then begin
print, “Error: can’t find var lon”
return
endif
lon_varinq = NCDF_VARINQ(fid, lon_varid)
if lon_varinq.Dim ne lon_dimid then begin
print, “Error: dim lon does not match var lon”
return
endif
lat_varid = NCDF_VARID(fid, ‘lat’)
if lat_varid eq -1 then begin
print, “Error: can’t find var lat”
return
endif
lat_varinq = NCDF_VARINQ(fid, lat_varid)
if lat_varinq.Dim ne lat_dimid then begin
print, “Error: dim lat does not match var lat”
return
endif
time_varid = NCDF_VARID(fid, ‘time’)
if time_varid eq -1 then begin
print, “Error: can’t find var time”
return
endif
time_varinq = NCDF_VARINQ(fid, time_varid)
if time_varinq.Dim ne time_dimid then begin
print, “Error: dim time does not match var time”
return
endif
data_varid = NCDF_VARID(fid, ‘data’)
if data_varid eq -1 then begin
print, “Error: can’t find var data”
return
endif
data_varinq = NCDF_VARINQ(fid, data_varid)
if data_varinq.Ndims ne 3 or data_varinq.Dim[0] ne lon_dimid $
or data_varinq.Dim[1] ne lat_dimid or data_varinq.Dim[2] ne time_dimid $
then begin
print, “Error: data’s dims are not structured as expected.”
return
endif
; read all vars
NCDF_VARGET, fid, lon_varid, lon
NCDF_VARGET, fid, lat_varid, lat
NCDF_VARGET, fid, time_varid, time
NCDF_VARGET, fid, data_varid, data
; read attrs
NCDF_ATTGET, fid, data_varid, ‘scale_factor’, scale
NCDF_ATTGET, fid, data_varid, ‘add_offset’, add
goal = (1 – add) / scale
sub = where(data eq 1)
if sub eq -1 then begin
print, ‘Error: no target is found!’
return
endif
ng = size(sub, /N_ELEMENTS)
NCDF_DIMINQ, fid, lon_dimid, name, nx
NCDF_DIMINQ, fid, lat_dimid, name, ny
NCDF_DIMINQ, fid, time_dimid, name, nt
for i=0, ng-1 do begin
ind = sub[i]
ti = ind / (nx * ny)
yi = (ind mod (nx * ny)) / nx
xi = ind – ti * (nx * ny) – yi * nx
print, FORMAT = ‘(“Find the target at lon-“, F6.3, ” lat-“, F6.3, ” time-“, I3)’, $
lon[xi], lat[yi], time[ti]
print, FORMAT = ‘(“Corresponding index: x-“, I3, ” y-“, I3, ” t-“, I3)’, xi, yi, ti
endfor
END
运行结果：
Find the target at lon-29.875 lat-14.925 time-740
Corresponding index: x-597 y-298 t-185
Elapsed time: 2.35 seconds
令人吃惊的2.35秒，这是10次的时间，也就是说单次耗时不到0.3秒，这当然跟代码的编写方式有关，我自认为在编写上述代码时有效利用了idl的优势，毫无疑问，idl已经是冠军了。
（8）awk是来打酱油的
解释环境：gawk
操作系统：Fedora 18 32位
处理代码：
Bash代码
t0=`date +%s`
ncdump test.nc | awk -f test_nc.awk
t1=`date +%s`
time=$((t1 – t0))
echo ‘Elapsed time: ‘${time}’ seconds’
awk代码
BEGIN {
mode = 0
FS = “,”
target = 1
nt = 186
nx = 600
ny = 300
c = 0
found = 0
}
function has_goal() {
for (i=1; i<=NF; ++i) {
if ($i == target) {
found = 1
return i
}
}
return 0
}
{
if ($0 ~ /^ lon =/) {
mode = 1
next
}
else if ($0 ~ /^ lat =/) {
mode = 2
next
}
else if ($0 ~ /^ time =/) {
mode = 3
next
}
else if ($0 ~ /^ data =/) {
mode = 4
n = split($1, group, “=”)
$1 = group[2]
ind = has_goal()
if ($NF !~ /^ *$/) {
print “rare case!”
if (ind > 0) {
c += ind
}
else {
c += NF
}
}
else {
if (ind > 0) {
c += ind
}
else {
c += NF – 1
}
}
next
}
if (mode == 4) {
if (found == 0) {
if ($NF !~ /^ *$/) {
n = split($NF, group, “;”)
$NF = group[1]
}
ind = has_goal()
if (ind > 0) {
c += ind
found = 1
}
else {
if ($NF !~ /^ *$/) {
c += NF
}
else {
c += NF – 1
}
}
}
next
}
}
END {
if (found == 1) {
ti = int(c / (nx * ny))
yi = int((c – ti * (nx * ny)) / nx)
xi = c – ti * (nx * ny) – yi * nx
++ti
++yi
print “Find the target at x-” xi ” y-” yi ” t-” ti
}
}
运行结果：
Find the target at x-598 y-299 t-186
Elapsed time: 19 seconds
运行结果正确，用时19秒，我觉得还不错了。
实际上这种处理方式是受我回家前帮所里一位师姐处理数据时的启发，她当时就是直接使用ncdump将nc数据直接输出成cdl格式，然后用fortran去读，实话说，我的第一反应是“好2”，但是后来觉得这种“土法炼钢”也无可厚非啊，完全可以“解决问题”了。所以，我就想如果用处理文本的脚本语言来做会是什么效率呢，then，我第一时间想到的就是awk了。
运行结果正确，用时19秒，我觉得还不错了。
实际上这种处理方式是受我回家前帮所里一位师姐处理数据时的启发，她当时就是直接使用ncdump将nc数据直接输出成cdl格式，然后用fortran去读，实话说，我的第一反应是“好2”，但是后来觉得这种“土法炼钢”也无可厚非啊，完全可以“解决问题”了。所以，我就想如果用处理文本的脚本语言来做会是什么效率呢，then，我第一时间想到的就是awk了。
赛后总结：
比赛成绩汇总如下，
| Java | C# | C | C++ | Fortran | Pyhon | IDL | awk |
| 0.754s | 0.495s | 0.254s | 0.257s | 0.490s | 10.022s | 0.235s | 19.000s |
IDL“意外”夺魁，awk意料之中的垫底，除了python有些慢之外，其他语言的表现都在可接受范围之内（for me）。
总结陈词：
在处理nc数据上，在我所使用的这些编译器之中，IDL、C和C++是最快的，其后是C#和Fortran，Java稍慢一些，Python在执行效率上表现较差。
在处理nc数据上，在我所使用的这些编译器之中，IDL、C和C++是最快的，其后是C#和Fortran，Java稍慢一些，Python在执行效率上表现较差。
需要指出的是，本次“比赛”中肯定有不“全面”、不“严谨”的地方，至少包括以下几方面：（1）没有尝试所有编译器以及编译优化；（2）不同的操作系统会对程序执行的效率产生影响；（3）netcdf库版本不同，执行效率可能（一定）也会不同；（4）编写的代码可能（一定）不是最高效的。
后记：In fact，此文有一定的误导性，在编写IDL的代码时，我“故意”使用了where函数，当时已然料到IDL会以“极快”的速度完成任务了，其实我个人觉得，IDL这种（类似）脚本语言更像是命令行界面的软件，它的优势在于“专有函数”，那些函数很有可能是使用C或者汇编语言实现的，这种“内置函数”在编写时定然是十分看重效率的，而函数之外，如果用户要自己写底层函数的话，效率想来也会较低，所以，从开发角度来看，还是C、C++更合适，我个人觉得这是很显然的（或许，我又“大放厥词”了）~
。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。
Language shapes the way we think, and determines what we can think about.
– B.L.Whorf
您好！我想咨询一下您关于C#读取nc格式数据的问题，您所用的netcdf.dll库是从哪里下载的？可否告知，谢谢！