c***************************************************************** c c Example 26 gauss_elimination.f c c c Objective : To solve the system of Linear Equations c using Gaussian Elimination without pivoting c on 'p' processors. c c Input : Read files (mdatgaus.inp) for Matrix A c and (vdatgaus.inp) for Vector b c c Output : The solution of matrix system of linear c equations Ax=b on processor 0. c c Description : Input matrix is stored in n by n format. c Columnwise cyclic distribution of input c matrix is used for partitioning of the c matrix. c c Necessary conditions : Number of Processes should be less than c or equal to 8.Matrix size for Matrix A c and vector size for vector b should be c equally striped, that is Matrix size and c Vector Size should be divisible by Number c of processes used. c c****************************************************************** c c c program main include 'mpif.h' PARAMETER ( NMAX = 100) integer MyRank,Numprocs,NoofRows,nlocal_proc,myrow, ierr integer recv_send,send_recv_tag integer idest, NoofCols integer i,j,k,STATUS(MPI_STATUS_SIZE) real Matrix_A,Input_B,Cols_MatA,b1,y,Buffer_Pivot, $ X_buffer,x1 dimension Matrix_A(NMAX,NMAX),Input_B(NMAX),Cols_MatA(NMAX), $ Buffer_Pivot(NMAX), X_buffer(NMAX) call MPI_INIT( ierr ) call MPI_COMM_RANK( MPI_COMM_WORLD, MyRank, ierr ) call MPI_COMM_SIZE( MPI_COMM_WORLD, Numprocs, ierr ) c c Read matrices A,B at process #0 c if ( MyRank .eq. 0) then open(unit=12, file = './data/mdatgaus.inp') read(12,*) NoofRows, NoofCols c print *, NoofRows, NoofCols write(6,*) 'Input Matrix' do i = 1,NoofRows read(12,*) (Matrix_A(i,j),j=1,NoofCols) write(6,75) (Matrix_A(i,j),j=1,NoofCols) enddo 75 format(8(2x,f8.3)) open(unit=13, file = './data/vdatgaus.inp') write(6,*) 'Input Vector' read(13,*) NoofRows read(13,*) (Input_B(i), i =1,NoofRows) write(6,75) (Input_B(i), i = 1, NoofRows) close(12) close(13) endif c Broadcast the number of rows of Matrix A c call MPI_BCAST(NoofRows,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) if(mod(NoofRows,Numprocs) .ne. 0 ) then if(MyRank .eq. 0) $ print*,"Matrix cannot be evenly striped among processes" goto 30 endif c Distribute A,B to other nodes c Cyclic rowwise distribution nlocal_proc = NoofRows/Numprocs if ( MyRank .eq. 0 ) then do i = 1,NoofRows do j = 1,NoofRows Cols_MatA(j) = Matrix_A(i,j) enddo b1 = Input_B(i) idest = mod(i,Numprocs) recv_send = (i/Numprocs+1)*2 if ( idest .ne. 0) then c c Node 0 send corresponding row of c matrix A and element of B to other nodes c c call MPI_SEND( Cols_MatA,NoofRows, MPI_REAL,idest, $ recv_send, MPI_COMM_WORLD,STATUS,ierr) send_recv_tag = recv_send + 1 c call MPI_SEND( b1,1, MPI_REAL,idest,send_recv_tag, $ MPI_COMM_WORLD,STATUS,ierr) endif enddo else c do recv_send = 2, nlocal_proc*2,2 c Receive information from node 0 call MPI_RECV( Cols_MatA,NoofRows, $ MPI_REAL,0,recv_send,MPI_COMM_WORLD,STATUS,ierr) send_recv_tag = recv_send + 1 call MPI_RECV( b1,1, $ MPI_REAL,0,send_recv_tag,MPI_COMM_WORLD,STATUS,ierr) myrow = MyRank + (recv_send/2-1)*Numprocs do k = 1, NoofRows Matrix_A(myrow,k) = Cols_MatA(k) enddo Input_B(myrow) = b1 enddo endif c c c Gaussian_elimination do k = 1, NoofRows c Decide whose turn it is to do normalize A(k,k) c if (MyRank .eq. mod(k,Numprocs)) then do j = k+1, NoofRows Matrix_A(k,j) = Matrix_A(k,j) / Matrix_A(k,k) Buffer_Pivot(j) = Matrix_A(k,j) enddo Input_B(k) = Input_B(k)/Matrix_A(k,k) y = Input_B(k) Matrix_A(k,k) = 1.0 endif c Broadcast the coefficients Matrix_A and result Input_B c call MPI_BCAST(y,1,MPI_REAL,mod(k,Numprocs), $ MPI_COMM_WORLD,ierr) call MPI_BCAST(Buffer_Pivot,NoofRows,MPI_REAL,mod(k,Numprocs), $ MPI_COMM_WORLD,ierr) do i = k+1,NoofRows if (MyRank .eq. mod(i,Numprocs)) then c update the coefficient Matrix_A and result Input_B do j = k+1,NoofRows Matrix_A(i,j) = Matrix_A(i,j) - $ Matrix_A(i,k)*Buffer_Pivot(j) enddo Input_B(i) = Input_B(i) - Matrix_A(i,k)*y Matrix_A(i,k) = 0.0 endif enddo enddo c c Triangle Solver do k = NoofRows, 1, -1 c Decide whose turn to solve the X_buffer(k) value if (MyRank .eq. mod(k,Numprocs)) then X_buffer(k) = Input_B(k)/Matrix_A(k,k) x1 = X_buffer(k) endif c Broadcast this X_buffer(k) value c Update result Input_B and eliminate X_buffer(k) term if ( k .ne. 0) $ call MPI_BCAST(x1,1,MPI_REAL,mod(k,Numprocs), $ MPI_COMM_WORLD,ierr) do i = k, 1, -1 if (MyRank .eq. mod(i,Numprocs)) then Input_B(i) = Input_B(i) - Matrix_A(i,k)*x1 endif c Accumulate X_buffer(k) value in node 0 if (MyRank .eq. 0) X_buffer(k) = x1 enddo enddo cc Node 0 output the final result X_buffer. c if (MyRank .eq. 0) then write(6,*) write(6,*) 'Results for solution of matrix system of equations' do i = 1,NoofRows write(6,97) i,X_buffer(i) 97 format(1x,'x[',i1,'] = ',f5.2) enddo endif 30 call MPI_FINALIZE(ierr) stop end c *************************************************** c