HPC/並列プログラミングポータルでは、HPC(High Performance Computing)プログラミングや並列プログラミングに関する情報を集積・発信しています。

新着トピックス

「最適化・並列化コードを生み出す最新コンパイラ「インテル コンパイラー」サンプルコード

このページでは『最適化・並列化コードを生み出す最新コンパイラ「インテル コンパイラー」』記事内で使用したソースコードをまとめています。

リスト1 「コンパイラを変えるだけでパフォーマンス向上、インテル コンパイラの実力を見る」記事内で使用した原始的な行列演算プログラム
/*
 * Matrix calculation test program
 * matrix1.c
 */

#include stdio.h
#include stdlib.h
#include time.h
#include string.h

#define N 400
#define REPEAT 100
#define AT(matrix,i,j) (matrix[i*dim+j])

static unsigned long next_rand = 1;
void print_matrix( double* matrix, int dim );
double my_rand(void);
void prod_matrix( double* matrixA, double* matrixB, double* matrixC, int dim );
void cp_matrix( double* matrixA, double* matrixB, int dim );
void mult_matrix( double* matrixA, double* matrixB, double x, int dim );

int main( int argc, const char **argv) {
  int i, j;
  int rep;
  int dim = N;
  double* matrixA;
  double* matrixB;
  double* matrixC;
  clock_t clock_start, clock_end;

  clock_start = clock();

  /* allocate matrix */
  matrixA = (double*)malloc( sizeof(double) * dim * dim );
  matrixB = (double*)malloc( sizeof(double) * dim * dim );
  matrixC = (double*)malloc( sizeof(double) * dim * dim );

  if( (matrixA == NULL) ||
      (matrixB == NULL) ||
      (matrixC == NULL) ) {
    printf( "cannot allocate memory. exit.\n" );
    free(matrixA);
    free(matrixB);
    free(matrixC);
    return -1;
  }

  /* initialize matrix */
  for( i = 0; i  dim; i++ ) {
    for( j = 0; j  dim; j++ ) {
      AT(matrixA,i,j) = my_rand();
      AT(matrixB,i,j) = my_rand();
    }
  }

  for( rep = 0; rep  REPEAT; rep++ ) {
    prod_matrix( matrixC, matrixA, matrixB, dim );
    mult_matrix( matrixA, matrixC, 2.0/(double)dim, dim );
  }


  /*  print_matrix(matrixC, dim);*/
  free(matrixA);
  free(matrixB);
  free(matrixC);

  clock_end = clock();

  printf( "start: %d end: %d elapse: %d.\n", clock_start, clock_end, clock_end - clock_start );
  printf( "elapse time: %lf sec.", ((double)(clock_end - clock_start)) / (double)CLOCKS_PER_SEC );

  return 0;
}


/* A = B; copy matrix B to matrix A */
void cp_matrix( double* matrixA, double* matrixB,int dim ) {
  memcpy( matrixA, matrixB, sizeof(double)*dim*dim );
}

/* A = B * C; matrixA = matrixB * matrixC */
void prod_matrix( double* matrixA, double* matrixB, double* matrixC, int dim ) {
  int i, j, m;

  memset( matrixA, 0, sizeof(double)*dim*dim );

  for( i = 0; i  dim; i++ ) {
    for( m = 0; m  dim; m++ ) {
      for( j = 0; j  dim; j++ ) {
        AT(matrixA,i,j) += AT(matrixB,i,m) * AT(matrixC, m, j);
      }
    }
  }
}

/* matrixA = matrixB * x */
void mult_matrix( double* matrixA, double* matrixB, double x, int dim ) {
  int i, j;

  for( i = 0; i  dim; i++ ) {
    for( j = 0; j  dim; j++ ) {
      AT(matrixA,i,j) = AT(matrixB,i,j) * x;
    }
  }
}

void print_matrix( double* matrix, int dim ) {
  int i,j;

  for( i = 0; i  dim; i++ ) {
    for( j = 0; j  dim; j++ ) {
      // output
      printf( "%lf ", AT(matrix,i,j) );
    }
    printf( "\n" );
  }

}

double my_rand(void) {
    next_rand = next_rand * 1103515245 + 12345;
    return (((double)((unsigned)(next_rand/65536) % 32768))/32768);
}