動態分配巨大的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; 
}