añadiendo benchmark rodinia

codes/lud/Makefile

+all: lud_cuda 
+
+lud_cuda:
+	cd cuda; make
+
+clean:
+	cd cuda; make clean

codes/lud/base/Makefile

+# Makefile for LU Decomposition
+#
+#
+
+# DEBUG can be set to YES to include debugging info, or NO otherwise
+DEBUG          := YES
+
+# PROFILE can be set to YES to include profiling info, or NO otherwise
+PROFILE        := NO
+
+# ------------  name of the executable  ----------------------------------------
+EXECUTABLE      := lud_base
+
+# ------------  list of all source files  --------------------------------------
+SOURCES         := lud.c lud_base.c ../common/common.c 
+
+# ------------  compiler  ------------------------------------------------------
+CC              := gcc
+CXX             := g++
+
+# ------------  compiler flags  ------------------------------------------------
+DEBUG_CFLAGS    := -Wall -O0 -g
+RELEASE_CFLAGS  := -Wall -O3 
+COMMON_CFLAGS   := 
+
+# ------------  compiler defs  ------------------------------------------------
+DEBUG_CDEFS    := -D_DEBUG_RESULT
+RELEASE_CDEFS  := 
+COMMON_CDEFS   := 
+
+# ------------  linker flags  --------------------------------------------------
+DEBUG_LDFLAGS    := -g
+RELEASE_LDFLAGS  := 
+COMMON_LDFLAGS   := 
+
+ifeq (YES, ${DEBUG})
+  CFLAGS       := ${DEBUG_CFLAGS} ${DEBUG_CDEFS}
+  CXXFLAGS     := ${DEBUG_CXXFLAGS} ${DEBUG_CDEFS}
+  LDFLAGS      := ${DEBUG_LDFLAGS}
+else
+  CFLAGS       := ${RELEASE_CFLAGS} ${RELEASE_CDEFS}
+  CXXFLAGS     := ${RELEASE_CXXFLAGS} ${RELEASE_CDEFS}
+  LDFLAGS      := ${RELEASE_LDFLAGS}
+endif
+
+CFLAGS         := ${CFLAGS} ${COMMON_CFLAGS}
+CXXFLAGS       := ${CXXFLAGS} ${COMMON_CFLAGS}
+LDFLAGS        := ${LDFLAGS} ${COMMON_LDFLAGS}
+
+ifeq (YES, ${PROFILE})
+  CFLAGS       := ${CFLAGS}   -pg -O3
+  CXXFLAGS     := ${CXXFLAGS} -pg -O3
+  LDFLAGS      := ${LDFLAGS}  -pg
+endif
+
+# ------------  additional system include directories  -------------------------
+GLOBAL_INC_DIR  =
+
+# ------------  private include directories  -----------------------------------
+LOCAL_INC_DIR   = ../common
+
+# ------------  system libraries  (e.g. -lm )  ---------------------------------
+SYS_LIBS        = 
+
+# ------------  additional system library directories  -------------------------
+GLOBAL_LIB_DIR  =
+
+# ------------  additional system libraries  -----------------------------------
+GLOBAL_LIBS     =
+
+# ------------  private library directories  -----------------------------------
+LOCAL_LIB_DIR   = 
+
+# ------------  private libraries  (e.g. libxyz.a )  ---------------------------
+LOCAL_LIBS      =
+
+
+C_SOURCES       = $(filter     %.c, $(SOURCES))
+CPP_SOURCES     = $(filter-out %.c, $(SOURCES))
+ALL_INC_DIR     = $(addprefix -I, $(LOCAL_INC_DIR) $(GLOBAL_INC_DIR))
+ALL_LIB_DIR     = $(addprefix -L, $(LOCAL_LIB_DIR) $(GLOBAL_LIB_DIR))
+GLOBAL_LIBSS    = $(addprefix $(GLOBAL_LIB_DIR)/, $(GLOBAL_LIBS))
+LOCAL_LIBSS     = $(addprefix $(LOCAL_LIB_DIR)/, $(LOCAL_LIBS))
+ALL_CFLAGS      = $(CFLAGS) $(ALL_INC_DIR)
+ALL_LFLAGS      = $(LDFLAGS) $(ALL_LIB_DIR)
+BASENAMES       = $(basename $(SOURCES))
+
+# ------------  generate the names of the object files  ------------------------
+OBJECTS         = $(addsuffix .o,$(BASENAMES))
+ 
+# ------------  make the executable (the default goal)  ------------------------
+$(EXECUTABLE):	$(OBJECTS)
+ifeq ($(strip $(CPP_SOURCES)),)
+	$(CC)  $(ALL_LFLAGS) -o $(EXECUTABLE) $(OBJECTS) $(LOCAL_LIBSS) $(GLOBAL_LIBSS) $(SYS_LIBS)
+else
+	$(CXX) $(ALL_LFLAGS) -o $(EXECUTABLE) $(OBJECTS) $(LOCAL_LIBSS) $(GLOBAL_LIBSS) $(SYS_LIBS)
+endif
+
+# ------------  make the objects  ----------------------------------------------
+%.o:		%.c
+				$(CC)  -c $(ALL_CFLAGS) $< -o $@
+
+%.o:		%.cc
+				$(CXX) -c $(ALL_CFLAGS) $< -o $@
+
+%.o:		%.cpp
+				$(CXX) -c $(ALL_CFLAGS) $< -o $@
+
+%.o:		%.C
+				$(CXX) -c $(ALL_CFLAGS) $< -o $@
+
+clean:
+	-rm --force $(EXECUTABLE) $(OBJECTS) *~
+
+.PHONY: clean
+
+# ==============================================================================
+# vim: set tabstop=2: set shiftwidth=2:

codes/lud/base/lud.c

+/*
+ * =====================================================================================
+ *
+ *       Filename:  suite.c
+ *
+ *    Description:  The main wrapper for the suite
+ *
+ *        Version:  1.0
+ *        Created:  10/22/2009 08:40:34 PM
+ *       Revision:  none
+ *       Compiler:  gcc
+ *
+ *         Author:  Liang Wang (lw2aw), lw2aw@virginia.edu
+ *        Company:  CS@UVa
+ *
+ * =====================================================================================
+ */
+
+#include <stdio.h>
+#include <unistd.h>
+#include <getopt.h>
+#include <stdlib.h>
+#include <assert.h>
+
+#include "common.h"
+
+static int do_verify = 0;
+
+static struct option long_options[] = {
+      /* name, has_arg, flag, val */
+      {"input", 1, NULL, 'i'},
+      {"size", 1, NULL, 's'},
+      {"verify", 0, NULL, 'v'},
+      {0,0,0,0}
+};
+
+extern void
+lud_base(float *m, int matrix_dim);
+
+int
+main ( int argc, char *argv[] )
+{
+  int matrix_dim = 32; /* default matrix_dim */
+  int opt, option_index=0;
+  func_ret_t ret;
+  const char *input_file = NULL;
+  float *m, *mm;
+  stopwatch sw;
+
+  while ((opt = getopt_long(argc, argv, "::vs:i:", 
+                            long_options, &option_index)) != -1 ) {
+      switch(opt){
+        case 'i':
+          input_file = optarg;
+          break;
+        case 'v':
+          do_verify = 1;
+          break;
+        case 's':
+          matrix_dim = atoi(optarg);
+          fprintf(stderr, "Currently not supported, use -i instead\n");
+          fprintf(stderr, "Usage: %s [-v] [-s matrix_size|-i input_file]\n", argv[0]);
+          exit(EXIT_FAILURE);
+        case '?':
+          fprintf(stderr, "invalid option\n");
+          break;
+        case ':':
+          fprintf(stderr, "missing argument\n");
+          break;
+        default:
+          fprintf(stderr, "Usage: %s [-v] [-s matrix_size|-i input_file]\n",
+                  argv[0]);
+          exit(EXIT_FAILURE);
+      }
+  }
+
+  if ( (optind < argc) || (optind == 1)) {
+      fprintf(stderr, "Usage: %s [-v] [-s matrix_size|-i input_file]\n", argv[0]);
+      exit(EXIT_FAILURE);
+  }
+  
+  if (input_file) {
+      printf("Reading matrix from file %s\n", input_file);
+      ret = create_matrix_from_file(&m, input_file, &matrix_dim);
+      if (ret != RET_SUCCESS) {
+          m = NULL;
+          fprintf(stderr, "error create matrix from file %s\n", input_file);
+          exit(EXIT_FAILURE);
+      }
+  } else {
+    printf("No input file specified!\n");
+    exit(EXIT_FAILURE);
+  } 
+
+  if (do_verify){
+    printf("Before LUD\n");
+    print_matrix(m, matrix_dim);
+    matrix_duplicate(m, &mm, matrix_dim);
+  }
+
+  stopwatch_start(&sw);
+  lud_base(m, matrix_dim);
+  stopwatch_stop(&sw);
+  printf("Time consumed(ms): %lf\n", 1000*get_interval_by_sec(&sw));
+
+  if (do_verify){
+    printf("After LUD\n");
+    print_matrix(m, matrix_dim);
+    printf(">>>Verify<<<<\n");
+    lud_verify(mm, m, matrix_dim); 
+    free(mm);
+  }
+
+  free(m);
+
+  return EXIT_SUCCESS;
+}				/* ----------  end of function main  ---------- */

codes/lud/base/lud_base.c

+void lud_base(float *a, int size)
+{
+     int i,j,k;
+     float sum;
+
+     for (i=0; i<size; i++){
+         for (j=i; j<size; j++){
+             sum=a[i*size+j];
+             for (k=0; k<i; k++) sum -= a[i*size+k]*a[k*size+j];
+             a[i*size+j]=sum;
+         }
+
+         for (j=i+1;j<size; j++){
+             sum=a[j*size+i];
+             for (k=0; k<i; k++) sum -=a[j*size+k]*a[k*size+i];
+             a[j*size+i]=sum/a[i*size+i];
+         }
+     }
+}

codes/lud/common/common.c

+#include <string.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <time.h>
+#include <math.h>
+
+#include "common.h"
+
+void stopwatch_start(stopwatch *sw){
+    if (sw == NULL)
+        return;
+
+    bzero(&sw->begin, sizeof(struct timeval));
+    bzero(&sw->end  , sizeof(struct timeval));
+
+    gettimeofday(&sw->begin, NULL);
+}
+
+void stopwatch_stop(stopwatch *sw){
+    if (sw == NULL)
+        return;
+
+    gettimeofday(&sw->end, NULL);
+}
+
+double 
+get_interval_by_sec(stopwatch *sw){
+    if (sw == NULL)
+        return 0;
+    return ((double)(sw->end.tv_sec-sw->begin.tv_sec)+(double)(sw->end.tv_usec-sw->begin.tv_usec)/1000000);
+}
+
+int 
+get_interval_by_usec(stopwatch *sw){
+    if (sw == NULL)
+        return 0;
+    return ((sw->end.tv_sec-sw->begin.tv_sec)*1000000+(sw->end.tv_usec-sw->begin.tv_usec));
+}
+
+func_ret_t 
+create_matrix_from_file(float **mp, const char* filename, int *size_p){
+  int i, j, size;
+  float *m;
+  FILE *fp = NULL;
+
+  fp = fopen(filename, "rb");
+  if ( fp == NULL) {
+      return RET_FAILURE;
+  }
+
+  fscanf(fp, "%d\n", &size);
+
+  m = (float*) malloc(sizeof(float)*size*size);
+  if ( m == NULL) {
+      fclose(fp);
+      return RET_FAILURE;
+  }
+
+  for (i=0; i < size; i++) {
+      for (j=0; j < size; j++) {
+          fscanf(fp, "%f ", m+i*size+j);
+      }
+  }
+
+  fclose(fp);
+
+  *size_p = size;
+  *mp = m;
+
+  return RET_SUCCESS;
+}
+
+
+func_ret_t
+create_matrix_from_random(float **mp, int size){
+  float *l, *u, *m;
+  int i,j,k;
+
+  srand(time(NULL));
+
+  l = (float*)malloc(size*size*sizeof(float));
+  if ( l == NULL)
+    return RET_FAILURE;
+
+  u = (float*)malloc(size*size*sizeof(float));
+  if ( u == NULL) {
+      free(l);
+      return RET_FAILURE;
+  }
+
+  for (i = 0; i < size; i++) {
+      for (j=0; j < size; j++) {
+          if (i>j) {
+              l[i*size+j] = GET_RAND_FP;
+          } else if (i == j) {
+              l[i*size+j] = 1;
+          } else {
+              l[i*size+j] = 0;
+          }
+      }
+  }
+
+  for (j=0; j < size; j++) {
+      for (i=0; i < size; i++) {
+          if (i>j) {
+              u[j*size+i] = 0;
+          }else {
+              u[j*size+i] = GET_RAND_FP; 
+          }
+      }
+  }
+
+  for (i=0; i < size; i++) {
+      for (j=0; j < size; j++) {
+          for (k=0; k <= MIN(i,j); k++)
+            m[i*size+j] = l[i*size+k] * u[j*size+k];
+      }
+  }
+
+  free(l);
+  free(u);
+
+  *mp = m;
+
+  return RET_SUCCESS;
+}
+
+void
+matrix_multiply(float *inputa, float *inputb, float *output, int size){
+  int i, j, k;
+
+  for (i=0; i < size; i++)
+    for (k=0; k < size; k++)
+      for (j=0; j < size; j++)
+        output[i*size+j] = inputa[i*size+k] * inputb[k*size+j];
+
+}
+
+func_ret_t
+lud_verify(float *m, float *lu, int matrix_dim){
+  int i,j,k;
+  float *tmp = (float*)malloc(matrix_dim*matrix_dim*sizeof(float));
+  int errores=0;
+
+  for (i=0; i < matrix_dim; i ++)
+    for (j=0; j< matrix_dim; j++) {
+        float sum = 0;
+        float l,u;
+        for (k=0; k <= MIN(i,j); k++){
+            if ( i==k)
+              l=1;
+            else
+              l=lu[i*matrix_dim+k];
+            u=lu[k*matrix_dim+j];
+            sum+=l*u;
+        }
+        tmp[i*matrix_dim+j] = sum;
+    }
+  /* printf(">>>>>LU<<<<<<<\n"); */
+  /* for (i=0; i<matrix_dim; i++){ */
+  /*   for (j=0; j<matrix_dim;j++){ */
+  /*       printf("%f ", lu[i*matrix_dim+j]); */
+  /*   } */
+  /*   printf("\n"); */
+  /* } */
+  /* printf(">>>>>result<<<<<<<\n"); */
+  /* for (i=0; i<matrix_dim; i++){ */
+  /*   for (j=0; j<matrix_dim;j++){ */
+  /*       printf("%f ", tmp[i*matrix_dim+j]); */
+  /*   } */
+  /*   printf("\n"); */
+  /* } */
+  /* printf(">>>>>input<<<<<<<\n"); */
+  /* for (i=0; i<matrix_dim; i++){ */
+  /*   for (j=0; j<matrix_dim;j++){ */
+  /*       printf("%f ", m[i*matrix_dim+j]); */
+  /*   } */
+  /*   printf("\n"); */
+  /* } */
+  cmp_res(m,tmp,matrix_dim);
+  free(tmp);
+}
+void
+cmp_res (float *m, float *tmp, int matrix_dim){
+  int i,j;
+	int errores=0;
+  for (i=0; i<matrix_dim; i++){
+      for (j=0; j<matrix_dim; j++){
+          if ( fabs(m[i*matrix_dim+j]-tmp[i*matrix_dim+j]) > 0.0001)
+	  { errores++;
+	    //printf("dismatch at (%d, %d): (o)%f (n)%f\n", i, j, m[i*matrix_dim+j], tmp[i*matrix_dim+j]);
+	  }
+      }
+  }
+  if (errores) printf("\n Results = FAIL \n");
+  else
+	       printf("\n Results = PASS \n");
+}
+
+void
+matrix_duplicate(float *src, float **dst, int matrix_dim) {
+    int s = matrix_dim*matrix_dim*sizeof(float);
+   float *p = (float *) malloc (s);
+   memcpy(p, src, s);
+   *dst = p;
+}
+
+void
+print_matrix(float *m, int matrix_dim) {
+    int i, j;
+    for (i=0; i<matrix_dim;i++) {
+      for (j=0; j<matrix_dim;j++)
+        printf("%f ", m[i*matrix_dim+j]);
+      printf("\n");
+    }
+}
+
+
+// Generate well-conditioned matrix internally  by Ke Wang 2013/08/07 22:20:06
+
+func_ret_t
+create_matrix(float **mp, int size){
+  float *m;
+  int i,j;
+  float lamda = -0.001;
+  float coe[2*size-1];
+  float coe_i =0.0;
+
+  for (i=0; i < size; i++)
+    {
+      coe_i = 10*exp(lamda*i); 
+      j=size-1+i;     
+      coe[j]=coe_i;
+      j=size-1-i;     
+      coe[j]=coe_i;
+    }
+
+  m = (float*) malloc(sizeof(float)*size*size);
+  if ( m == NULL) {
+      return RET_FAILURE;
+  }
+
+  for (i=0; i < size; i++) {
+      for (j=0; j < size; j++) {
+	m[i*size+j]=coe[size-1-i+j];
+      }
+  }
+
+  *mp = m;
+
+  return RET_SUCCESS;
+}

codes/lud/common/common.h

+#ifndef _COMMON_H
+#define _COMMON_H
+
+#include <time.h>
+#include <sys/time.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+
+#define GET_RAND_FP ( (float)rand() /   \
+                     ((float)(RAND_MAX)+(float)(1)) )
+
+#define MIN(i,j) ((i)<(j) ? (i) : (j))
+
+typedef enum _FUNC_RETURN_CODE {
+    RET_SUCCESS,
+    RET_FAILURE
+}func_ret_t;
+
+typedef struct __stopwatch_t{
+    struct timeval begin;
+    struct timeval end;
+}stopwatch;
+
+void 
+stopwatch_start(stopwatch *sw);
+
+void 
+stopwatch_stop (stopwatch *sw);
+
+double 
+get_interval_by_sec(stopwatch *sw);
+
+int 
+get_interval_by_usec(stopwatch *sw);
+
+func_ret_t
+create_matrix_from_file(float **mp, const char *filename, int *size_p);
+
+func_ret_t
+create_matrix_from_random(float **mp, int size);
+
+func_ret_t
+create_matrix(float **mp, int size);
+
+func_ret_t
+lud_verify(float *m, float *lu, int size);
+
+void 
+cmp_res (float *m, float *tmp, int matrix_dim);
+void
+matrix_multiply(float *inputa, float *inputb, float *output, int size);
+
+void
+matrix_duplicate(float *src, float **dst, int matrix_dim);
+
+void
+print_matrix(float *mm, int matrix_dim);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

codes/lud/lud.conf

+[DEFAULT]
+
+debug =True 
+
+# Name of the gdb executable
+gdbExecName = /usr/local/cuda-10.1.243/bin/cuda-gdb 
+
+# Which fault model to use, 0 -> single; 1 -> double;
+# 2 -> random; 3 -> zeros; 4 -> least 16 significant bits (LSB);
+# 5 -> least 8 significant bits (LSB)
+# If you want multiple fault models, place them separated by ','
+# faultModel = 0,2,3
+faultModel = 0
+
+# Injection site
+# Can be:
+# RF -> Register File
+# INST_OUT -> Instruction Output (NOT IMPLEMENTED YET)
+# INST_composed -> Instruction Adress (NOT IMPLEMENTED YET)
+injectionSite = RF 
+
+# Max time factor to finish the app, this will be multiplied by the application running time
+# For example if your app spend 2s, and the maxWaitTimes is 5, the max running time before it is
+# Considered as a crash is 10s
+maxWaitTimes = 5
+
+# binary file of the application
+# Must be full path
+benchmarkBinary = /home/badia/rodinia_3.1/cuda/lud/cuda/lud_cuda
+
+#benchmarkBinary_noverificar = /home/badia/rodinia_3.1/cuda/lavaMD/lavaMD 
+
+# Commands to set the session inside GDB environment
+benchmarkArgs =  -s 10000 -v 
+
+benchmarkArgs_noverificar =  -s 10000
+
+# CSV output file. It will be overwrite at each injection
+csvFile = results/lud_RF.csv
+
+# You should create a script on the benchmark source folder to verify GOLD_OUTPUT x INJ_OUTPUT
+goldenCheckScript = ../rodinia_3.1/cuda/lud/sdc_check.sh
+
+# Number of signals that will be sent to the application
+seqSignals = 20
+# Initial sleep time in seconds before start sending signals
+# Generally the memory setup time
+initSleep = 0.9  
+
+kernels = lud_internal,lud_perimeter,lud_diagonal
+section_begin=lud.cu:151
+section_end=lud.cu:155
+
+#initSleep = 0.6 8192 

codes/lud/tools/Makefile

+gen_input: gen_input.c
+	gcc -o gen_input -O3 gen_input.c -fopenmp

codes/lud/tools/README

+gen_input will genearte the input matrix with given dimension. Running it in
+following way:
+
+./gen_input [matrix_dimension]
+
+e.g.:
+./gen_input 32
+will generate the input matrix of 32x32. The ouput contains three files:
+${matrix_dim}.dat, l-${matrix_dim}.dat, u-${matrix_dim}.dat. In the above
+example, the output should be 32.dat, l-32.dat, u-32.dat .

codes/lud/tools/gen_input.c

+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#ifdef FP_NUMBER
+typedef double FP_NUMBER;
+#else
+typedef float FP_NUMBER;
+#endif
+
+
+#define GET_RAND_FP ((FP_NUMBER)rand()/((FP_NUMBER)(RAND_MAX)+(FP_NUMBER)(1)))
+char L_FNAME[32], U_FNAME[32], A_FNAME[32];
+
+int main (int argc, char **argv){
+    int i,j,k,MatrixDim;
+    FP_NUMBER sum, **L, **U, **A;
+    FILE *fl,*fu,*fa;
+
+    if ( argc < 2) {
+        printf("./gen_input [Matrix_Dimension_size]\n");
+        return 1;
+    }
+
+    MatrixDim = atoi(argv[1]);
+    L = (FP_NUMBER **) malloc(sizeof(FP_NUMBER*)*MatrixDim);
+    U = (FP_NUMBER **) malloc(sizeof(FP_NUMBER*)*MatrixDim);
+    A = (FP_NUMBER **) malloc(sizeof(FP_NUMBER*)*MatrixDim);
+
+    if ( !L || !U || !A){
+        printf("Can not allocate memory\n");
+        if (L) free(L);
+        if (U) free(U);
+        if (A) free(A);
+        return 1;
+    }
+
+    srand(time(NULL));
+
+    sprintf(L_FNAME, "l-%d.dat", MatrixDim);
+    fl = fopen(L_FNAME, "wb");
+    if (fl == NULL) {
+        printf("Cannot open file %s\n", L_FNAME);
+        return 1;
+    }
+
+    sprintf(U_FNAME, "u-%d.dat", MatrixDim);
+    fu = fopen(U_FNAME, "wb");
+    if (fu == NULL) {
+        printf("Cannot open file %s\n", U_FNAME);
+        return 1;
+    }
+
+    sprintf(A_FNAME, "%d.dat", MatrixDim);
+    fa = fopen(A_FNAME, "wb");
+    if (!fa) {
+        printf("Cannot open file %s\n", A_FNAME);
+        return 1;
+    }
+
+    for (i=0; i < MatrixDim; i ++){
+        L[i]=(FP_NUMBER*)malloc(sizeof(FP_NUMBER)*MatrixDim);
+        U[i]=(FP_NUMBER*)malloc(sizeof(FP_NUMBER)*MatrixDim);
+        A[i]=(FP_NUMBER*)malloc(sizeof(FP_NUMBER)*MatrixDim);
+    }
+#if 1
+#pragma omp parallel for default(none)\
+    private(i,j) shared(L,U,MatrixDim)
+#endif
+    for (i=0; i < MatrixDim; i ++){
+        for (j=0; j < MatrixDim; j++){
+            if ( i == j) {
+                L[i][j] = 1.0;
+                U[i][j] = GET_RAND_FP;
+            } else if (i < j){
+                L[i][j] = 0;
+                U[i][j] = GET_RAND_FP;
+            } else { // i > j
+                L[i][j] = GET_RAND_FP;
+                U[i][j] = 0;
+            }
+        }
+    }
+
+#if 1
+#pragma omp parallel for default(none) \
+    private(i,j,k,sum) shared(L,U,A,MatrixDim)
+#endif
+    for (i=0; i < MatrixDim; i++ ) {
+        for (j=0; j < MatrixDim; j++){
+            sum = 0;
+            for(k=0; k < MatrixDim; k++)
+                sum += L[i][k]*U[k][j];
+            A[i][j] = sum;
+        }
+    }
+
+    for (i=0; i < MatrixDim; i ++) {
+        for (j=0; j < MatrixDim; j++)
+            fprintf(fl, "%f ", L[i][j]);
+        fprintf(fl, "\n");
+    }
+    fclose(fl);
+
+    for (i=0; i < MatrixDim; i ++) {
+        for (j=0; j < MatrixDim; j++)
+            fprintf(fu, "%f ", U[i][j]);
+        fprintf(fu, "\n");
+    }
+    fclose(fu);
+
+    fprintf(fa, "%d\n", MatrixDim);
+    for (i=0; i < MatrixDim; i ++) {
+        for (j=0; j < MatrixDim; j++)
+            fprintf(fa, "%f ", A[i][j]);
+        fprintf(fa, "\n");
+    }
+    fclose(fa);
+
+    for (i = 0; i < MatrixDim; i ++ ){
+        free(L[i]);
+        free(U[i]);
+        free(A[i]);
+    }
+    free(L);
+    free(U);
+    free(A);
+
+    return 0;
+}