diff --git a/header.h b/header.h
index ce8bc5a66a5634c3897d9f6ae5ef13d402779762..df2af685171541080c1ea6f5cfad2952aa519748 100644
--- a/header.h
+++ b/header.h
@@ -3,6 +3,7 @@
#include <string.h>
#include <math.h>
#include <errno.h>
+#include <time.h>
#define ARGS_NUM 9
@@ -13,6 +14,8 @@ double Pipi = M_PI * M_PI; // Pipi = Pi * Pi
void print_errno(void);
void getParams(int argc, char* argv[], FILE *fp);
-double f(double x, double y);
-void u(int i, int j, double **u, double **f);
+//double f(double x, double y);
+//void u(int i, int j, double **u, double **f);
double residueNorm(double *r);
+double timestamp(void);
+//void sor(double *A, double *b, double *x, double w);
diff --git a/pdeSolver.c b/pdeSolver.c
index afd011be23a9e5616b07599227fd9bcf4630377f..c23fca311adc569defc0a6f491ace6b8bfe2a6b4 100644
--- a/pdeSolver.c
+++ b/pdeSolver.c
@@ -4,11 +4,38 @@
u(i,j) = f(x,y) + (u(i+1,j) + u(i-1,j))/Δx² + (u(i,j+1) + u(i,j-1))/Δy² + (-u(i+1,j)+u(i-1,j))/2Δx + (-u(i,j+1)+u(i,j-1))/2Δy
--------------------------------------------------------------------------------------------------------------------
2/Δx² + 2/Δy² + 4π²
-
+Resíduo:
+f(x,y) = (2/Δx² + 2/Δy² + 4π²) * u(i,j) - ((u(i+1,j) + u(i-1,j))/Δx² + (u(i,j+1) + u(i,j-1))/Δy² + (-u(i+1,j)+u(i-1,j))/2Δx + (-u(i,j+1)+u(i,j-1))/2Δy)
*/
//#define MAX_SIZE 100000000 // 100 MB.
//int inMemory = 0;
+int N = 9;
+
+void print_vector(double *x) {
+ int i;
+ for(i=0; i<N; ++i) {
+ printf("%f ", x[i]);
+ }
+ printf("\n");
+}
+
+void print_matrix(double *A) {
+ int i,j;
+ for(i=0; i<Nx * Ny; ++i) {
+ for(j=0; j<Nx * Ny; ++j) {
+ printf("%.0f ", A[in(i,j)]);
+ }
+ printf("\n");
+ }
+ printf("\n\n\n");
+}
+
+double timestamp(void) {
+ struct timeval tp;
+ gettimeofday(&tp, NULL);
+ return((double)(tp.tv_sec + tp.tv_usec/1000000.0));
+}
void print_errno(void) {
printf("%s",strerror(errno));
@@ -35,16 +62,22 @@ void getParams(int argc, char* argv[], FILE *fp) {
}
}
}
-
+/*
double f(double x, double y) { // If stored in memory, it will need Nx * Ny * 8 bytes of memory. Should we do it?
return 4 * M_PI * ( ( sin(2 * M_PI * x) ) * ( sinh(M_PI * y) ) + ( sin(2 * Pipi - M_PI * x) ) * ( sinh(Pipi - M_PI * y) ));
}
+*/
+double f(int n) { // If stored in memory, it will need Nx * Ny * 8 bytes of memory. Should we do it?
+ int i = n / Nx, j = n % Nx;
+ double x = i * Hx, y = j * Hy;
+ return 4 * M_PI * ( ( sin(2 * M_PI * x) ) * ( sinh(M_PI * y) ) + ( sin(2 * Pipi - M_PI * x) ) * ( sinh(Pipi - M_PI * y) ));
+}
inline int in(int i, int j) {
// Calcula o indice do vetor, como se fosse uma matriz.
return i*Ny + j;
}
-
+/*
void calcU(int i, int j, double *u) {
double res = 0;
res += f(i,j) + (u[ in(i+1,j) ] + u[ in(i-1,j) ] ) / Hx * Hx + (u[ in(i,j+1) ] + u[ in(i,j-1) ]) / Hy * Hy;
@@ -52,6 +85,29 @@ void calcU(int i, int j, double *u) {
res = res / UDivisor;
u[ in(i,j) ] = res;
}
+*/
+double calcU(int n, double *u) {
+/*
+u(i,j) = f(x,y) + (u(i+1,j) + u(i-1,j))/Δx² + (u(i,j+1) + u(i,j-1))/Δy² + (-u(i+1,j)+u(i-1,j))/2Δx + (-u(i,j+1)+u(i,j-1))/2Δy
+ --------------------------------------------------------------------------------------------------------------------
+ 2/Δx² + 2/Δy² + 4π²
+*/
+ double res = 0;
+ res += f(n) + (u[n+Ny] + u[n-Ny] ) / Hx * Hx + (u[n+1] + u[n-1]) / Hy * Hy;
+ res += (u[n-Ny] - u[n+Ny]) / 2 * Hx + (u[n-1] - u[n+1]) / 2 * Hy;
+ res = res / UDivisor;
+ //u[n] = res;
+ return res;
+}
+
+double calcRes(int n, double *u) {
+//f(x,y) = (2/Δx² + 2/Δy² + 4π²) * u(i,j) - ((u(i+1,j) + u(i-1,j))/Δx² + (u(i,j+1) + u(i,j-1))/Δy² + (-u(i+1,j)+u(i-1,j))/2Δx + (-u(i,j+1)+u(i,j-1))/2Δy)
+ double res = 0;
+ //res = (2/Hx*Hx + 2/Hy*Hy + 4*M_PI*M_PI) * u[n];
+ res = UDivisor * u[n];
+ res -= ((u[n+Ny] + u[n-Ny] ) / Hx * Hx + (u[n+1] + u[n-1]) / Hy * Hy + (u[n-Ny] - u[n+Ny]) / 2 * Hx + (u[n-1] - u[n+1]) / 2 * Hy);
+ return res;
+}
double residueNorm(double *r) {
// u is only a row and r is the last residue
@@ -63,64 +119,145 @@ double residueNorm(double *r) {
return sqrt(res);
}
+void generate_matrix(double *A, double *b) {
+ int i,j;
+ for(i=0; i<N*N; ++i) {
+ A[i] = 0;
+ }
+ for(i=0; i<N; ++i) {
+ for(j=0; j<N; ++j) {
+ if(i == j)
+ A[in(i,j)] = -4;
+ }
+ }
+ A[in(0,1)] = 1;
+ A[in(0,3)] = 1;
+ A[in(1,0)] = 1;
+ A[in(1,2)] = 1;
+ A[in(1,4)] = 1;
+ A[in(2,1)] = 1;
+ A[in(2,5)] = 1;
+ A[in(3,0)] = 1;
+ A[in(3,4)] = 1;
+ A[in(3,6)] = 1;
+ A[in(4,1)] = 1;
+ A[in(4,3)] = 1;
+ A[in(4,5)] = 1;
+ A[in(4,7)] = 1;
+ A[in(5,2)] = 1;
+ A[in(5,4)] = 1;
+ A[in(5,8)] = 1;
+ A[in(6,3)] = 1;
+ A[in(6,7)] = 1;
+ A[in(7,4)] = 1;
+ A[in(7,6)] = 1;
+ A[in(7,8)] = 1;
+ A[in(8,5)] = 1;
+ A[in(8,7)] = 1;
+ b[0] = -100;
+ b[1] = -20;
+ b[2] = -20;
+ b[3] = -80;
+ b[4] = 0;
+ b[5] = 0;
+ b[6] = -260;
+ b[7] = -180;
+ b[8] = -180;
+}
+
+
+/*
+a11 * x1 + a12 * x2 + a13 * x3 = b1 -> x1 = (b1 - (a12 * x2 + a13 * x3)) / a11
+x1 = u11;
+x2 = u12;
+x3 = u13;
+x4 = u21;
+x5 = u22;
+x6 = u23;
+
+Res[1] = b1 - (a11 * x1 + a12 * x2 + a13 * x3)
+*/
+
+
+void sor(double *b, double *x, double *r, double *timeSor, double *timeResNorm) {
+ int i,j,k; // Ax = b -> ?u = f
+ double sigma, now, fxy, res, maxRes = 0, tRes = 0; // maxRes is the biggest residue, tRes is total residue in this iteration.
+ for(k=0; k<MaxI; k++) {
+ now = timestamp();
+ for(i=0; i<Nx * Ny; ++i) {
+ //b[i] = f();
+ /*sigma = 0;
+ for(j=0; j<N; ++j) {
+ sigma = sigma + A[in(i,j)] * x[j];
+ }
+ sigma = sigma - A[in(i,i)] * x[i];*/
+ // b[1] - somatorio / a11
+ //x[i] = x[i] + W * ((b[i] - sigma)/A[in(i,i)] - x[i]);
+ x[i] = x[i] + W * (calcU(i,x) - x[i]);
+ //res = b[i] - x[i];
+ }
+ *timeSor += timestamp() - now;
+ now = timestamp();
+ res = f(k);
+ printf("Res = %lf - ",res);
+ /*for(i = 0; i < Nx * Ny; ++i)
+ res -= A[in(i,i)] * x[i];*/
+ res -= calcRes(k,x);
+ printf("Res = %lf\n",res);
+ if(res > maxRes)
+ maxRes = res;
+ tRes += res * res;
+ //residueNorm(r);
+ r[k] = sqrt(tRes);
+ tRes = 0;
+ *timeResNorm += timestamp() - now;
+ }
+ *timeSor = *timeSor / MaxI;
+ *timeResNorm = *timeResNorm / MaxI;
+}
+
int main(int argc, char *argv[]) {
int i, j, k;
FILE *fpExit;
getParams(argc,argv,fpExit);
- Nx = (1/Hx);
- Ny = (1/Hy);
+ Nx = (M_PI/Hx);
+ Ny = (M_PI/Hy);
Nx++;
Ny++;
- if(Nx != Ny) {
- puts("Not a square matrix. Result will be undetermined or impossible.");
- exit(-1);
- }
W = (2 - (Hx + Hy)) / 2;
- if(W < 1) {
- puts("Relaxation factor less than 1 - program may not work as you expect.");
- }
UDivisor = (2 / Hx * Hx) + (2 / Hy * Hy) + 4 * Pipi;
- /*if((Nx * Ny * 8) <= MAX_SIZE) { // Should we save a matrix with f(x,y) values?
- inMemory = 1;
- }*/
double sigma;
- //double **u,**f;
+ //double *A = calloc(sizeof(double), Nx * Ny * Nx * Ny);
+ double *b, *x, *r, *timeSor, *timeResNorm;
- /*for(i=0; i<Nx; i++) {
- for(j=0; i<Ny; i++) {
- u[i][j] = calloc(sizeof(double));
- f[i][j] = malloc(sizeof(double));
- f[i][j] = f(i*Hx,j*Hy);
- }
- }*/
+ //generate_matrix(A,b);
- double *u = calloc(sizeof(double), Nx * Ny);
+ b = malloc(Nx * Ny * sizeof(double));
+ x = malloc(Nx * Ny * sizeof(double));
+ r = malloc(MaxI * sizeof(double));
+ timeSor = calloc(1,sizeof(double));
+ timeResNorm = calloc(1,sizeof(double));
sigma = sinh(M_PI * M_PI);
- for(i=0; i<Nx; i++) { // Calculate side extremities (top and bottom are always 0).
- u[ in(i,0) ] = sin(2 * M_PI * (M_PI - (i * Hx))) * sigma;
- u[ in(i,Nx) ] = sin(2 * M_PI * (i * Hx));
- }
+ /*for(i=0; i<Nx; ++i) { // Calculate side extremities (top and bottom are always 0, so calloc already took care of it).
+ A[ in(i,0) ] = sin(2 * M_PI * (M_PI - (i * Hx))) * sigma;
+ A[ in(i,Nx) ] = sin(2 * M_PI * (i * Hx));
+ }*/
- // Nx columns and Ny rows.
+ //print_matrix(A);
- /*for(k=0; k<MaxI; k++) {
- for(i=1; i<Ny; i++) {
- sigma = 0;
- for(j=1; j<Nx; j++) {
- /*if(j != i) {
- sigma += a[i][j] * fi[i];
- }*/
- /*sigma += a[i][j] * fi[i];
- }
- sigma -= a[i][i] * fi[i];
- //fi[i] = (1 - W) * fi[i] + W/a[i][i] * (b[i] - sigma);
- fi[i] = fi[i] + W * ((b[i] - sigma) / a[i][i] - fi[i])
- }
- }*/
+ sor(b,x,r,timeSor,timeResNorm);
+
+ print_vector(x);
+
+ printf("TimeSor: %lf\nTimeResNorm: %lf\n\nNorma do Resíduo\n",*timeSor,*timeResNorm);
+ for(i=0; i<MaxI; ++i) {
+ printf("i=%d: %lf\n",i,r[i]);
+ }
+ // Nx columns and Ny rows.
return 0;
}
\ No newline at end of file
diff --git a/sor.c b/sor.c
new file mode 100644
index 0000000000000000000000000000000000000000..779bf5aa0bf5b38f3199e87b8a9d5d801bec22f3
--- /dev/null
+++ b/sor.c
@@ -0,0 +1,124 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <errno.h>
+
+#define SIZE 9
+
+int N = SIZE;
+int MaxI = 10;
+
+inline int in(int i, int j) {
+ return i*N + j;
+}
+
+void generate_matrix(double *A, double *b) {
+/* A[in(0,0)] = 10;
+ A[in(0,1)] = 4;
+ A[in(0,2)] = 2;
+ A[in(1,0)] = -2;
+ A[in(1,1)] = 15;
+ A[in(1,2)] = 3;
+ A[in(2,0)] = 2;
+ A[in(2,1)] = -2;
+ A[in(2,2)] = 13;
+ b[0] = 5;
+ b[1] = 3;
+ b[2] = 4;*/
+ int i,j;
+ for(i=0; i<N*N; ++i) {
+ A[i] = 0;
+ }
+ for(i=0; i<N; ++i) {
+ for(j=0; j<N; ++j) {
+ if(i == j)
+ A[in(i,j)] = -4;
+ }
+ }
+ A[in(0,1)] = 1;
+ A[in(0,3)] = 1;
+ A[in(1,0)] = 1;
+ A[in(1,2)] = 1;
+ A[in(1,4)] = 1;
+ A[in(2,1)] = 1;
+ A[in(2,5)] = 1;
+ A[in(3,0)] = 1;
+ A[in(3,4)] = 1;
+ A[in(3,6)] = 1;
+ A[in(4,1)] = 1;
+ A[in(4,3)] = 1;
+ A[in(4,5)] = 1;
+ A[in(4,7)] = 1;
+ A[in(5,2)] = 1;
+ A[in(5,4)] = 1;
+ A[in(5,8)] = 1;
+ A[in(6,3)] = 1;
+ A[in(6,7)] = 1;
+ A[in(7,4)] = 1;
+ A[in(7,6)] = 1;
+ A[in(7,8)] = 1;
+ A[in(8,5)] = 1;
+ A[in(8,7)] = 1;
+ b[0] = -100;
+ b[1] = -20;
+ b[2] = -20;
+ b[3] = -80;
+ b[4] = 0;
+ b[5] = 0;
+ b[6] = -260;
+ b[7] = -180;
+ b[8] = -180;
+}
+
+void print_vector(double *x) {
+ int i;
+ for(i=0; i<N; ++i) {
+ printf("%f ", x[i]);
+ }
+ printf("\n");
+}
+
+void print_matrix(double *A) {
+ int i,j;
+ for(i=0; i<N; ++i) {
+ for(j=0; j<N; ++j) {
+ printf("%.0f ", A[in(i,j)]);
+ }
+ printf("\n");
+ }
+ printf("\n\n\n");
+}
+
+void sor(double *A, double *b, double *x, double w) {
+ int i,j,k;
+ double sigma;
+ for(k=0; k<MaxI; k++) {
+ for(i=0; i<N; ++i) {
+ sigma = 0;
+ for(j=0; j<N; ++j) {
+ sigma = sigma + A[in(i,j)] * x[j];
+ }
+ sigma = sigma - A[in(i,i)] * x[i];
+ x[i] = x[i] + w * ((b[i] - sigma)/A[in(i,i)] - x[i]);
+ }
+ }
+}
+
+int main(int argc, char* argv[]) {
+ double *A,*b,*x,w;
+ int i,j;
+ A = malloc(N * N * sizeof(double));
+ b = malloc(N * sizeof(double));
+ x = malloc(N * sizeof(double));
+ generate_matrix(A,b);
+ for(i=0; i<N; ++i) {
+ x[i] = 1;
+ }
+ print_vector(x);
+ print_matrix(A);
+ print_vector(b);
+ w = 1.25;
+ sor(A,b,x,w);
+
+ print_vector(x);
+}
\ No newline at end of file