Matrices/matrices.c at master · joanmawo/Matrices · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

int main(int argc, char **argv){


  float *datos;
  gsl_vector *posYgsl;
  gsl_matrix *Original;
  gsl_matrix *Ggsl;
  gsl_matrix *GTgsl;
  gsl_matrix *GtimeGTgsl;
  gsl_matrix *inverse;
  gsl_vector *Vgsl;

  int lines = 0;
  int i;
  float temp1, temp2, temp3;

  FILE *file;
  file = fopen(argv[1], "r");

  int n_lines = 0;
  int columnas = 2;
  int c=0;

  do{
    c = fgetc(file);
    if(c=='\n'){
      n_lines++;
    }
  }while(c!=EOF);

  rewind(file);

  Original = gsl_matrix_calloc (n_lines, 2);
  gsl_matrix_fscanf (file, Original);
  lines = n_lines;

  fprintf(stdout, "n_lines %d\n", n_lines);

  posYgsl = gsl_vector_calloc(lines);
  Ggsl = gsl_matrix_alloc(n_lines, 3);

  for(i = 0; i<n_lines; i++){
    gsl_vector_set (posYgsl, i, gsl_matrix_get(Original, i,0));
  }

  for(i = 0; i<n_lines; i++){

    temp1 = 1;
    temp2 = gsl_matrix_get(Original,i, 0);
    temp3 = (1/2)*gsl_matrix_get(Original,i, 0)*gsl_matrix_get(Original,i, 0);

    gsl_matrix_set(Ggsl, i, 0, temp1);
    gsl_matrix_set(Ggsl, i, 1, temp2);
    gsl_matrix_set(Ggsl, i, 2, temp3);
  }


  GTgsl = gsl_matrix_alloc(3, lines);

  gsl_matrix_transpose_memcpy(GTgsl,Ggsl);

  GtimeGTgsl = gsl_matrix_alloc(3, 3);
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, GTgsl, Ggsl, 0.0, GtimeGTgsl);

  fprintf(stdout, "n_lines %d\n", n_lines);

  gsl_matrix *LU;
  gsl_permutation *p=gsl_permutation_alloc (3);
  int signum;

  gsl_linalg_LU_decomp(GtimeGTgsl, p, &signum);
  gsl_linalg_LU_invert(GtimeGTgsl, p, inverse);

  Vgsl = gsl_vector_alloc (lines);
  gsl_matrix *temp;
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inverse, GTgsl, 0.0, temp);

  gsl_blas_dgemv (CblasNoTrans, 1.0, temp, posYgsl, 0.0, Vgsl);

  FILE *output;
  output = fopen("parametros_movimiento.dat", "w");

  fprintf(output, "%f %f %f" ,gsl_vector_get (Vgsl, 0), gsl_vector_get (Vgsl, 1) ,gsl_vector_get (Vgsl, 2));


  return 1;

}