动态分配巨大的3D阵列

Question

我正在一个项目中，我需要创建一个3D数组，一些2D和1D数组。 3D阵列表示空间中的离散坐标，我需要很多关于我的问题的要点。数组大小将在2000 * 2000 * 2000左右。我需要在这些数组中存储“double”值。任何人都可以提出一个有效的方案来实现这个在C？

在此先感谢

/*********************************************************** 
* Copyright Univ. of Texas M.D. Anderson Cancer Center 
* 1992. 
* 
* Some routines modified from Numerical Recipes in C, 
* including error report, array or matrix declaration 
* and releasing. 
****/ 
#include <stdlib.h> 
#include <stdio.h> 
#include <math.h> 
#include <malloc.h> 

/*********************************************************** 
* Report error message to stderr, then exit the program 
* with signal 1. 
****/ 
void nrerror(char error_text[]) 

{ 
    fprintf(stderr,"%s\n",error_text); 
    fprintf(stderr,"...now exiting to system...\n"); 
    exit(1); 
} 

/*********************************************************** 
* Allocate an array with index from nl to nh inclusive. 
* 
* Original matrix and vector from Numerical Recipes in C 
* don't initialize the elements to zero. This will 
* be accomplished by the following functions. 
****/ 
double *AllocVector(short nl, short nh) 
{ 
    double *v; 
    short i; 

    v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double)); 
    if (!v) nrerror("allocation failure in vector()"); 

    v -= nl; 
    for(i=nl;i<=nh;i++) v[i] = 0.0; /* init. */ 
    return v; 
} 

/*********************************************************** 
* Allocate a matrix with row index from nrl to nrh 
* inclusive, and column index from ncl to nch 
* inclusive. 
****/ 
double **AllocMatrix(short nrl,short nrh, 
        short ncl,short nch) 
{ 
    short i,j; 
    double **m; 

    m=(double **) malloc((unsigned) (nrh-nrl+1) 
         *sizeof(double*)); 
    if (!m) nrerror("allocation failure 1 in matrix()"); 
    m -= nrl; 

    for(i=nrl;i<=nrh;i++) { 
    m[i]=(double *) malloc((unsigned) (nch-ncl+1) 
         *sizeof(double)); 
    if (!m[i]) nrerror("allocation failure 2 in matrix()"); 
    m[i] -= ncl; 
    } 

    for(i=nrl;i<=nrh;i++) 
    for(j=ncl;j<=nch;j++) m[i][j] = 0.0; 
    return m; 
} 

/*********************************************************** 
* Allocate a 3D array with x index from nxl to nxh 
* inclusive, y index from nyl to nyh and z index from nzl to nzh 
* inclusive. 
****/ 
double ***Alloc3D(short nxl,short nxh, 
        short nyl,short nyh, 
         short nzl, short nzh) 
{ 
    double ***t; 
    short i,j,k; 


    t=(double ***) malloc((unsigned) (nxh-nxl+1)*sizeof(double **)); 
    if (!t) nrerror("allocation failure 1 in matrix()"); 
    t -= nxl; 

    for(i=nxl;i<=nxh;i++) { 
    t[i]=(double **) malloc((unsigned) (nyh-nyl+1)*sizeof(double *)); 
    if (!t[i]) nrerror("allocation failure 2 in matrix()"); 
    t[i] -= nyl; 
    for(j=nyl;j<=nyh;j++) { 
     t[i][j]=(double *) malloc((unsigned) (nzh-nzl+1)*sizeof(double)); 
     if (!t[i][j]) nrerror("allocation failure 3 in matrix()"); 
     t[i][j] -= nzl;} 

    } 

    for(i=nxl;i<=nxh;i++) 
    for(j=nyl;j<=nyh;j++) 
     for(k=nzl; k<=nzh;k++) t[i][j][k] = 0.0; 
    return t; 
} 
/*********************************************************** 
*Index to 3D array. 
****/ 
long index(int x, int y, int z, int Size) 
{ 
    return (Size*Size*x + Size*y + z); 
} 
/*********************************************************** 
* Release the memory. 
****/ 
void FreeVector(double *v,short nl,short nh) 
{ 
    free((char*) (v+nl)); 
} 

/*********************************************************** 
* Release the memory. 
****/ 
void FreeMatrix(double **m,short nrl,short nrh, 
       short ncl,short nch) 
{ 
    short i; 

    for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl)); 
    free((char*) (m+nrl)); 
} 


/*********************************************************** 
* Release the memory. 
****/ 
void Free3D(double ***t,short nxl,short nxh, 
       short nyl,short nyh, short nzl, short nzh) 
{ 
    short i,j; 

    for(i=nxh;i>=nxl;i--) 
    {for(j=nyl;j>=nyl;j--) free((char*) (t[i][j]+nzl)); 
    free((char*) (t[i]+nyl)); 
    } 
    free((char*) (t+nxl)); 
} 

*********************************************************************************** 

void InitOutputData(InputStruct In_Parm, OutStruct * Out_Ptr) 
{ 
    short nz = In_Parm.nz; 
    short nr = In_Parm.nr; 
    short na = In_Parm.na; 
    short nl = In_Parm.num_layers; 
    short size = nr/2*nr/2*nz; 
    /* remember to use nl+2 because of 2 for ambient. */ 

    if(nz<=0 || nr<=0 || na<=0 || nl<=0) 
    nrerror("Wrong grid parameters.\n"); 

    /* Init pure numbers. */ 
    Out_Ptr->Rsp = 0.0; 
    Out_Ptr->Rd = 0.0; 
    Out_Ptr->A = 0.0; 
    Out_Ptr->Tt = 0.0; 

    /* Allocate the arrays and the matrices. */ 
    //Out_Ptr->Rd_ra = AllocMatrix(0,nr-1,0,na-1); 
    //Out_Ptr->Rd_r = AllocVector(0,nr-1); 
    //Out_Ptr->Rd_a = AllocVector(0,na-1); 

    Out_Ptr->A_xyz1 = AllocVector(0,size-1); 
    Out_Ptr->A_xyz2 = AllocVector(0,size-1); 
    Out_Ptr->A_xyz3 = AllocVector(0,size-1); 
    Out_Ptr->A_xyz4 = AllocVector(0,size-1); 
    Out_Ptr->A_xz = AllocMatrix(0,nr-1,0,nz-1); 
    Out_Ptr->A_z = AllocVector(0,nz-1); 
    Out_Ptr->A_l = AllocVector(0,nl+1); 

    Out_Ptr->Tt_ra = AllocMatrix(0,nr-1,0,na-1); 
    Out_Ptr->Tt_r = AllocVector(0,nr-1); 
    Out_Ptr->Tt_a = AllocVector(0,na-1); 
} 

以上是用于分配的阵列和函数来调用它们的代码。失败的调用是'Out_Ptr-> A_xyz1 = AllocVector（0，size-1）;'当尺寸大于约。 7000.

Answer 1

如果它们在运行时的固定大小（或者至少是固定的最大大小），并且它们大于物理RAM，那么您也可以使用内存映射文件。访问至少与RAM +交换一样快，并且您可以将数据序列化到磁盘上免费。您还可以映射总体上大于您的地址空间的映射文件的区域（即窗口）的视图。

如果您需要大量的单元格，因为您需要某些区域的高细节，但不是统一的，可以考虑八叉树。然后，您可以在某些部分存储逐渐更高的分辨率，并且您可以选择重新排序以优化对3D附近区域的访问 - 这在CFD或医疗成像等领域非常常见。

Answer 2

这是mmap的完美用例！您的阵列为64 GB - 太大而无法一次装入RAM。幸运的是，我们可以迫使内核为我们完成所有繁重的工作。

这里的一些（当然是未经测试）代码：

#include <sys/mman.h> 
#include <sys/stat.h> 
#include <fcntl.h> 
#include <unistd.h> 
#include <string.h> 
#include <stdio.h> 

#define ARRAY_SIZE ((off_t)2000*2000*2000*8) 

static inline off_t idx(off_t x, off_t y, off_t z) 
{ 
    return 2000*2000*x + 2000*y + z; 
} 

int main() 
{ 
    int fd = open("my_huge_array.dat", O_RDWR | O_CREAT | O_TRUNC); 

    // We have to write to the end of the file for the size to be set properly. 
    lseek(fd, ARRAY_SIZE - 1, SEEK_SET); 
    write(fd, "", 1); 

    double* my_huge_array = mmap(NULL, ARRAY_SIZE, PROT_READ | PROT_WRITE, MAP_SHARED, fd, 0); 

    // Now play with the array! 
    my_huge_array[idx(4, 7, 9)] = 12; 
    my_huge_array[idx(0, 0, 0)] = my_huge_array[idx(4, 7, 9)] * 2; 

    // Close it up. Don't leak your fd's! 
    munmap(my_huge_array, ARRAY_SIZE); 
    close(fd); 
    remove("my_huge_array.dat"); 
    return 0; 
}