matrixcpu.c - 行列で遊んでみた

matrixcpu.c

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

// If you want OMP Power
#define OMP_ENABLE
#ifdef OMP_ENABLE
#include <omp.h>
#endif

/*n正方行列のサイズを定義*/
#define MATRIX_SIZE 1024

int** alc(){
  int** rmat;
  rmat = malloc(MATRIX_SIZE * sizeof(int*));

  #pragma omp parallel for
  for(int i = 0 ; i < MATRIX_SIZE; i++){
    rmat[i] = malloc(MATRIX_SIZE * sizeof(int));
    #pragma omp parallel for
    for(int j = 0; j < MATRIX_SIZE; j++)
      rmat[i][j] = 0;
  }

  return rmat;
}

int** conv2S(int* smat, int n){
  int** rmat;
  rmat = malloc(n * sizeof(int*));

  #pragma omp parallel for
  for(int i = 0 ; i < n; i++){
    rmat[i] = malloc(n * sizeof(int));
    #pragma omp parallel for
    for(int j = 0; j < n; j++)
      rmat[i][j] = smat[i * j];
  }

  return rmat;
}

void fray(int** ary){
  for(int i = 0; i < MATRIX_SIZE; i++)
      free(ary[i]);
  free(ary);
}

void print(int** mat){
  #pragma omp parallel for
  for(int i = 0; i < MATRIX_SIZE; i++){
    #pragma omp parallel for
    for(int j = 0; j < MATRIX_SIZE; j++)
      printf("%d ", mat[i][j]);
    printf("\n");
  }
}

int main(int argc, char** argv){
  unsigned int col_idx, row_idx, scan_idx;
  int* matA;
  int* matB;
  int* matC;
  int** matD;
  int** matE;

  /*タイマー作成*/
  time_t Start, Stop;

  /*int型のn×n領域をメモリに確保*/
  matA = (int*)malloc(sizeof(int) * MATRIX_SIZE * MATRIX_SIZE);
  matB = (int*)malloc(sizeof(int) * MATRIX_SIZE * MATRIX_SIZE);
  matC = (int*)malloc(sizeof(int) * MATRIX_SIZE * MATRIX_SIZE);

  #pragma omp parallel for
  for (col_idx = 0; col_idx < MATRIX_SIZE; col_idx++)
    #pragma omp parallel for
    for (row_idx = 0; row_idx < MATRIX_SIZE; row_idx++) {
      matA[col_idx * MATRIX_SIZE + row_idx] = rand() % (MATRIX_SIZE * MATRIX_SIZE);
      matB[col_idx * MATRIX_SIZE + row_idx] = rand() % (MATRIX_SIZE * MATRIX_SIZE);
      matC[col_idx * MATRIX_SIZE + row_idx] = 0;
    }

  time(&Start);

  #pragma omp parallel for
  for (col_idx = 0; col_idx < MATRIX_SIZE; col_idx++)
    #pragma omp parallel for
    for (row_idx = 0; row_idx < MATRIX_SIZE; row_idx++)
      #pragma omp parallel for
      for (scan_idx = 0; scan_idx < MATRIX_SIZE; scan_idx++)
        matC[col_idx * MATRIX_SIZE + row_idx] += matA[col_idx * MATRIX_SIZE + scan_idx] * matB[scan_idx * MATRIX_SIZE + row_idx];

  matD = conv2S(matC, MATRIX_SIZE);
  //print(matD);
  //printf("-------\n");
  matE = alc();

  #pragma omp parallel for
  for(int i = 0; i < MATRIX_SIZE; i++)
    #pragma omp parallel for
    for(int j = 0; j < MATRIX_SIZE; j++)
      matE[i][j] = matD[MATRIX_SIZE - j - 1][i];

  //print(matE);

  time(&Stop);
  printf("Processing time: %ld (sec)\n", Stop - Start);

  /*メモリを解放*/
  free(matA);
  free(matB);
  free(matC);
  fray(matD);
  fray(matE);
  return 0;
}