Προς το περιεχόμενο

Μαθηματικα και προγραμματισμος


NiaZ

Προτεινόμενες αναρτήσεις

Δημοσ.

Εχω μια εργασια στη σχολη και χρειαζομαι μια μικρη βοηθεια τα παντα ειναι ευπροσδεκτα(LINKS,CODE,etc..)

Οποιος γνωριζει κατι σχετικα με Jacobi σε παραλληλο κωδικα. Αχ=Β.

ειτε σε Fortran 77 ειτε σε C.

 

Ευχαριστω εκ των προτερων.

Δημοσ.

to vrika kai se psifiaki morfh

C PROGRAM JACOB

C

C JACOBI'S METHOD

C

PARAMETER(ND=10)

DOUBLE PRECISION EPS,A,B,X,Y,E

DIMENSION A(ND,ND),B(ND),X(ND),Y(ND)

READ(5,*) N,NIT,EPS

READ(5,*) ((A(I,J),J=1,N),I=1,N)

READ(5,*) (B(I),I=1,N)

READ(5,*) (X(I),I=1,N)

WRITE(6,*) 'INPUT:'

WRITE(6,*) 'N=',N,' NIT=',NIT,' EPS=',EPS

WRITE(6,*) 'A='

DO 10 I=1,N

WRITE(6,*) (A(I,J),J=1,N)

10 WRITE(6,*) ' '

WRITE(6,*) 'B=',(B(I),I=1,N)

WRITE(6,*) 'X=',(X(I),I=1,N)

WRITE(6,*) ' '

WRITE(6,*) 'OUTPUT:'

CALL JACOBI(ND,N,A,B,NIT,EPS,X,KIT,E,ITEST,Y)

WRITE(6,*) 'ITEST=',ITEST,' KIT=',KIT

WRITE(6,*) 'X=',(X(I),I=1,N)

WRITE(6,*) 'E=',E

STOP

END

C-----------------------------------------------------------------

SUBROUTINE JACOBI(ND,N,A,B,NIT,EPS,X,KIT,E,ITEST,Y)

C

C SOLUTION OF THE LINEAR SYSTEM AX=B BY JACOBI'S METHOD

C

C INPUT:

C ND: MAXIMUM COLUMN DIMENSION OF ARRAYS

C N: DIMENSION OF ARRAYS

C A: SYSTEM MATRIX (NXN)

C B: RIGHT-HAND SIDE VECTOR (N)

C NIT: MAXIMUM NUMBER OF ITERATIONS

C EPS: CONVERGENCE TEST NUMBER

C

C OUTPUT:

C X: APPROXIMATE SOLUTION VECTOR (N)

C KIT: NUMBER OF EXECUTED ITERATIONS

C E: ERROR BOUND

C ITEST = 1, CONVERGENCE TEST SATISFIED

C = 2, CONVERGENCE TEST NOT SATISFIED

C

DOUBLE PRECISION EPS,A,B,X,Y,CN,S,E

DIMENSION A(ND,N),B(ND),X(ND),Y(ND)

ITEST=1

C

C ITERATION MATRIX NORM

C

CN=0.D0

DO 10 I=1,N

S=0.D0

DO 20 J=1,N

IF(J.EQ.I) GO TO 20

S=S+DABS(A(I,J))

20 CONTINUE

S=S/DABS(A(I,I))

IF(S.GT.CN) CN=S

10 CONTINUE

C

C ITERATIONS

C

DO 30 K=1,NIT

KIT=K

DO 60 I=1,N

60 Y(I)=X(I)

E=0.D0

DO 40 I=1,N

X(I)=B(I)

DO 50 J=1,N

IF(J.EQ.I) GO TO 50

X(I)=X(I)-A(I,J)*Y(J)

50 CONTINUE

X(I)=X(I)/A(I,I)

S=DABS(X(I)-Y(I))

IF(S.GT.E) E=S

40 CONTINUE

WRITE(6,*) ' X=',(X(I),I=1,N)

C

C ERROR TEST

C

E=E*CN/(1.D0-CN)

IF(E.LE.EPS) RETURN

30 CONTINUE

ITEST=2

RETURN

END

Δημοσ.
powerfty exeis tipota na vohthiseis?

 

Agapite NiaZ syggnwmi gia tin kathysterisi alla etyxe kati epeigon, h apofoitisi mou! :D Exw na sou dwsw pros to paron parallilo kwdika MPI se gia tin methodo Gauss - Jordan, kapou exw kai tin Jacobi ylopoihmeni alla den tin vriskw, wstoso einai paromoies autes des tin mia na pareis mia idea.

 

#include "mpi.h"

#include <stdio.h>

#include <stdlib.h>

 

void

main(int argc, char *argv[])

{

int i, n, j, k;

int type = 1;

int me, numprocs;

int count = 0;

int st;

int temp;

double c, d;

 

MPI_Status status;

 

double** ap;

double* a;

 

FILE* finput;

FILE* foutput;

 

int namelen;

char processor_name[MPI_MAX_PROCESSOR_NAME];

 

MPI_Init(&argc,&argv);

MPI_Comm_size(MPI_COMM_WORLD,&numprocs);

MPI_Comm_rank(MPI_COMM_WORLD,&me);

MPI_Get_processor_name(processor_name,&namelen);

 

fprintf(stderr,"Process %d on %s\n",me, processor_name);

fflush(stderr);

 

if(me == 0) {

finput = fopen("input.txt", "r");

foutput = fopen("output.txt", "w");

 

fscanf(finput, "%d", &st);

}

 

MPI_Bcast(&st, 1, MPI_INT, 0, MPI_COMM_WORLD);

 

if(me == 0) {

 

n = st;

 

ap = (double**)malloc(n / numprocs * sizeof(double));

 

for(i = 0; i < n / numprocs; ++i)

ap = (double*)malloc((n + 1) * sizeof(double));

 

a = (double*)malloc((n + 1) * sizeof(double));

 

i = 0;

while((fscanf(finput, "%lg", &d)) && (count < st * (st + 1))) {

fflush(stdin);

temp = count / (st + 1);

i = temp / (n / numprocs);

a[count % (st + 1)] = d;

 

if(i == me)

ap[temp][count % (st + 1)] = d;

 

if(count % (st + 1) == st) {

 

if(i != me)

MPI_Send(a, n + 1, MPI_DOUBLE, i, type, MPI_COMM_WORLD);

}

 

++count;

}

 

for(k = 0; k < n; ++k) {

temp = k / (n / numprocs);

if(temp == me) {

for(j = 1; j < numprocs; ++j) {

MPI_Send(ap[k], n + 1, MPI_DOUBLE, j, type, MPI_COMM_WORLD);

}

}

else {

MPI_Recv(a, n + 1, MPI_DOUBLE, temp, type, MPI_COMM_WORLD, &status);

}

if(k < n / numprocs)

d = ap[k][k];

else

d = a[k];

 

if(k < n / numprocs) {

for(j = 0; j < k; ++j) {

if(d != 0) {

c = -ap[j][k] / d;

ap[j][k] = 0.0;

for(i = k + 1; i < n + 1; ++i) {

ap[j] += c * ap[k];

}

}

}

 

for(j = k + 1; j < n / numprocs; ++j) {

if(d != 0) {

c = -ap[j][k] / d;

ap[j][k] = 0.0;

for(i = k + 1; i < n + 1; ++i) {

ap[j] += c * ap[k];

}

}

}

}

 

else {

for(j = n / numprocs; j < n; ++j) {

if(d != 0) {

temp = j / (n / numprocs);

c = -ap[j - temp * n / numprocs][k] / d;

ap[j - temp * n / numprocs][k] = 0.0;

for(i = k + 1; i < n + 1; ++i) {

ap[j - temp * n / numprocs] += c * a;

}

}

}

}

}

 

for(i = 0; i < n / numprocs; ++i) {

printf("Process %d --> ", me);

fflush(stdout);

for(j = 0; j < n + 1; ++j) {

printf("%g ", ap[j]);

fflush(stdout);

}

printf("\n");

fflush(stdout);

}

 

for(i = 0; i < n / numprocs; ++i) {

for(j = 0; j < n + 1; ++j) {

fprintf(foutput, "%g ", ap[j]);

fflush(foutput);

}

fprintf(foutput, "\t %g\n", ap[n] / ap);

fflush(foutput);

}

for(k = 1; k < numprocs; ++k) {

for(i = 0; i < n / numprocs; ++i) {

MPI_Recv(a, n + 1, MPI_DOUBLE, k, type, MPI_COMM_WORLD, &status);

for(j = 0; j < n + 1; ++j) {

fprintf(foutput, "%g ", a[j]);

fflush(foutput);

}

fprintf(foutput, "\t %g\n", a[n] / a[k * n / numprocs + i]);

fflush(foutput);

}

}

}

else {

 

n = st;

a = (double*)malloc((n + 1) * sizeof(double));

 

ap = (double**)malloc(n / numprocs * sizeof(double));

 

for(i = 0; i < n / numprocs; ++i)

ap = (double*)malloc((n + 1) * sizeof(double));

 

for(i = 0; i < n / numprocs; ++i) {

MPI_Recv(ap, n + 1, MPI_DOUBLE, 0, type, MPI_COMM_WORLD, &status);

}

 

for(k = 0; k < n; ++k) {

temp = k / (n / numprocs);

if(temp == me) {

for(j = 0; j < numprocs; ++j) {

if(j != me) {

MPI_Send(ap[k - me * n / numprocs], n + 1, MPI_DOUBLE, j, type, MPI_COMM_WORLD);

}

}

}

else {

MPI_Recv(a, n + 1, MPI_DOUBLE, temp, type, MPI_COMM_WORLD, &status);

}

 

if((k < (me + 1) * n / numprocs) && (k >= me * n / numprocs))

d = ap[k - me * n / numprocs][k];

else

d = a[k];

 

if((k < (me + 1) * n / numprocs) && (k >= me * n / numprocs)) {

for(j = me * n / numprocs; j < k; ++j) {

if(d != 0) {

temp = j / (n / numprocs);

c = -ap[j - temp * n / numprocs][k] / d;

ap[j - temp * n / numprocs][k] = 0.0;

for(i = k + 1; i < n + 1; ++i) {

ap[j - temp * n / numprocs] += c * ap[k - temp * n / numprocs];

}

}

}

 

for(j = k + 1; j < (me + 1) * n / numprocs; ++j) {

if(d != 0) {

temp = j / (n / numprocs);

c = -ap[j - temp * n / numprocs][k] / d;

ap[j - temp * n / numprocs][k] = 0.0;

 

for(i = k + 1; i < n + 1; ++i) {

ap[j - temp * n / numprocs] += c * ap[k - temp * n / numprocs];

}

}

}

}

 

else if(k < me * n / numprocs) {

for(j = 0; j < me * n / numprocs; ++j) {

if(d != 0) {

temp = j / (n / numprocs);

c = -ap[j - temp * n / numprocs][k] / d;

ap[j - temp * n / numprocs][k] = 0.0;

 

for(i = k + 1; i < n + 1; ++i) {

ap[j - temp * n / numprocs] += c * a;

}

}

}

}

 

else if(k >= (me + 1) * n / numprocs) {

 

for(j = (me + 1) * n / numprocs; j < n; ++j) {

if(d != 0) {

temp = j / (n / numprocs);

c = -ap[j - temp * n / numprocs][k] / d;

ap[j - temp * n / numprocs][k] = 0.0;

 

for(i = k + 1; i < n + 1; ++i) {

ap[j - temp * n / numprocs] += c * a;

}

}

}

}

}

 

for(i = 0; i < n / numprocs; ++i) {

printf("Process %d --> ", me);

fflush(stdout);

for(j = 0; j < n + 1; ++j) {

printf("%g ", ap[j]);

fflush(stdout);

}

printf("\n");

fflush(stdout);

}

 

for(i = 0; i < n / numprocs; ++i)

MPI_Send(ap, n + 1, MPI_DOUBLE, 0, type, MPI_COMM_WORLD);

 

}

 

MPI_Finalize();

}

 

An exeis kapoia aporia min distaseis na me rwtiseis :) Kales giortes!!

 

P.S.: Stin parousa ylopoihsh o arithmos twn grammwn tou pinaka A prepei na einai pollaplasio tou arithmou twn epeksergastwn gia na min exoyme apwleies stous ypologismous.

Δημοσ.

Να πώ και γώ ότι υπάρχει το numerical recipes in fortran ή C στην διεύθυνση:

 

http://www.nr.com/

 

Έχει θεωρία και έτοιμες υπορουτίνες για χρήση. Τα βιβλία αντίστοιχα fοrtran ή C είναι κάπου 20ΜΒ το καθένα (δωρεάν πάντα).

 

Επίσης, ξέρει κανείς πού κυκλοφορούν έτοιμες υπορουτίνες για το MPI σε fortran?

Αρχειοθετημένο

Αυτό το θέμα έχει αρχειοθετηθεί και είναι κλειστό για περαιτέρω απαντήσεις.

  • Δημιουργία νέου...