gen_input.c 3.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130
  1. #include <time.h>
  2. #include <stdlib.h>
  3. #include <stdio.h>
  4. #ifdef FP_NUMBER
  5. typedef double FP_NUMBER;
  6. #else
  7. typedef float FP_NUMBER;
  8. #endif
  9. #define GET_RAND_FP ((FP_NUMBER)rand()/((FP_NUMBER)(RAND_MAX)+(FP_NUMBER)(1)))
  10. char L_FNAME[32], U_FNAME[32], A_FNAME[32];
  11. int main (int argc, char **argv){
  12. int i,j,k,MatrixDim;
  13. FP_NUMBER sum, **L, **U, **A;
  14. FILE *fl,*fu,*fa;
  15. if ( argc < 2) {
  16. printf("./gen_input [Matrix_Dimension_size]\n");
  17. return 1;
  18. }
  19. MatrixDim = atoi(argv[1]);
  20. L = (FP_NUMBER **) malloc(sizeof(FP_NUMBER*)*MatrixDim);
  21. U = (FP_NUMBER **) malloc(sizeof(FP_NUMBER*)*MatrixDim);
  22. A = (FP_NUMBER **) malloc(sizeof(FP_NUMBER*)*MatrixDim);
  23. if ( !L || !U || !A){
  24. printf("Can not allocate memory\n");
  25. if (L) free(L);
  26. if (U) free(U);
  27. if (A) free(A);
  28. return 1;
  29. }
  30. srand(time(NULL));
  31. sprintf(L_FNAME, "l-%d.dat", MatrixDim);
  32. fl = fopen(L_FNAME, "wb");
  33. if (fl == NULL) {
  34. printf("Cannot open file %s\n", L_FNAME);
  35. return 1;
  36. }
  37. sprintf(U_FNAME, "u-%d.dat", MatrixDim);
  38. fu = fopen(U_FNAME, "wb");
  39. if (fu == NULL) {
  40. printf("Cannot open file %s\n", U_FNAME);
  41. return 1;
  42. }
  43. sprintf(A_FNAME, "%d.dat", MatrixDim);
  44. fa = fopen(A_FNAME, "wb");
  45. if (!fa) {
  46. printf("Cannot open file %s\n", A_FNAME);
  47. return 1;
  48. }
  49. for (i=0; i < MatrixDim; i ++){
  50. L[i]=(FP_NUMBER*)malloc(sizeof(FP_NUMBER)*MatrixDim);
  51. U[i]=(FP_NUMBER*)malloc(sizeof(FP_NUMBER)*MatrixDim);
  52. A[i]=(FP_NUMBER*)malloc(sizeof(FP_NUMBER)*MatrixDim);
  53. }
  54. #if 1
  55. #pragma omp parallel for default(none)\
  56. private(i,j) shared(L,U,MatrixDim)
  57. #endif
  58. for (i=0; i < MatrixDim; i ++){
  59. for (j=0; j < MatrixDim; j++){
  60. if ( i == j) {
  61. L[i][j] = 1.0;
  62. U[i][j] = GET_RAND_FP;
  63. } else if (i < j){
  64. L[i][j] = 0;
  65. U[i][j] = GET_RAND_FP;
  66. } else { // i > j
  67. L[i][j] = GET_RAND_FP;
  68. U[i][j] = 0;
  69. }
  70. }
  71. }
  72. #if 1
  73. #pragma omp parallel for default(none) \
  74. private(i,j,k,sum) shared(L,U,A,MatrixDim)
  75. #endif
  76. for (i=0; i < MatrixDim; i++ ) {
  77. for (j=0; j < MatrixDim; j++){
  78. sum = 0;
  79. for(k=0; k < MatrixDim; k++)
  80. sum += L[i][k]*U[k][j];
  81. A[i][j] = sum;
  82. }
  83. }
  84. for (i=0; i < MatrixDim; i ++) {
  85. for (j=0; j < MatrixDim; j++)
  86. fprintf(fl, "%f ", L[i][j]);
  87. fprintf(fl, "\n");
  88. }
  89. fclose(fl);
  90. for (i=0; i < MatrixDim; i ++) {
  91. for (j=0; j < MatrixDim; j++)
  92. fprintf(fu, "%f ", U[i][j]);
  93. fprintf(fu, "\n");
  94. }
  95. fclose(fu);
  96. fprintf(fa, "%d\n", MatrixDim);
  97. for (i=0; i < MatrixDim; i ++) {
  98. for (j=0; j < MatrixDim; j++)
  99. fprintf(fa, "%f ", A[i][j]);
  100. fprintf(fa, "\n");
  101. }
  102. fclose(fa);
  103. for (i = 0; i < MatrixDim; i ++ ){
  104. free(L[i]);
  105. free(U[i]);
  106. free(A[i]);
  107. }
  108. free(L);
  109. free(U);
  110. free(A);
  111. return 0;
  112. }