fork download 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. #define N 10
  6. #define TMAX 100
  7. #define EPS (10e-6)
  8.  
  9. void init_vector();
  10. double cg_method(double **matrix_a, double *vector_x, double *vector_b);
  11. void init_vector(double *vector_x);
  12. void vector_x_matrix(double *vector_y, double **matrix_a,double *vector_x);
  13. double dot_product(double *vector_x, double *vector_y);
  14. double vector_norm(double *vector_x);
  15.  
  16. int main (void){
  17.     FILE *fin, *fout;
  18.     static double **matrix_a;
  19.     static double *vector_b;
  20.     static double *vector_x;
  21.     int i,j;
  22.  
  23.     matrix_a = (double**)malloc(sizeof(double*)*N);
  24.     for(i=0;i<N;i++){
  25.         matrix_a[i] = (double*)malloc(sizeof(double)*N);
  26.     }
  27.  
  28.     vector_b = (double*)malloc(sizeof(double)*N);
  29.     vector_x = (double*)malloc(sizeof(double)*N);
  30.  
  31.     if((fin = fopen("input_sp.dat","r"))==NULL)
  32.     {
  33.         printf("no file :input_sp.dat\n");
  34.     }       
  35.     if((fout = fopen("output_sp.dat","w"))==NULL)
  36.     {
  37.         printf("do not maike file : output_sp.dat\n");
  38.         exit(1);
  39.     }
  40.  
  41.     fprintf( fout, "A\n");
  42.     printf("matrix_a\n");
  43.     for( i = 0 ; i < N ; i++)
  44.     {
  45.         for( j = 0 ; j < N ; j++)
  46.         {
  47.             fscanf(fin, "%lf", &matrix_a[i][j]);
  48.             fprintf(fout, "%5.2f\t", matrix_a[i][j]);
  49.             printf("%lf\n",matrix_a[i][j]);
  50.         }
  51.         printf("\n");
  52.         fprintf( fout, "\n");
  53.     }
  54.  
  55.     fprintf( fout, "b\n");
  56.     printf("vector_b\n");
  57.     for( i = 0 ; i < N ; i++)
  58.     {
  59.         fscanf(fin, "%lf", &vector_b[i]);
  60.         fprintf(fout, "%5.2f\t", vector_b[i]);
  61.         fprintf( fout, "\n");
  62.         printf("%lf\n",vector_b[i]);
  63.     }
  64.  
  65.     printf("\n");
  66.  
  67.     fprintf( fout, "x\n");
  68.     printf("vector_x\n");
  69.     for( i = 0 ; i < N ; i++)
  70.     {
  71.         fscanf(fin, "%lf", &vector_x[i]);
  72.         fprintf(fout, "%5.2f\t", vector_x[i]);
  73.         fprintf( fout, "\n");
  74.         printf("%lf\n",vector_x[i]);
  75.     }
  76.  
  77.     printf("\n");
  78.  
  79.     cg_method(matrix_a, vector_x, vector_b);
  80.  
  81.     fprintf( fout, "Ax=b is solved\n");
  82.     for( i = 1 ; i <= N ; i++ )
  83.     {
  84.         fprintf(fout, "%f\n", vector_x[i]);
  85.     }
  86.     fclose(fin); 
  87.     fclose(fout);
  88.  
  89.  
  90.     return 0;
  91. }
  92.  
  93. void init_vector(double *vector_x)
  94. {
  95.     int i;
  96.     for(i = 0; i < N; i++){
  97.         vector_x[i] = 0;
  98.     }
  99. }
  100.  
  101. void vector_x_matrix(double *vector_y,double **matrix_a,double *vector_x)
  102. {
  103.     int    i, j;
  104.     double vxm;
  105.     for(i = 0; i < N; i++){
  106.         vxm = 0;
  107.         for(j = 0; j < N; j++){
  108.             vxm += matrix_a[i][j] * vector_x[i];
  109.         }
  110.         vector_y[i] = vxm;
  111.     }
  112. }
  113.  
  114. double dot_product(double *vector_x,double *vector_y)
  115. {
  116.     int    i;
  117.     double dot_p = 0;
  118.     for(i = 0; i < N; i++){
  119.         dot_p += vector_x[i] * vector_y[i];
  120.     }
  121.     return dot_p;
  122. }
  123.  
  124. double vector_norm(double *vector_x)
  125. {
  126.     int i;
  127.     double norm = 0;
  128.     for(i = 0; i < N; i++){
  129.         norm += fabs(vector_x[i]);
  130.     }
  131.     return norm;
  132. }
  133.  
  134. double cg_method(double **matrix_a, double *vector_x, double *vector_b)
  135. {
  136.     double *vector_p;
  137.     double *vector_r;
  138.     double *vector_ax;
  139.     double *vector_ap;
  140.     int i, iter;
  141.  
  142.     vector_p = (double*)malloc(sizeof(double)*N); 
  143.     vector_r = (double*)malloc(sizeof(double)*N);
  144.     vector_ax = (double*)malloc(sizeof(double)*N);
  145.     vector_ap = (double*)malloc(sizeof(double)*N);
  146.  
  147.     vector_x_matrix(vector_ax, matrix_a, vector_x);
  148.  
  149.     for(i = 0; i < N; i++){
  150.         vector_p[i] = vector_b[i] - vector_ax[i];
  151.         vector_r[i] = vector_p[i];
  152.     }
  153.  
  154.     for(iter = 1; iter < TMAX; iter++){
  155.         double alpha, beta, err = 0;
  156.  
  157.         vector_x_matrix(vector_ap, matrix_a, vector_p);
  158.         alpha = dot_product(vector_p, vector_r)/dot_product(vector_p, vector_ap);
  159.  
  160.         for(i = 0; i < N; i++){
  161.             vector_x[i] += +alpha * vector_p[i];
  162.             vector_r[i] += -alpha * vector_ap[i];
  163.         }
  164.  
  165.         err = vector_norm(vector_r);
  166.         printf("LOOP : %d\t Error : %g\n", iter, err);
  167.         if(EPS > err) break;
  168.  
  169.         beta = -dot_product(vector_r, vector_ap)/dot_product(vector_p, vector_ap);
  170.         for(i = 0; i < N; i++){
  171.             vector_p[i] = vector_r[i] + beta * vector_p[i];
  172.         }
  173.     }
  174. }
Runtime error #stdin #stdout 0.01s 1852KB
comments (?)
stdin 
Standard input is empty
stdout 
no file :input_sp.dat
do not maike file : output_sp.dat
https://ideone.com/OWamE
language:
C (gcc 8.3)
created:
12 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!