3

私は、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);' です。サイズが約30cm以上の場合 7000。

4

2 に答える 2

2

実行時に固定サイズ (または少なくとも固定最大サイズ) であり、物理 RAM よりも大きい場合は、メモリ マップ ファイルを使用することもできます。アクセスは少なくとも RAM+swap と同じくらい速く、ディスクにシリアル化されたデータを無料で取得できます。アドレス空間より全体的に大きいマップされたファイルの領域 (つまり、ウィンドウ) のビューをマップすることもできます。

一部の領域で高詳細が必要であるが、均一ではないために多数のセルが必要な場合は、octree を検討できます。次に、いくつかの部分で段階的に細かい解像度を保存できます。また、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 に答える