API tools faq
paste
Guest User

Untitled

a guest
May 12th, 2018
67
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 22.43 KB | None | 0 0
raw download clone embed print report
  1. program Parallel_Bratu_2D
  2. implicit none
  3. include 'mpif.h'
  4. !     2-dimensional non linear problem
  5. !     biquadratic finite element basis functions
  6. integer,parameter::nemax=5635,nnmax=23011,m=20
  7. integer*8:: nzeros
  8. integer:: nex=0,ney=0,nnx=0,nny=0,ne=0,np=0,nop(nemax,9),nell=0
  9. real:: xorigin=0,yorigin=0,xlast=0,ylast=0,deltax=0,deltay=0,xpt(nnmax),ypt(nnmax)
  10. real*8:: lamda=0,alfa=0,wtime_i=0,wtime_f=0,wtime=0
  11. real*8:: r1(nnmax),u(nnmax),d(nnmax),rest
  12. real*8,allocatable::AA(:)
  13. integer,allocatable:: JA(:)
  14. integer:: i,j,N=0,Nz=0,ncod(nnmax),p,q,iterat,dim=0,dim_buf=0
  15. integer:: k=0,l=0,indexAA(2),init_el=0,fin_el=0
  16. integer,allocatable:: IA(:),inner(:)
  17. real*8:: value=0,time=0
  18. !MPI variables
  19. integer:: nump,rank,ierr,tag1=0,tag2=0,right=0,left=0
  20. integer:: stats(MPI_status_size,8),reqs(8)
  21. real*8,allocatable:: buffer1(:),buffer2(:),buffer3(:),buffer4(:)
  22. call MPI_init(ierr)
  23. if (ierr .ne. MPI_success) then
  24.     write(*,*) 'Error starting MPI program'
  25.     call MPI_finalize(ierr)
  26. endif
  27. wtime_i=MPI_wtime()
  28. call MPI_comm_rank(MPI_comm_world,rank,ierr)
  29. call MPI_comm_size(MPI_comm_world,nump,ierr)
  30. !ORISMOS alfa,lamda
  31. alfa=0
  32. lamda=2
  33. !Mesh construction
  34. call xydiscr(nex,ney,xorigin,yorigin,xlast,ylast,deltax,deltay)
  35. call nodnumb(nemax,nnmax,ne,nnx,nny,nex,ney,np,nop)
  36. call xycoord(nemax,nnmax,xpt,ypt,xorigin,yorigin,nnx,nny,deltax,deltay)
  37. if (nex<nump) then
  38.     write(*,*) 'Condition nex>nump not met.Parallel solving not possible'
  39.     call MPI_finalize(ierr)
  40. endif
  41. !orismos N,m gia GMRES,dim_buf gia size_of buffer
  42. N=np
  43. nzeros=0.008*N*(N/nump)
  44. dim_buf=nny+2
  45. allocate(AA(nzeros))
  46. allocate(JA(nzeros))
  47. allocate(inner(np))
  48. allocate(buffer1(dim_buf))
  49. allocate(buffer2(dim_buf))
  50. allocate(buffer3(dim_buf))
  51. allocate(buffer4(dim_buf))
  52. !orismos diastasis IA gia CSR
  53. dim=np+1
  54. allocate(IA(dim))
  55. !Prepare for Dirichlet Boundary conditions
  56. do i=1,np
  57.     inner(i)=-1
  58.     ncod(i)=0
  59. enddo
  60. !markarisma aristera
  61. do i=1,nny
  62.     ncod(i)=1
  63. enddo
  64. !markarisma katw
  65. do i=1,(np-nny+1),nny
  66.     ncod(i)=1
  67. enddo
  68. !markarisma panw
  69. do i=nny,np,nny
  70.     ncod(i)=1
  71. enddo
  72. !markarisma dexia
  73. do i=(np-nny+1),np
  74.     ncod(i)=1
  75. enddo
  76.  
  77. iterat=1
  78. do while (iterat<=4)
  79.     iterat=iterat+1
  80. !arxikopoiisi
  81.     do i=1,dim
  82.         IA(i)=i
  83.     enddo
  84.     do i=1,nzeros
  85.         AA(i)=0
  86.         JA(i)=0
  87.     enddo
  88.     do i=1,np
  89.         r1(i)=0
  90.         d(i)=0
  91.     enddo
  92.     Nz=0
  93. !KSEKINIMA PARALLEL
  94. !kathe process sarwnei ena kommati tou mesh->sugekrimena elements
  95. if(mod(ne,nump)==0) then
  96.     init_el=(rank*int(ne/(nump))+1)
  97.     fin_el=(rank+1)*int(ne/nump)
  98.     do nell=init_el,fin_el
  99.              call abfind(nemax,nnmax,nell,r1,nop,xpt,ypt,u,lamda,alfa,value,AA,JA,IA,Nz,ncod,np,dim,nzeros)
  100.     enddo
  101. else
  102.     if(rank/=(nump-1)) then
  103.         init_el=(rank*int(ne/nump)+rank+1)
  104.         fin_el=((rank+1)*(int(ne/nump)+1))
  105.     do nell=init_el,((rank+1)*(int(ne/nump)+1))
  106.                  call abfind(nemax,nnmax,nell,r1,nop,xpt,ypt,u,lamda,alfa,value,AA,JA,IA,Nz,ncod,np,dim,nzeros)
  107.     enddo
  108.     else if(rank==(nump-1)) then
  109.         init_el=(rank*int(ne/nump)+rank+1)
  110.         fin_el=((rank+1)*(int(ne/nump))+mod(ne,nump))
  111.     do nell=init_el,fin_el
  112.                  call abfind(nemax,nnmax,nell,r1,nop,xpt,ypt,u,lamda,alfa,value,AA,JA,IA,Nz,ncod,np,dim,nzeros)
  113.     enddo
  114.     endif
  115. endif
  116. !indexAA periexei tin prwti kai teleutaia mi mideniki grammi gia ton AAlocal
  117. if(iterat==2) then
  118. i=0
  119. if (AA(1)/=0) then
  120.     indexAA(1)=1
  121.     else
  122.         do while (AA(IA(i))==0)
  123.             i=i+1
  124.         enddo
  125.         indexAA(1)=i
  126. endif
  127. do i=indexAA(1),np
  128.     if (AA(IA(i))/=0) then
  129.     indexAA(2)=i
  130.     endif
  131. enddo
  132. !Markarisma sunoriakwn nodes
  133. !INNER NODES->0
  134. !LEFT BOUNDARY NODES->1
  135. !RIGHT BOUNDARY NODES->2
  136. if(rank==0) then
  137.         do i=1,nop(fin_el,6)
  138.             inner(i)=0
  139.         enddo
  140.     else if(rank==nump-1) then
  141.         do i=nop(init_el,4),np
  142.             inner(i)=0
  143.         enddo
  144.     else
  145.         do i=nop(init_el,4),nop(fin_el,6)
  146.             inner(i)=0
  147.         enddo
  148.     endif
  149.    
  150.     if(rank/=0) then
  151.     !markarisma sunoriakwn aristera nodes,inner(boundary_node)=1
  152.         do i=0,2*(ney-mod((init_el-1),ney))
  153.             inner(nop(init_el,1)+i)=1
  154.         enddo
  155.         if(mod((init_el-1),ney)/=0) then
  156.             inner(nop(init_el,4))=1
  157.             do i=0,2*mod((init_el-1),ney)
  158.                 inner(nop(init_el+ney-mod((init_el-1),ney),1)+i)=1
  159.             enddo
  160.         endif
  161.     endif
  162.     if(rank/=(nump-1)) then
  163.     !markarisma sunoriakwn dexia nodes,inner(boundary_node)=2
  164.         if (mod(fin_el,ney)==0) then
  165.             do i=0,2*(ney)
  166.                 inner(nop(fin_el,9)-i)=2
  167.             enddo
  168.         endif
  169.         if(mod(fin_el,ney)/=0) then
  170.             do i=0,2*(mod(fin_el,ney))
  171.                 inner(nop(fin_el,9)-i)=2
  172.             enddo      
  173.             inner(nop(fin_el,6))=2
  174.             do i=0,2*(ney-mod(fin_el,ney))
  175.                 inner(nop(fin_el-mod(fin_el,ney),9)-i)=2
  176.             enddo
  177.         endif
  178. endif
  179. endif
  180. if(nump>1) then
  181. do i=1,dim_buf
  182.     buffer1=0
  183.     buffer2=0
  184.     buffer3=0
  185.     buffer4=0
  186. enddo
  187. tag1=10
  188. tag2=20
  189. left=rank-1
  190. right=rank+1
  191. !to kathe process exei to kommati tou r1 pou upologise+olokliro sta sunoriaka
  192. k=1
  193. l=1
  194. do i=indexAA(1),indexAA(2)
  195.     if (inner(i)==1) then
  196.         buffer2(k)=r1(i)
  197.     k=k+1
  198.     else if(inner(i)==2) then
  199.         buffer1(l)=r1(i)
  200.     l=l+1
  201.     endif
  202. enddo
  203. if (rank==0) then
  204.     call MPI_IRECV(buffer3,dim_buf,MPI_double_precision,right,tag1,MPI_COMM_WORLD,reqs(1),ierr)
  205.     call MPI_ISSEND(buffer1,dim_buf,MPI_double_precision,right,tag2,MPI_COMM_WORLD,reqs(2),ierr)
  206. else if(rank==nump-1) then
  207.     call MPI_IRECV(buffer4,dim_buf,MPI_double_precision,left,tag2,MPI_COMM_WORLD,reqs(1),ierr)
  208.     call MPI_ISSEND(buffer2,dim_buf,MPI_double_precision,left,tag1,MPI_COMM_WORLD,reqs(2),ierr)
  209. else
  210.     call MPI_IRECV(buffer4,dim_buf,MPI_double_precision,left,tag2,MPI_COMM_WORLD,reqs(1),ierr)
  211.     call MPI_ISSEND(buffer2,dim_buf,MPI_double_precision,left,tag1,MPI_COMM_WORLD,reqs(2),ierr)
  212.     call MPI_IRECV(buffer3,dim_buf,MPI_double_precision,right,tag1,MPI_COMM_WORLD,reqs(3),ierr)
  213.     call MPI_ISSEND(buffer1,dim_buf,MPI_double_precision,right,tag2,MPI_COMM_WORLD,reqs(4),ierr)
  214. endif  
  215. do i=1,2
  216.     call MPI_wait(reqs(i),stats(1,i),ierr)
  217. enddo
  218. if ((rank>0) .and. (rank<(nump-1))) then
  219.     do i=3,4
  220.         call MPI_wait(reqs(i),stats(1,i),ierr)
  221.     enddo
  222. endif
  223. k=1
  224. l=1
  225. do i=indexAA(1),indexAA(2)
  226.     if (inner(i)==1) then
  227.     r1(i)=r1(i)+buffer4(k)
  228.     k=k+1
  229.     else if(inner(i)==2) then
  230.     r1(i)=r1(i)+buffer3(l)
  231.     l=l+1
  232.     endif
  233. enddo
  234. endif
  235. !impose essential boundary conditions gia ton pinaka r1 / gia ton sk einai stin CSR mes
  236. do i=1,np
  237.     if (ncod(i).eq.1) then
  238.     r1(i)=0
  239.     r1(i)=-(u(i)-0)
  240.     endif
  241. enddo
  242. call GMRES_m(nemax,nnmax,r1,N,m,d,AA,JA,IA,dim,nzeros,indexAA,rank,nump,inner,reqs,stats,dim_buf)
  243. !KATHE PROCESS exei to kommati tis lusis tou+olo to sunoriako
  244. do i=indexAA(1),indexAA(2)
  245.     u(i)=u(i)+d(i)
  246. enddo
  247. enddo
  248. !   open(unit=1,file='./results1.txt',status='unknown')
  249. !   do i=1,np+1
  250. !       write(1,5618) xpt(i),ypt(i),u(i)
  251. !5618 format(3(e12.5,2x))
  252.  !    enddo
  253. !   close(1)
  254. wtime_f=MPI_wtime()
  255. wtime=wtime_f-wtime_i
  256.  
  257. if (rank==0) then
  258. write(*,*) '          ###############################################################'
  259. write(6,1050)
  260. 1050 format(/,15x,'Parallel solving of 2-D non linear Bratu Problem')
  261. write(*,1051)
  262. 1051 format(8x,'Solver: GMRES iterative   Matrix format: CSR   Parallel Library: MPI',/)
  263. write(*,231) rank,nump
  264. 231 format('Writing process:',3x,I2,3x,'out of',3x,I2,3x,'processes!',//)
  265. write(6,3002) nex,ney,ne,np
  266. 3002 format(5x,'nex=',i5,3x,'ney=',i5,3x,'ne=',i6,3x,'np=',i10,/)
  267. write(*,*) ''
  268. write(*,612)
  269. 612 format(/,'Problem Solved, Newton Converged',/)
  270. endif
  271. call MPI_barrier(mpi_comm_world,ierr)
  272. if (u(int(0.5*np+0.5))/=0) then
  273. write(*,*)  'u=',u(int(0.5*np+0.5)),'λ=',lamda,'rank=',rank
  274. endif
  275. call MPI_barrier(mpi_comm_world,ierr)
  276. call MPI_allreduce(wtime,time,1,MPI_double_precision,MPI_max,MPI_comm_world,ierr)
  277. if (rank==0) then
  278. write(*,222) time
  279. 222 format(/,'Maximum Time=',f30.5,'s',//)
  280. endif
  281. if(rank==0) then
  282.     write(*,*) 'NUMP=',nump,'nzeros=',nzeros
  283. endif
  284. !telos MPI world
  285. call MPI_finalize(ierr)
  286. end
  287.  
  288. subroutine xydiscr(nex,ney,xorigin,yorigin,xlast,ylast,deltax,deltay)
  289. implicit none
  290. !orismos arithmou elements-suntetagmenwn
  291. integer,intent(OUT):: nex,ney
  292. real,intent(OUT):: xorigin,yorigin,xlast,ylast,deltax,deltay
  293.       nex=25
  294.       ney=25
  295.       xorigin=0.
  296.       yorigin=0.
  297.       xlast=1
  298.       ylast=1
  299.       deltax=(xlast-xorigin)/nex
  300.       deltay=(ylast-yorigin)/ney
  301.       return
  302.       end
  303.  
  304.  
  305. subroutine nodnumb(nemax,nnmax,ne,nnx,nny,nex,ney,np,nop)
  306. implicit none
  307. !orismos arithmou nodes-node numbering
  308. integer,intent(IN)::nemax,nnmax,nex,ney
  309.       integer,intent(OUT)::ne,nnx,nny,np,nop(nemax,9)
  310.      integer:: nel=0,l=0,k=0,i=0,j=0
  311.       ne=nex*ney
  312.       nnx=2*nex+1
  313.       nny=2*ney+1
  314.       np=nnx*nny
  315. ! nodal numbering
  316.       nel=0
  317.       do 10 i=1,nex
  318.       do 10 j=1,ney
  319.       nel=nel+1
  320.       do 10 k=1,3
  321.       l=3*k-2
  322.       nop(nel,l)=nny*(2*i+k-3)+2*j-1
  323.       nop(nel,l+1)=nop(nel,l)+1
  324.  10   nop(nel,l+2)=nop(nel,l)+2
  325.       return
  326.       end
  327.  
  328.  
  329. subroutine xycoord(nemax,nnmax,xpt,ypt,xorigin,yorigin,nnx,nny,deltax,deltay)
  330. implicit none
  331. integer,intent(IN)::nemax,nnmax,nnx,nny
  332. real,intent(OUT):: xpt(nnmax),ypt(nnmax)
  333. real,intent(IN):: xorigin,yorigin,deltax,deltay
  334.   integer:: nnode=0,i=0,j=0
  335.       xpt(1)=xorigin
  336.       ypt(1)=yorigin
  337.       do 10 i=1,nnx
  338.       nnode=(i-1)*nny+1
  339.       xpt(nnode)=xpt(1)+(i-1)*deltax/2
  340.       ypt(nnode)=ypt(1)
  341.       do 20 j=2,nny
  342.       xpt(nnode+j-1)=xpt(nnode)
  343.       ypt(nnode+j-1)=ypt(nnode)+(j-1)*deltay/2
  344.  20   continue
  345.  10   continue
  346.       return
  347. end
  348.  
  349.  
  350. subroutine tsfun(x,y,phi,phic,phie)
  351.       real*8:: x,y
  352.       real*8,intent(OUT):: phi(9),phic(9),phie(9)
  353.       real*8:: l1,l2,l3,c
  354.       l1(c)=2.*c**2-3.*c+1.
  355.       l2(c)=-4.*c**2+4.*c
  356.       l3(c)=2.*c**2-c
  357.       dl1(c)=4.*c-3.
  358.       dl2(c)=-8.*c+4.
  359.       dl3(c)=4.*c-1.
  360.       phi(1)=l1(x)*l1(y)
  361.       phi(2)=l1(x)*l2(y)
  362.       phi(3)=l1(x)*l3(y)
  363.       phi(4)=l2(x)*l1(y)
  364.       phi(5)=l2(x)*l2(y)
  365.       phi(6)=l2(x)*l3(y)
  366.       phi(7)=l3(x)*l1(y)
  367.       phi(8)=l3(x)*l2(y)
  368.       phi(9)=l3(x)*l3(y)
  369.       phic(1)=l1(y)*dl1(x)
  370.       phic(2)=l2(y)*dl1(x)
  371.       phic(3)=l3(y)*dl1(x)
  372.       phic(4)=l1(y)*dl2(x)
  373.       phic(5)=l2(y)*dl2(x)
  374.       phic(6)=l3(y)*dl2(x)
  375.       phic(7)=l1(y)*dl3(x)
  376.       phic(8)=l2(y)*dl3(x)
  377.       phic(9)=l3(y)*dl3(x)
  378.       phie(1)=l1(x)*dl1(y)
  379.       phie(2)=l1(x)*dl2(y)
  380.       phie(3)=l1(x)*dl3(y)
  381.       phie(4)=l2(x)*dl1(y)
  382.       phie(5)=l2(x)*dl2(y)
  383.       phie(6)=l2(x)*dl3(y)
  384.       phie(7)=l3(x)*dl1(y)
  385.       phie(8)=l3(x)*dl2(y)
  386.       phie(9)=l3(x)*dl3(y)
  387.       return
  388.       end
  389.  
  390.  
  391.      subroutine abfind(nemax,nnmax,nell,r1,nop,xpt,ypt,u,lamda,alfa,value,AA,JA,IA,Nz,ncod,np,dim,nzeros)
  392. implicit none
  393.      integer*8,intent(IN):: nzeros
  394.      integer,intent(IN)::dim,nemax,nnmax
  395.      integer,intent(IN):: nop(nemax,9),ncod(nnmax)
  396.      real*8,intent(OUT)::r1(nnmax),value,AA(nzeros)
  397.      integer,intent(IN):: nell,np
  398.      integer,intent(INOUT)::JA(nzeros),IA(dim),Nz
  399.      real,intent(IN):: xpt(nnmax),ypt(nnmax)
  400.      real*8,intent(IN)::lamda,alfa,u(nnmax)
  401.      real*8:: phi(9),phic(9),phie(9),usol,usol1,usol2,x1,x2,y1,y2,dett
  402.      integer:: ngl(10)
  403.      real*8:: tphx(9),tphy(9)
  404.      real*8:: rest,w(3),gp(3),estifm(9,9)
  405.      integer i,j,k,p,q,t,r,in,l,m,n
  406.      data w,gp/0.27777777777778,0.444444444444,0.27777777777778,0.1127016654,0.5,0.8872983346/
  407.     do i=1,9
  408.         do j=1,9
  409.             estifm(i,j)=0
  410.         enddo
  411.     enddo
  412.       do 5 i=1,9
  413.  5    ngl(i)=nop(nell,i)
  414.       do 100 j=1,3
  415.       do 100 k=1,3
  416.       call tsfun(gp(j),gp(k),phi,phic,phie)
  417. ! isoparametric transformation
  418.       x1=0.
  419.       x2=0.
  420.       y1=0.
  421.       y2=0.
  422.     usol=0.
  423.     usol1=0.
  424.     usol2=0.
  425.       do 103 n=1,9
  426.       x1=x1+xpt(ngl(n))*phic(n)
  427.       x2=x2+xpt(ngl(n))*phie(n)
  428.       y1=y1+ypt(ngl(n))*phic(n)
  429.  103  y2=y2+ypt(ngl(n))*phie(n)
  430.       dett=x1*y2-x2*y1
  431.       do in=1,9
  432.       tphx(in)=(y2*phic(in)-y1*phie(in))/dett
  433.       tphy(in)=(x1*phie(in)-x2*phic(in))/dett
  434.     enddo
  435.  
  436. do n=1,9
  437. !athroisma sto oloklirwma
  438. usol=usol+u(ngl(n))*phi(n)
  439. usol1=usol1+u(ngl(n))*tphx(n)
  440. usol2=usol2+u(ngl(n))*tphy(n)
  441. enddo
  442. ! residuals
  443.       do 100 l=1,9
  444.       r1(ngl(l))=r1(ngl(l))+w(j)*w(k)*dett*(usol1*tphx(l)+usol2*tphy(l))
  445.       r1(ngl(l))=r1(ngl(l))-w(j)*w(k)*dett*lamda*phi(l)*exp(usol/(1+alfa*usol))
  446.       do 100 m=1,9
  447.       estifm(l,m)=estifm(l,m)-w(j)*w(k)*dett*((tphx(l)*tphx(m))+(tphy(l)*tphy(m)))
  448.       estifm(l,m)=estifm(l,m)+w(j)*w(k)*dett*lamda*phi(l)*phi(m)*exp(usol/(1+(alfa*usol)))*(1/(1+alfa*usol)**2)
  449.       100  continue
  450.       do t=1,9
  451.         do r=1,9
  452.             value=estifm(t,r)
  453.             p=ngl(t)
  454.             q=ngl(r)
  455.         call CSR(nemax,nnmax,AA,JA,IA,p,q,value,Nz,ncod,np,dim,nzeros)
  456.         enddo
  457.       enddo
  458. end
  459.  
  460.  
  461.  
  462. subroutine GMRES_m(nemax,nnmax,r1,N,m,d,AA,JA,IA,dim,nzeros,indexAA,rank,nump,inner,reqs,stats,dim_buf)
  463. implicit none
  464. include 'mpif.h'
  465. integer*8,intent(IN):: nzeros
  466. integer,intent(IN)::N,m,dim,rank,nump,dim_buf,nemax,nnmax
  467. integer,intent(INOUT)::reqs(8),stats(MPI_status_size,8)
  468. integer,intent(IN)::IA(dim),JA(nzeros),indexAA(2),inner(dim)
  469. real*8,intent(IN)::r1(nnmax),AA(nzeros)
  470. real*8,intent(OUT)::d(nnmax)
  471. integer:: i=0,j=0,k=0,l=0,k1=0,k2=0,iter=0,tag1=0,tag2=0,right=0,left=0,ierr
  472. real*8:: tot_rest=0,vita=0,eps=0,buffer1(dim_buf),buffer2(dim_buf),buffer3(dim_buf),buffer4(dim_buf),buff=0
  473. real*8:: b(N),Hm(m+1,m),x(N),x0(N),r0(N),e(m+1),y(m),u_base(N,m),uj(N),g(m+1),rest(N)
  474. !arxikopoiisi pinakwn kai eisagwgi dedomenwn
  475. do i=indexAA(1),indexAA(2)
  476.     b(i)=r1(i)
  477. enddo
  478. do i=1,N
  479.     x(i)=0
  480. enddo
  481. tot_rest=1
  482. iter=0
  483. !epithumiti akriveia
  484. iter=1
  485. do while(iter<=30)
  486. do i=indexAA(1),indexAA(2)
  487.     x0(i)=x(i)
  488. enddo
  489. iter=iter+1
  490. !midenismos pinakwn-anusmatwn Krylov
  491. do i=1,N
  492.     uj(i)=0
  493.     r0(i)=0
  494.     do j=1,m
  495.         u_base(i,j)=0
  496.     enddo
  497. enddo
  498. do i=1,m+1
  499.     e(i)=0
  500.     g(i)=0
  501.     do j=1,m
  502.         Hm(i,j)=0
  503.         y(j)=0
  504.     enddo
  505. enddo
  506. e(1)=1
  507. !upologismos dianusmatos r0=b-Ax0 me r0=b-MATMUL(A,x0) se CSR format
  508. do i=indexAA(1),indexAA(2)
  509.     k1=IA(i)
  510.     k2=IA(i+1)-1
  511.     do j=k1,k2
  512.         r0(i)=r0(i)+AA(j)*x0(JA(j))
  513.     enddo
  514. enddo
  515. !to kathe process exei to kommati tou r1 pou upologise+olokliro sta sunoriaka
  516. if(nump>1) then
  517. do i=1,dim_buf
  518.     buffer1(i)=0
  519.     buffer2(i)=0
  520.     buffer3(i)=0
  521.     buffer4(i)=0
  522. enddo
  523. tag1=11
  524. tag2=21
  525. right=rank+1
  526. left=rank-1
  527. k=1
  528. l=1
  529. do i=indexAA(1),indexAA(2)
  530.     if (inner(i)==1) then
  531.         buffer2(k)=r0(i)
  532.     k=k+1
  533.     else if(inner(i)==2) then
  534.         buffer1(l)=r0(i)
  535.     l=l+1
  536.     endif
  537. enddo
  538. if (rank==0) then
  539.     call MPI_IRECV(buffer3,dim_buf,MPI_double_precision,right,tag1,MPI_COMM_WORLD,reqs(1),ierr)
  540.     call MPI_ISSEND(buffer1,dim_buf,MPI_double_precision,right,tag2,MPI_COMM_WORLD,reqs(2),ierr)
  541. else if(rank==nump-1) then
  542.     call MPI_IRECV(buffer4,dim_buf,MPI_double_precision,left,tag2,MPI_COMM_WORLD,reqs(1),ierr)
  543.     call MPI_ISSEND(buffer2,dim_buf,MPI_double_precision,left,tag1,MPI_COMM_WORLD,reqs(2),ierr)
  544. else
  545.     call MPI_IRECV(buffer4,dim_buf,MPI_double_precision,left,tag2,MPI_COMM_WORLD,reqs(1),ierr)
  546.     call MPI_ISSEND(buffer2,dim_buf,MPI_double_precision,left,tag1,MPI_COMM_WORLD,reqs(2),ierr)
  547.     call MPI_IRECV(buffer3,dim_buf,MPI_double_precision,right,tag1,MPI_COMM_WORLD,reqs(3),ierr)
  548.     call MPI_ISSEND(buffer1,dim_buf,MPI_double_precision,right,tag2,MPI_COMM_WORLD,reqs(4),ierr)
  549. endif  
  550. do i=1,2
  551.     call MPI_wait(reqs(i),stats(1,i),ierr)
  552. enddo
  553. if ((rank>0) .and. (rank<(nump-1))) then
  554.     do i=3,4
  555.         call MPI_wait(reqs(i),stats(1,i),ierr)
  556.     enddo
  557. endif
  558. k=1
  559. l=1
  560. do i=indexAA(1),indexAA(2)
  561.     if (inner(i)==1) then
  562.         r0(i)=r0(i)+buffer4(k)
  563.     k=k+1
  564.     else if(inner(i)==2) then
  565.         r0(i)=r0(i)+buffer3(l)
  566.     l=l+1
  567.     endif
  568. enddo
  569. endif
  570. do i=1,N
  571. r0(i)=-r0(i)+b(i)
  572. enddo
  573. !end CSR multi
  574. vita=0
  575. !ypologismos normas tou r0 k apothikeusi sto vita
  576. !prosoxi na min vazei 2 fores ta sunoriaka
  577. do i=indexAA(1),indexAA(2)
  578.     if((inner(i)==0) .or. (inner(i)==1)) then
  579.         vita=vita+r0(i)**2
  580.     endif
  581. enddo
  582. buff=0
  583. call MPI_allreduce(vita,buff,1,MPI_double_precision,MPI_sum,MPI_comm_world,ierr)
  584. vita=sqrt(buff)
  585. !upologismos arxikou dianusmatos uj k apothikeusi ston uj kai u_base(i,1)
  586. uj=r0/vita
  587. !apothikeusi tou u1 ston pinaka vasewn tou upoxwrou Krylov(to kathena apothikeuei to kommati tou)
  588. do i=1,N
  589.     u_base(i,1)=uj(i)
  590. enddo
  591. !dianusma g=b*e1
  592. g=vita*e
  593. !sunartisi pou paragei to xwro Krylov
  594. call Krylov(AA,JA,IA,dim,nzeros,Hm,u_base,uj,N,m,indexAA,rank,nump,inner,reqs,stats,dim_buf)
  595.  
  596. call Least_sq(Hm,g,m)
  597.  
  598. call Back_sub(Hm,y,g,m)
  599.  
  600. x=x0+MATMUL(u_base,y)
  601. enddo
  602. do i=indexAA(1),indexAA(2)
  603.     d(i)=x(i)
  604. enddo
  605. end subroutine GMRES_m
  606.  
  607. subroutine Krylov(AA,JA,IA,dim,nzeros,Hm,u_base,uj,N,m,indexAA,rank,nump,inner,reqs,stats,dim_buf)
  608. implicit none
  609. include 'mpif.h'
  610. integer:: i=0,j=0,k=0,k1=0,k2=0,p=0,tag1=0,tag2=0,right=0,left=0,count=0,ierr,l=0
  611. integer*8,intent(IN):: nzeros
  612. integer,intent(IN):: N,m,dim,rank,nump,dim_buf
  613. integer,intent(INOUT)::reqs(8),stats(MPI_status_size,8)
  614. integer,intent(IN):: JA(nzeros),IA(dim),indexAA(2),inner(dim)
  615. real*8,intent(OUT)::Hm(m+1,m),u_base(N,m),uj(N)
  616. real*8,intent(IN):: AA(nzeros)
  617. real*8:: w(N),buffer1(dim_buf),buffer2(dim_buf),buffer3(dim_buf),buffer4(dim_buf),buff,buff1
  618. do j=1,m
  619. !to j einai count gia mas!
  620. !vazei to kommati tou uj pou antistoixei sto process       
  621.     do k=indexAA(1),indexAA(2)
  622.     if(inner(k)/=(-1)) then
  623.         uj(k)=u_base(k,j)
  624.     endif
  625.     enddo
  626.  
  627. !MATMUL me CSR format w=MATMUL(A,uj)
  628. do i=1,N
  629.     w(i)=0
  630. enddo
  631. do i=indexAA(1),indexAA(2)
  632. if(inner(i)/=(-1)) then
  633.     k1=IA(i)
  634.     k2=IA(i+1)-1
  635.     do k=k1,k2
  636.         w(i)=w(i)+AA(k)*uj(JA(k))
  637.     enddo
  638. endif
  639. enddo
  640. if(nump>1) then
  641. do i=1,dim_buf
  642.     buffer1(i)=0
  643.     buffer2(i)=0
  644.     buffer3(i)=0
  645.     buffer4(i)=0
  646. enddo
  647. tag1=31
  648. tag2=41
  649. right=rank+1
  650. left=rank-1
  651. k=1
  652. l=1
  653. do i=indexAA(1),indexAA(2)
  654.     if (inner(i)==1) then
  655.         buffer2(k)=w(i)
  656.     k=k+1
  657.     else if(inner(i)==2) then
  658.         buffer1(l)=w(i)
  659.     l=l+1
  660.     endif
  661. enddo
  662. if (rank==0) then
  663.     call MPI_IRECV(buffer3,dim_buf,MPI_double_precision,right,tag1,MPI_COMM_WORLD,reqs(1),ierr)
  664.     call MPI_ISSEND(buffer1,dim_buf,MPI_double_precision,right,tag2,MPI_COMM_WORLD,reqs(2),ierr)
  665. else if(rank==nump-1) then
  666.     call MPI_IRECV(buffer4,dim_buf,MPI_double_precision,left,tag2,MPI_COMM_WORLD,reqs(1),ierr)
  667.     call MPI_ISSEND(buffer2,dim_buf,MPI_double_precision,left,tag1,MPI_COMM_WORLD,reqs(2),ierr)
  668. else
  669.     call MPI_IRECV(buffer4,dim_buf,MPI_double_precision,left,tag2,MPI_COMM_WORLD,reqs(1),ierr)
  670.     call MPI_ISSEND(buffer2,dim_buf,MPI_double_precision,left,tag1,MPI_COMM_WORLD,reqs(2),ierr)
  671.     call MPI_IRECV(buffer3,dim_buf,MPI_double_precision,right,tag1,MPI_COMM_WORLD,reqs(3),ierr)
  672.     call MPI_ISSEND(buffer1,dim_buf,MPI_double_precision,right,tag2,MPI_COMM_WORLD,reqs(4),ierr)
  673. endif  
  674. do i=1,2
  675.     call MPI_wait(reqs(i),stats(1,i),ierr)
  676. enddo
  677. if ((rank>0) .and. (rank<(nump-1))) then
  678.     do i=3,4
  679.         call MPI_wait(reqs(i),stats(1,i),ierr)
  680.     enddo
  681. endif
  682. k=1
  683. l=1
  684. do i=indexAA(1),indexAA(2)
  685.     if (inner(i)==1) then
  686.         w(i)=w(i)+buffer4(k)
  687.     k=k+1
  688.     else if(inner(i)==2) then
  689.         w(i)=w(i)+buffer3(l)
  690.     l=l+1
  691.     endif
  692. enddo
  693. endif
  694. !end CSR multi
  695. do i=1,j
  696.         do k=indexAA(1),indexAA(2)
  697.         if(inner(k)/=(-1)) then
  698.             uj(k)=u_base(k,i)
  699.         endif
  700.         enddo
  701. !       Hm(i,j)=DOT_PRODUCT(w,uj)
  702.     Hm(i,j)=0
  703.     do k=indexAA(1),indexAA(2)
  704. !einai swsto to IF ALLA THA TO VGALW META ti suglisi wste na sugrinw me to test
  705.         if((inner(k)==0) .or. (inner(k)==1)) then
  706.             Hm(i,j)=Hm(i,j)+w(k)*uj(k)
  707.         endif
  708.     enddo  
  709. buff=0
  710. buff1=Hm(i,j)
  711. call MPI_allreduce(buff1,buff,1,MPI_double_precision,MPI_sum,MPI_comm_world,ierr)
  712. Hm(i,j)=buff
  713. !END DOT PRODUCT
  714. do k=indexAA(1),indexAA(2)
  715.     if(inner(k)/=(-1)) then
  716.         w(k)=w(k)-Hm(i,j)*uj(k)
  717.     endif
  718. enddo
  719. enddo
  720.     do k=indexAA(1),indexAA(2)
  721.         if((inner(k)==0) .or. (inner(k)==1)) then
  722.             Hm(j+1,j)=Hm(j+1,j)+(w(k))**2
  723.         endif
  724.     enddo
  725. buff=0
  726. call MPI_allreduce(Hm(j+1,j),buff,1,MPI_double_precision,MPI_sum,MPI_comm_world,ierr)  
  727.     Hm(j+1,j)=sqrt(buff)
  728. if(j<m) then
  729.         do k=indexAA(1),indexAA(2)
  730.             if(inner(k)/=(-1)) then
  731.                 u_base(k,j+1)=w(k)/Hm(j+1,j)
  732.             endif
  733.         enddo
  734.     endif
  735. enddo
  736. end subroutine Krylov  
  737.    
  738. subroutine Least_sq(Hm,g,m)
  739. implicit none
  740. integer:: i,j,k
  741. integer,intent(IN):: m
  742. real*8,intent(OUT)::Hm(m+1,m),g(m+1)
  743. real*8:: W(m+1,m+1),si,ci
  744. !m pollaplasiasmoi pinakwn
  745. do i=1,m
  746. !upologismos si kai ci
  747.     si=Hm(i+1,i)/sqrt((Hm(i,i)**2)+(Hm(i+1,i)**2))
  748.     ci=Hm(i,i)/sqrt((Hm(i,i)**2)+(Hm(i+1,i)**2))
  749. !kataskeui pinaka W
  750.     do j=1,m+1
  751.         do k=1,m+1
  752.             if (j==k) then
  753.                 W(j,k)=1
  754.                 else
  755.                 W(j,k)=0
  756.             endif
  757.         enddo
  758.     enddo
  759. !eisagwgi timwn pinaka stin i,i+1 grammi
  760.     W(i,i)=ci
  761.     W(i,i+1)=si
  762.     W(i+1,i)=-si
  763.     W(i+1,i+1)=ci
  764. Hm=MATMUL(W,Hm)
  765. g=MATMUL(W,g)
  766. enddo
  767. end subroutine Least_sq            
  768.  
  769. subroutine Back_sub(Hm,y,g,m)
  770. implicit none
  771. integer:: i,j,k
  772. integer,intent(IN):: m
  773. real*8,intent(IN)::Hm(m+1,m),g(m+1)
  774. real*8,intent(OUT):: y(m+1)
  775. real*8:: sum=0,help(m+1)
  776. sum=0
  777. do i=1,m+1
  778.     help(i)=g(i)
  779. enddo
  780. do i=m,1,-1
  781.     sum=help(i)
  782.     if (i.LT.m) then
  783.         do j=i+1,m
  784.             sum=sum-Hm(i,j)*help(j)
  785.         enddo
  786.     endif
  787.     help(i)=sum/Hm(i,i)
  788. enddo
  789. do i=1,m
  790.     y(i)=help(i)
  791. enddo
  792. end subroutine Back_sub
  793.    
  794. subroutine CSR(nemax,nnmax,AA,JA,IA,i,j,value,Nz,ncod,np,dim,nzeros)
  795. implicit none
  796.  integer*8,intent(IN):: nzeros
  797.  integer,intent(IN):: i,j,ncod(nnmax),np,dim,nemax,nnmax
  798.  real*8,intent(INOUT):: AA(nzeros),value
  799.  real*8:: BB(nzeros)
  800.  integer,intent(INOUT)::JA(nzeros),IA(dim),Nz
  801. integer:: k=0,l=0,t=0,update=0
  802. !edw ksekiname
  803. update=0
  804. !Dirichlet Boundary Condition
  805. if (j/=np+1) then
  806. if ((ncod(i)==1) .and. (i/=j)) then
  807.     value=0
  808.     else if ((ncod(i)==1) .and. (i==j) .and. (AA(IA(i))==0)) then
  809.     value=1
  810.     else if ((ncod(i)==1) .and. (i==j) .and. (AA(IA(i))/=0)) then
  811.     value=0
  812. endif
  813. endif
  814. !Ksekinaei i eisagwgi
  815. if (abs(value)>1E-10) then
  816. !to stoixeio IA(i) deixnei ti thesi ston AA pou einai to prwto mi mideniko tou pinaka i
  817.                 if(AA(IA(i))==0) then
  818.                     AA(IA(i))=value
  819.                     JA(IA(i))=j
  820.                     Nz=Nz+1
  821.                 else
  822.                 t=IA(i)
  823.                 do while (t<IA(i+1))
  824.                     if(JA(t)==j) then
  825.                         update=1
  826.                         exit
  827.                     endif
  828.                     t=t+1
  829.                 enddo
  830.                 endif  
  831.                 if (update==1 .and. ncod(i)/=1) then
  832.                     AA(t)=AA(t)+value
  833.                     else if (update==1 .and. (j==np+1)) then
  834.                         AA(t)=AA(t)+value
  835.                     else
  836.                     if ((AA(IA(i))/=0) .and. (JA(IA(i))>j)) then
  837.                     do k=(IA(dim)+2),IA(i),-1
  838.                         AA(k+1)=AA(k)
  839.                         JA(k+1)=JA(k)
  840.                     enddo
  841.                     AA(IA(i))=value
  842.                     JA(IA(i))=j
  843.                     do k=i+1,dim
  844.                         IA(k)=IA(k)+1
  845.                     enddo
  846.                     Nz=Nz+1
  847.                     else if  ((AA(IA(i))/=0) .and. (JA(IA(i))<j)) then
  848.                     k=IA(i)
  849.                     do while ((JA(k)<j) .and. (k<IA(i+1)))
  850.                         k=k+1
  851.                     enddo
  852.                         do l=(IA(dim)+2),k,-1
  853.                             AA(l+1)=AA(l)
  854.                             JA(l+1)=JA(l)
  855.                         enddo
  856.                         AA(k)=value
  857.                         JA(k)=j
  858.                         do k=i+1,dim
  859.                             IA(k)=IA(k)+1
  860.                         enddo
  861.                         Nz=Nz+1
  862.                 endif
  863.             endif
  864. endif
  865.     l=IA(dim-1)+1
  866.     do while (AA(l)/=0)
  867.     l=l+1
  868.     enddo
  869.     IA(dim)=l
  870.     end subroutine CSR
  871.    
  872. real function get_CSR(AA,JA,IA,dim,nzeros,i,j)
  873. implicit none
  874. integer,intent(IN):: dim,nzeros,i,j
  875. integer,intent(IN):: JA(nzeros),IA(dim)
  876. real*8,intent(IN):: AA(nzeros)
  877. integer:: k=0,l=0,test=0
  878. !Default value
  879. get_CSR=0
  880. k=IA(i)
  881. if (JA(IA(i))==j) then
  882.     get_CSR=AA(IA(i))
  883.     else
  884.         do k=IA(i),(IA(i+1)-1)
  885.             if (JA(k)==j) then
  886.                 get_CSR=AA(k)
  887.                 exit
  888.             else
  889.                 get_CSR=0
  890.             endif
  891.         enddo
  892. endif
  893. end function get_CSR
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!