Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <stdio.h>
- #include <stdlib.h>
- #define DEBUG 0
- void swap(double *s, int M, int N, int i, int l) {
- double t;
- for(int j = 0; j < N; ++j) {
- t = s[i*N+j];
- s[i*N+j] = s[l*N+j];
- s[l*N+j] = t;
- }
- }
- int non0pivot(double *s, int M, int N, int i, int j) {
- for(int i = 0; i < N; ++i)
- if(s[i*N+j] != 0)
- return i;
- return -1;
- }
- void readS(FILE *in, double **s, int *M, int *N) {
- int x;
- x = fscanf(in,"%d",M);
- *N = *M+1;
- #if DEBUG
- printf("%d %d\n",*M,*N);
- #endif
- *s = calloc((*M)*(*N),sizeof(double));
- for(int i = 0; i < *M**N; ++i) {
- #if DEBUG
- printf("i: %d\n",i);
- #endif
- x = fscanf(in,"%lf",&((*s)[i]));
- }
- }
- void printS(FILE *out, double *s, int M, int N) {
- for(int i = 0; i < M; ++i) {
- for(int j = 0; j < N; ++j) {
- if(s[i*N+j] >= 0) {
- printf(" ");
- fprintf(out," ");
- }
- printf("%.3lf ",s[i*N+j]);
- fprintf(out,"%.3lf ",s[i*N+j]);
- }
- putchar('\n');
- }
- putchar('\n');
- }
- void mult(double *s, double *x, double *c, int M, int N) {
- for(int i = 0; i < M; ++i)
- for(int j = 0; j < N-1; ++j)
- for(int k = 0; k < M; ++k)
- c[k] += s[i*N+k]*x[k];
- }
- int main(int argc, char **argv) {
- if(argc != 3)
- return -1;
- FILE *in = fopen(argv[1],"r");
- FILE *out = fopen(argv[2],"w");
- int M, N;
- double *s;
- readS(in,&s,&M,&N);
- printS(out,s,M,N);
- double *x = calloc(M,sizeof(double));
- double *b = calloc(M,sizeof(double));
- double *c = calloc(M,sizeof(double));
- #if DEBUG
- printf("linha: %d\n",__LINE__);
- #endif
- for(int i = 0, j = 0; i < M; ++i, ++j) {
- #if DEBUG
- printf("linha: %d, i = %d\n",__LINE__,i);
- #endif
- if(s[i*N+j] == 0) {
- #if DEBUG
- printf("linha: %d, i = %d\n",__LINE__,i);
- #endif
- int l = non0pivot(s,M,N,i+1,j);
- if(l == -1) {
- ++j;
- break;
- }
- swap(s,M,N,i,l);
- }
- #if DEBUG
- printf("linha: %d, i = %d\n",__LINE__,i);
- #endif
- for(int k = 0; k < M; ++k) {
- if(i == k)
- continue;
- register double x = s[k*N+j]/s[i*N+j];
- s[k*N+i] = 0.;
- for(int h = i+1; h < N; ++h)
- s[k*N+h] -= x*s[i*N+h];
- }
- printS(out,s,M,N);
- }
- for(int i = 0, j = 0; i < M; ++i, ++j) {
- if(s[i*N+j] == 0) {
- ++j;
- continue;
- }
- register double x = s[i*N+j];
- for(int k = j; k < N; ++k)
- s[i*N+k] = s[i*N+k]/x;
- }
- printS(out,s,M,N);
- for(int i = 0, j = 0; i < M; ++i, ++j) {
- if(s[i*N+j] == 0) {
- ++j;
- continue;
- }
- x[i] = s[i*N+j];
- b[i] = s[i*N+N-1];
- }
- mult(s,x,c,M,N);
- for(int i = 0; i < M; ++i)
- printf("err(%d): %lf\n", i, b[i]-c[i]);
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
