C++读取带头部的Fortran输出二进制矩阵文件

核心前提说明 #

Fortran DIRECT访问模式的无格式二进制文件没有记录边界的额外标记，所有记录按固定长度连续存储，且Linux下gfortran默认以字节为RECL（记录长度）单位。读取的关键是严格匹配Fortran写入时的记录长度、数据类型和存储顺序。

步骤1：明确Fortran写入规则 #

先从原Fortran代码中确认以下信息：

• 打开文件时指定的RECL值（每条记录的字节长度）
• 7条头部记录各自的内容、数据类型（如字符、整数、双精度浮点数）及对应字节数
• 矩阵数据的存储方式（列主序是Fortran默认，写入时是单条大记录还是分多条列记录）

步骤2：高效读取头部信息 #

用C++低级IO（std::ifstream::read()）直接按二进制块读取，避免字符串流等中间层。通过结构体映射Fortran的记录结构，禁止结构体字节填充确保匹配。

示例代码框架 #

#include <iostream>
#include <fstream>
#include <vector>
#include <cstdint>
#include <algorithm>
#include <cctype>
#include <fcntl.h> // 用于posix_fadvise优化

using namespace std;

// 按Fortran头部记录定义结构体，禁止字节填充
#pragma pack(push, 1)
// 第1条记录：256字节描述文本
struct HeaderRec1 {
    char desc[256];
};
// 第2条记录：行列数（4字节整数）
struct HeaderRec2 {
    int32_t rows;
    int32_t cols;
};
// 第3条记录：双精度参数
struct HeaderRec3 {
    double alpha;
    double beta;
};
// 补充剩余4条头部记录的结构体
#pragma pack(pop)

int main() {
    // 关闭C++与C标准IO同步，加速读取
    ios::sync_with_stdio(false);
    cin.tie(nullptr);

    ifstream infile("data.bin", ios::binary | ios::in);
    if (!infile.is_open()) {
        cerr << "文件打开失败" << endl;
        return 1;
    }

    // 告诉操作系统预读文件，优化顺序读取效率
    int fd = infile.native_handle();
    posix_fadvise(fd, 0, 0, POSIX_FADV_SEQUENTIAL);

    // 读取第1条头部记录
    HeaderRec1 rec1;
    infile.read(reinterpret_cast<char*>(&rec1), sizeof(rec1));
    if (!infile) { cerr << "读取头部记录1失败" << endl; return 1; }
    // 清理Fortran字符串末尾的空格
    string desc(rec1.desc);
    desc.erase(find_if(desc.rbegin(), desc.rend(), [](int ch) {
        return !isspace(ch);
    }).base(), desc.end());
    cout << "描述：" << desc << endl;

    // 读取第2条头部记录（获取行列数）
    HeaderRec2 rec2;
    infile.read(reinterpret_cast<char*>(&rec2), sizeof(rec2));
    if (!infile) { cerr << "读取头部记录2失败" << endl; return 1; }
    const int32_t rows = rec2.rows;
    const int32_t cols = rec2.cols;
    cout << "行数：" << rows << "，列数：" << cols << endl;

    // 读取第3条及剩余4条头部记录
    HeaderRec3 rec3;
    infile.read(reinterpret_cast<char*>(&rec3), sizeof(rec3));
    if (!infile) { cerr << "读取头部记录3失败" << endl; return 1; }
    cout << "Alpha：" << rec3.alpha << "，Beta：" << rec3.beta << endl;
    // ... 读取剩余4条头部记录的代码 ...

    // 后续读取矩阵数据

步骤3：高效读取大矩阵（20-60G） #

矩阵读取的核心是减少IO次数，优先大块读取，内存不足时分块处理。

方案1：内存足够时一次性读取 #

// 预分配矩阵内存（列主序，与Fortran一致）
    const size_t total_elements = static_cast<size_t>(rows) * cols;
    vector<double> mat(total_elements);

    // 一次性读取整个矩阵
    infile.read(reinterpret_cast<char*>(mat.data()), total_elements * sizeof(double));
    if (!infile) { cerr << "读取矩阵失败" << endl; return 1; }

    // 访问方式：mat[i + rows*j] 对应Fortran的 mat(i+1, j+1)（Fortran是1索引）
    cout << "第一个元素：" << mat[0] << endl;
    cout << "最后一个元素：" << mat[total_elements - 1] << endl;

方案2：内存不足时分块读取 #

// 每次读取100列（可根据内存调整块大小）
    const size_t block_cols = 100;
    const size_t block_size = static_cast<size_t>(rows) * block_cols;
    vector<double> block(block_size);
    vector<double> mat(total_elements);

    for (size_t j = 0; j < cols; j += block_cols) {
        const size_t current_cols = min(block_cols, cols - j);
        const size_t current_size = static_cast<size_t>(rows) * current_cols;
        
        infile.read(reinterpret_cast<char*>(block.data()), current_size * sizeof(double));
        if (!infile) { cerr << "读取块" << j/block_cols << "失败" << endl; return 1; }

        // 将块数据复制到矩阵对应位置
        copy(block.begin(), block.begin() + current_size, mat.begin() + static_cast<size_t>(rows)*j);
    }

编译与优化 #

用g++开启最高优化编译：

g++ -O3 -std=c++11 read_data.cpp -o read_data

内容的提问来源于stack exchange，提问作者Adam Ciesielski