3

我正在做一个需要创建 3D 数组、一些 2D 和 1D 数组的项目。3D 数组表示空间中的离散坐标，我需要很多点来解决我的问题。数组大小约为 2000*2000*2000。我需要在这些数组中存储“双”值。任何人都可以提出一个有效的方案来在 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。

4

2 回答 2

2

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

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

于 2012-08-02T13:57:48.153 回答
0

这是 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;
}
于 2012-08-02T15:02:34.193 回答