program PCMRH !*************************************************************** ! Parallel CMRH program with MPI !*************************************************************** use mpi external SSWAP,DCOPY,DSCAL,DGEMV,DSWAP intrinsic dabs,dsqrt !******************************************** ! VARIABLES !******************************************** integer,parameter ::Nprocs_max=16 integer ::rang,Nprocs,n,nl,code,i,j,k,choix,solchoice,iprint real(kind=8) ::init(5) integer,dimension(MPI_STATUS_SIZE) ::statut real(kind=8),allocatable,dimension(:,:) ::AL real(kind=8) ::btemp,bmax,btemp2 real(kind=8),allocatable,dimension(:) ::b,sol,rs,temp,x,workl,work,tempn,c,s,xl,cl,temp2,xx,yy,valeurs integer,allocatable,dimension(:) ::p real*8 ::dnorm2,D_inv,beta,dnorm,dnrm2,tol,alpha,H_k1,dgam,r integer ::indice,idamax,k1,n_k,knl,knlp,inl,inlp real*8 dtime, delta, delta0, dt1, dt2, dt3 real*8 etime, eelta, eelta0, et1, et2, et3 real*8 t1, t2, t3, tm1, tm2, tm3 real*8 dtime1, delta1, delta01 real*8 etime1, eelta1, eelta01,l,p1,p2 integer nb_ligne real(kind=8) ::temps,tempsd,tempsf,temps2 integer*4,dimension(3) ::timearray real ::rand real*4 da(2),ea(2),da1(2),ea1(2) !******************************************** ! MPI CONFIG !******************************************** call MPI_INIT(code) call MPI_COMM_RANK(MPI_COMM_WORLD,rang,code) call MPI_COMM_SIZE(MPI_COMM_WORLD,Nprocs,code) call MPI_BARRIER(MPI_COMM_WORLD,code) call itime(timearray) ! Get the current time i = rand ( timearray(1)+timearray(2)+timearray(3) ) !************************************ ! Matrix Initialisation !************************************ if (rang==0) then open(unit=10,file='inputfile.dat',status='old') read(10,*)init close(10) iprint=int(init(1)) n=int(init(2)) choix=int(init(3)) solchoice=int(init(4)) tol=init(5) end if tempsd=MPI_WTIME() call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,code) call MPI_BCAST(choix,1,MPI_INTEGER,0,MPI_COMM_WORLD,code) call MPI_BCAST(tol,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,code) if (mod(n,Nprocs)==0)then nl=n/Nprocs else nl=(n/Nprocs)+1 end if allocate(b(n),rs(n),x(n),AL(nl,n),workl(nl),sol(n),temp(n),temp2(n)) if (choix<3)then call matrice(AL,nl,n,2.d0,1.d0,choix,rang) else allocate(valeurs(n)) open(unit=30,file='matrixfile.dat',status='old') do i=1,n call position(i,nl,inl,inlp) read(30,*)valeurs if (rang.eq.inl)then AL(inlp,:)=valeurs(:) endif enddo close(30) endif !******************************************** ! End Matrix initialisation !******************************************** !******************************************** ! Solution initialisation !******************************************** call MPI_BCAST(solchoice,1,MPI_INTEGER,0,MPI_COMM_WORLD,code) call solution(AL,workl,sol,nl,n,solchoice) call MPI_ALLGATHER(workl,nl,MPI_DOUBLE_PRECISION,b,nl,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,code) !******************************************** ! End solution initialisation !******************************************** !******************************************** ! Permutation Vector initialisation ! and computation of initial residual !******************************************** if (rang==0) then allocate(p(n)) do i=1,n p(i)=i end do end if dnorm2=dnrm2(n,b,1) if (rang==0)then if (iprint==2)then print*,'' print*,'la norme 2 de r0=',dnorm2 print*,'' end if end if !******************************************** ! Computation of max b(i) !******************************************** indice = idamax(n,b,1) beta = b(indice) D_inv = 1.d0/beta !******************************************** ! First permutation !******************************************** if (indice .ne.1) then call dswap(1,b(indice),1,b(1),1) call dswap(nl,AL(1,indice),1,AL(1,1),1) call position(indice,nl,inl,inlp) if (inl==0) then if (rang==inl)then call dswap(n,AL(inlp,1),nl,AL(1,1),nl) end if else if (rang==inl) then call dcopy(n,AL(inlp,1),nl,temp,1) end if call MPI_BCAST(temp,n,MPI_DOUBLE_PRECISION,inl,MPI_COMM_WORLD,code) if (rang==0)then call dcopy(n,AL(1,1),nl,temp2,1) end if call MPI_BCAST(temp2,n,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,code) if (rang==inl)then call dcopy(n,temp2,1,AL(inlp,1),nl) end if if (rang==0)then call dcopy(n,temp,1,AL(1,1),nl) end if end if if (rang==0)then call sswap(1,p(indice),1,p(1),1) end if end if !******************************************** ! End permutation !******************************************** ! call dscal(n,0.d0,x,1) call dscal(n,D_inv,b,1) rs(1)=beta dnorm=dabs(beta) allocate(tempn(nl),work(n),c(n),s(n),xl(nl),cl(nl)) call dscal(nl,0.d0,xl,1) call MPI_BARRIER(MPI_COMM_WORLD,code) k=0 if (rang==0) then if (iprint==2)then write(*,*)'iter., pseudo residu ' end if open(unit=30,file='presidual.dat',status='unknown') write(30,*)'Iter., Pseudo Residual Norm' end if !tol=1.d-7 delta0 = dtime(da) eelta0 = etime(ea) tempsf=MPI_WTIME() !******************************************** ! Loop !******************************************** do while((dnorm.gt.tol).and.(k.lt.n)) k=k+1; n_k=n-k; k1=k+1; if (rang==0) then if (iprint==2)then write(*,10)k-1,dnorm end if write(30,10)k-1,dnorm end if 10 format(i5,2x,e12.4) !******************************************** ! workl = A*b !******************************************** call dgemv('N',nl,n_k+1,1.d0,AL(1,k),nl,b(k),1,0.d0,workl,1) ! ! call dscal(nl,0.d0,workl,1) ! ! do j=k1,n ! ! do i=1,nl ! ! workl(i)=workl(i)+AL(i,j)*b(j) ! ! end do ! ! end do !******************************************** ! A(:,k)=b !******************************************** call position(k1,nl,knl,knlp) if (rang>knl)then call dcopy(nl,b(rang*nl+1),1,AL(1,k),1) else if (rang==knl) then call dcopy(nl-knlp+1,b(k1),1,AL(knlp,k),1) end if end if !******************************************** ! j Loop !******************************************** do j=1,k !******************************************** ! A(j,k)=workl(j) ! workl(j)=0 !******************************************** call position(j,nl,knl,knlp) if (rang==knl) then AL(knlp,k)=workl(knlp) workl(knlp)=0.d0 alpha=-AL(knlp,k) end if call MPI_BCAST(alpha,1,MPI_DOUBLE_PRECISION,knl,MPI_COMM_WORLD,code) !******************************************** ! workl=workl - alpha A(:,j) !******************************************** call position(j+1,nl,knl,knlp) if (rang >knl)then call daxpy(nl,alpha,AL(1,j),1,workl,1) end if if (rang==knl) then call daxpy(nl-knlp+1,alpha,AL(knlp,j),1,workl(knlp),1) end if end do !******************************************** ! End j Loop !******************************************** if (k.lt.n)then call MPI_ALLGATHER(workl,nl,MPI_DOUBLE_PRECISION,b,nl,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,code) !******************************************** ! Search indice i_0 !******************************************** indice=k+idamax(n_k,b(k1),1) H_k1=b(indice) !******************************************** ! Permutation !******************************************** if (indice.ne.k1) then if (rang==0) then call sswap(1,p(indice),1,p(k1),1) end if call dswap(1,b(indice),1,b(k1),1) call dswap(nl,AL(1,indice),1,AL(1,k1),1) call MPI_BARRIER(MPI_COMM_WORLD,code) call position(k1,nl,knl,knlp) call position(indice,nl,inl,inlp) if (inl==knl) then if (rang==inl)then call dswap(n,AL(inlp,1),nl,AL(knlp,1),nl) end if else if (rang==inl) then call dcopy(n,AL(inlp,1),nl,temp,1) end if call MPI_BCAST(temp,n,MPI_DOUBLE_PRECISION,inl,MPI_COMM_WORLD,code) if (rang==knl)then call dcopy(n,AL(knlp,1),nl,temp2,1) end if call MPI_BCAST(temp2,n,MPI_DOUBLE_PRECISION,knl,MPI_COMM_WORLD,code) if (rang==inl)then call dcopy(n,temp2,1,AL(inlp,1),nl) end if if (rang==knl)then call dcopy(n,temp,1,AL(knlp,1),nl) end if end if end if !******************************************** ! End Permutation !******************************************** if (H_k1.ne.0.d0)then !******************************************** ! b=b/(b)_(i_0) !******************************************** D_inv=1.d0/H_k1 ! call dcopy(n,work,1,b,1) call dscal(n,D_inv,b,1) end if else H_k1=0.d0 call dscal(n,0.d0,b,1) end if !******************************************** ! Application of Givens rotations !******************************************** if (k.gt.1)then do i=2,k call position(i-1,nl,knl,knlp) call position(i,nl,inl,inlp) if (inl==knl)then if (rang==inl)then btemp=AL(knlp,k) AL(knlp,k)=c(i-1)*btemp+s(i-1)*AL(inlp,k) AL(inlp,k)=-s(i-1)*btemp+c(i-1)*AL(inlp,k) end if else if (rang==inl)then btemp=AL(inlp,k) end if if (rang==knl)then btemp2=AL(knlp,k) end if call MPI_BCAST(btemp,1,MPI_DOUBLE_PRECISION,inl,MPI_COMM_WORLD,code) call MPI_BCAST(btemp2,1,MPI_DOUBLE_PRECISION,knl,MPI_COMM_WORLD,code) net=net+2 if (rang==inl)then AL(inlp,k)=-s(i-1)*btemp2+c(i-1)*btemp end if if (rang==knl)then AL(knlp,k)=c(i-1)*btemp2+s(i-1)*btemp end if end if end do end if if (k.le.n)then call position(k,nl,knl,knlp) if (rang==knl)then btemp=AL(knlp,k) end if call MPI_BCAST(btemp,1,MPI_DOUBLE_PRECISION,knl,MPI_COMM_WORLD,code) dgam=dsqrt(btemp**2+H_k1**2) c(k)=dabs(btemp)/dgam s(k)=(H_k1*dabs(btemp))/(btemp*dgam) rs(k1)=-s(k)*rs(k) rs(k)=c(k)*rs(k) dnorm=dabs(rs(k1)) if (rang==knl)then AL(knlp,k)=c(k)*btemp+s(k)*H_k1 end if end if !******************************************** ! End Givens rotations !******************************************** end do !******************************************** ! Solve H_k d_k = alpha e_1 !******************************************** call position(k,nl,knl,knlp) do j=k,1,-1 call position(j,nl,inl,inlp) if (rang==inl) then rs(j)=rs(j)/AL(inlp,j) btemp=rs(j) end if call MPI_BCAST(btemp,1,MPI_DOUBLE_PRECISION,inl,MPI_COMM_WORLD,code) rs(j)=btemp work(j)=btemp call position(j-1,nl,inl,inlp) if (rang==inl)then call daxpy(inlp,-btemp,AL(1,j),1,rs(rang*nl+1),1) end if if (rang.lt.inl)then call daxpy(nl,-btemp,AL(1,j),1,rs(rang*nl+1),1) end if end do do j=k,1,-1 btemp=rs(j) call position(k,nl,knl,knlp) call position(j+1,nl,inl,inlp) if (rang==inl)then if (rang==knl)then call daxpy((knlp-inlp+1),btemp,AL(inlp,j),1,rs(rang*nl+inlp),1) end if if (rang.lt.knl)then call daxpy(nl-inlp+1,btemp,AL(inlp,j),1,rs(rang*nl+inlp),1) end if end if if (rang.gt.inl)then if (rang==knl)then call daxpy(knlp,btemp,AL(1,j),1,rs(rang*nl+1),1) end if if (rang.lt.knl)then call daxpy(nl,btemp,AL(1,j),1,rs(rang*nl+1),1) end if end if end do !******************************************** ! End Loop !******************************************** !******************************************** ! Update x_k=x_0 + L_k d_k !******************************************** call position(k1,nl,knl,knlp) if (rang.gt.knl)then call dgemv('N',nl,k,1.d0,AL(1,1),nl,work,1,0.d0,c(rang*nl+1),1) end if if (rang==knl)then call dgemv('N',nl-knlp+1,k,1.d0,AL(knlp,1),nl,work,1,0.d0,c(rang*nl+knlp),1) end if call position(k,nl,knl,knlp) if (rang.lt.knl)then call daxpy(nl,1.d0,rs(rang*nl+1),1,xl,1) end if if (rang==knl)then call daxpy(knlp,1.d0,rs(rang*nl+1),1,xl,1) call daxpy(nl-knlp,1.d0,c(rang*nl+knlp+1),1,xl(knlp+1),1) end if if (rang.gt.knl)then call daxpy(nl,1.d0,c(rang*nl+1),1,xl,1) end if !******************************************** ! Reorder the components of x_k !******************************************** call MPI_GATHER(xl,nl,MPI_DOUBLE_PRECISION,x,nl,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,code) call MPI_BARRIER(MPI_COMM_WORLD,code) if (rang==0)then do j=1,n b(p(j))=x(j) end do end if tempsf=(MPI_WTIME()-tempsf) tempsd=(MPI_WTIME()-tempsd) !write(*,*)'rang et temps',rang,tempsf,tempsd call MPI_REDUCE(tempsf,temps,1,MPI_DOUBLE_PRECISION,MPI_MAX,0,MPI_COMM_WORLD,code) call MPI_REDUCE(tempsd,temps2,1,MPI_DOUBLE_PRECISION,MPI_MAX,0,MPI_COMM_WORLD,code) if (rang==0)then if (iprint >= 2) then print*,'Results : iter, norm residual, time' write(*,40)k,dnorm,temps,temps2 end if write(30,10),k,dnorm close(30) open(unit=10,file='temps.txt',status='unknown',access='append') write(unit=10,FMT=*)'Proc. Number :',Nprocs,'Iter. Number:',k,'Time:',temps,temps2 close(10) end if if (rang==0)then if (iprint==2) then write(*,*) write(*,*)' la perm. et les solutions CMRH et EXACTE sont : ' write(*,*)' ______________________________________________ ' if (n.lt.10)then do i=1,n write(*,20)p(i), b(i), sol(i) end do else do i=1,10 write(*,20)p(i), b(i), sol(i) end do end if end if end if 30 format(A37,4(f14.7,1x)) 20 format(1x,i5,2(1x,f14.8)) 40 format(1x,i5,4(1x,e12.4)) ! deallocate(AL,b,sol,workl,temp,tempn,work,c,s,rs,xl,x,cl) if (rang==0) then deallocate(p) end if call MPI_FINALIZE(code) end program PCMRH !******************************************** ! End Parallel CMRH !******************************************** !******************************************** ! Matrix Choice !******************************************** subroutine matrice(A,nl,n,alpha,beta,choix,rang) integer ::n,i,j,choix,rang,nl,k real(kind=8),dimension(nl,n)::A real*8 ::alpha,beta if (choix==1) then do i = 1, nl k=rang*nl+i do j = 1, k A(i,j) = (2.d0*dble(j)-1.d0)/dble(n-k+j) end do do j = k+1, n A(i,j) = (2.d0*dble(k)-1.d0)/dble(n-k+j) end do end do end if if (choix==2)then do i=1,nl k=rang*nl+i do j=1,k A(i,j)=dble(k) end do do j=k+1,n A(i,j)=dble(j) end do end do end if end subroutine matrice !******************************************** ! Solution Choice !******************************************** subroutine solution(A,b,sol,nl,n,solchoice) integer ::n,i,nl,solchoice real(kind=8),dimension(n,n) ::A real(kind=8),dimension(n) ::sol real(kind=8),dimension(nl) ::b real ::rand do i=1,n if (solchoice==1) then sol(i)=1.d0 end if if (solchoice==2) then sol(i)=dble(n-i+1) end if if (solchoice==3)then sol(i)=rand(0) end if end do call dgemv('N',nl,n,1.d0,A,nl,sol,1,0.d0,b,1) end subroutine solution !******************************************** ! Where are we ? !******************************************** subroutine position(k,nl,knl,knlp) integer :: k,nl,knl,knlp knl=(k-1)/nl knlp=k-knl*nl end subroutine position