DOUG 0.2

Vect.F90

Go to the documentation of this file.
00001 ! DOUG - Domain decomposition On Unstructured Grids
00002 ! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
00003 ! Department of Mathematics, University of Bath
00004 !
00005 ! This library is free software; you can redistribute it and/or
00006 ! modify it under the terms of the GNU Lesser General Public
00007 ! License as published by the Free Software Foundation; either
00008 ! version 2.1 of the License, or (at your option) any later version.
00009 !
00010 ! This library is distributed in the hope that it will be useful,
00011 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013 ! Lesser General Public License for more details.
00014 !
00015 ! You should have received a copy of the GNU Lesser General Public
00016 ! License along with this library; if not, write to the Free Software
00017 ! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00018 ! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
00019 ! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
00020 ! mailto:info(at)dougdevel.org)
00021 
00022 !!----------------------------------------
00023 !! Useful subroutines to work with vectors
00024 !!----------------------------------------
00025 module Vect_mod
00026 
00027   use DOUG_utils
00028   use RealKind
00029   use Mesh_class
00030   use DenseMtx_mod
00031 
00032   implicit none
00033 
00034 #include<doug_config.h>
00035 
00036 #ifdef D_COMPLEX
00037 #define float complex
00038 #else
00039 #define float real
00040 #endif
00041 
00042   ! = = = = = = = = = =
00043   ! Indexes of freedoms to be exchanged between neighbours :
00044   ! fexchindxV[maxval(Mesh%nfreesend_map), Mesh%nnghbrs]
00045   integer, dimension(:,:), pointer :: fexchindxV
00046   ! Auxiliary array for indexing 'fexchindxV(:,pid2indx(pid))'
00047   integer,   dimension(:), pointer :: pid2indxV
00048 
00049   type farrayV
00050      float(kind=rk), dimension(:), pointer :: arr
00051   end type farrayV
00052   ! Input/output bufers for send/receiv freedoms
00053   type(farrayV), dimension(:), pointer :: inbufsV
00054   type(farrayV), dimension(:), pointer :: outbufsV
00055 
00056   ! Whether auxiliary arrays for assisting in
00057   ! communications initialised or not
00058   logical :: D_COMMSTRUCTS_INITEDV = .false.
00059 
00060   ! = = = = = = = = = =
00061   ! Points to the last element where interface freedoms
00062   ! end in a permuted (oldToNew) RHS or solution vector
00063   !    inerf   v        inner
00064   ! /[..........]..................../
00065   integer                        :: intf_end = -1
00066   ! Mask to perform dot product in case of copied values
00067   ! of interface freedoms : dot_intf_fmask[intf_end]
00068   integer, dimension(:), pointer :: dot_intf_fmask
00069   ! Map : dot_intf_fmap[<=intf_end]
00070   integer, dimension(:), pointer :: dot_intf_fmap
00071   ! Aux. parameters
00072   integer, parameter :: D_DOTMASK_MINE    = 1
00073   integer, parameter :: D_DOTMASK_NOTMINE = 0
00074   integer :: ninner
00075   integer, parameter :: D_RHS_TEXT   = 0
00076   integer, parameter :: D_RHS_BINARY = 1
00077   integer, parameter :: D_RHS_XDR    = 2
00078 
00079   ! = = = = = = = = = =
00080   interface Vect_Print
00081      module procedure DVect_Print, IVect_Print, I1Vect_Print;
00082   end interface
00083 
00084   ! = = = = = = = = = =
00085   private :: &
00086        CommStructs_init,    &
00087        CommStructs_destroy, &
00088        DVect_Print,         &
00089        IVect_Print,         &
00090        I1Vect_Print
00091 
00092 contains
00093 
00094 
00095   !--------------------------
00096   ! Deallocate allocated data
00097   !--------------------------
00098   subroutine Vect_cleanUp()
00099     implicit none
00100 
00101     if (associated(fexchindxV)) deallocate(fexchindxV)
00102     if (associated(pid2indxV)) deallocate(pid2indxV)
00103     if (associated(inbufsV)) deallocate(inbufsV)
00104     if (associated(outbufsV)) deallocate(outbufsV)
00105     if (associated(dot_intf_fmask)) deallocate(dot_intf_fmask)
00106     if (associated(dot_intf_fmap)) deallocate(dot_intf_fmap)
00107   end subroutine Vect_cleanUp
00108 
00109   ! subroutine Vect_permute(x, perm)
00110   ! subroutine Vect_assembleFromElem(x, E, M, perm)
00111 
00112   ! ===============================================
00113 
00114   !----------------------------------------------
00115   ! Sets value to this module's local 'intf_end'
00116   !----------------------------------------------
00117   subroutine Vect_setIntfEnd(end)
00118     implicit none
00119     integer, intent(in) :: end
00120     intf_end = end
00121   end subroutine Vect_setIntfEnd
00122 
00123 
00124   !--------------------------------------------------------------
00125   ! Build 'dot_intf_fmask' mask to assist in parallel dot product
00126   ! where interface freedoms are copied among all processors
00127   !--------------------------------------------------------------
00128   subroutine Vect_buildDotMask(M)
00129     implicit none
00130 
00131     type(Mesh),            intent(in) :: M
00132 
00133     integer :: lf, gf, h, ge, s
00134 
00135 ! Demonstrates unoptimality of the hash construct:
00136 !integer :: i
00137 !write(stream,*)' hash: #####################################################'
00138 !do i =1,M%hashsize
00139 !! if (M%hash(i,1)>0) write(stream,*)i,' hash:',M%hash(i,1),M%hash(i,2)
00140 !enddo
00141 !write(stream,*)' .....hash: ################################################'
00142 !write(stream,*)'M%hscale=',M%hscale
00143 
00144     if (sctls%input_type==DISTRIBUTION_TYPE_ASSEMBLED.or.&
00145          sctls%input_type==DISTRIBUTION_TYPE_STRUCTURED) then
00146       ninner=M%ninner
00147       return
00148     endif
00149     if (intf_end == -1) then
00150       s=sum(M%inner_interf_fmask)
00151       call Vect_setIntfEnd(s)
00152     endif
00153     allocate(dot_intf_fmask(intf_end))
00154     ! Mark all freedoms as to perform dot product on
00155     dot_intf_fmask = D_DOTMASK_MINE
00156 
00157     do lf = 1,M%nlf ! cycle through local freedoms
00158        ! interface freedoms
00159        if (M%inner_interf_fmask(lf) == D_FREEDOM_INTERF) then
00160           gf = M%lg_fmap(lf) ! map index of freedom to global numbering
00161           h = M%hashlook(int(gf/M%hscale)+1)
00162           do while (M%hash(h,1) > 0)
00163              if (M%hash(h,1) == gf) then
00164                 ge = M%hash(h,2)
00165                 ! Don't operate on those whos rank is lower than mine
00166                 if (M%eptnmap(ge) < myrank+1) then
00167                    !dot_intf_fmask(newToOldPerm(lf)) = D_DOTMASK_NOTMINE
00168                    dot_intf_fmask(lf) = D_DOTMASK_NOTMINE
00169                 end if
00170              end if
00171              h = h + 1
00172           end do ! do while
00173        end if
00174     end do
00175   end subroutine Vect_buildDotMask
00176 
00177 
00178   !------------------------------------
00179   ! Builds map to assist in dot product
00180   !------------------------------------
00181   subroutine Vect_buildDotMap()
00182     implicit none
00183 
00184     integer :: i, j
00185 
00186     if (sctls%input_type==DISTRIBUTION_TYPE_ASSEMBLED) then
00187       return
00188     endif
00189     if (.not.associated(dot_intf_fmask)) &
00190          call DOUG_abort('[Vect_buildDotMap] :SEVERE: dot_intf_fmask must'//&
00191          ' be built first.',-1)
00192     allocate(dot_intf_fmap(sum(dot_intf_fmask)))
00193     j = 1
00194     do i = 1,intf_end
00195        if (dot_intf_fmask(i) == D_DOTMASK_MINE) then
00196             dot_intf_fmap(j) = i
00197             j = j + 1
00198          end if
00199     end do
00200   end subroutine Vect_buildDotMap
00201 
00202 
00203   !--------------------------------------------
00204   ! Global dot product
00205   !--------------------------------------------
00206   function Vect_dot_product(x1, x2) result(res)
00207     implicit none
00208 
00209     float(kind=rk), dimension(:), intent(in) :: x1, x2
00210     float(kind=rk) :: res
00211 
00212     float(kind=rk) :: local_intf, local_res
00213     integer        :: ierr,i
00214 
00215     if (size(x1) /= size(x2)) &
00216          call DOUG_abort('[Vect_dot_product] :SEVERE: sizes of vectors'//&
00217          ' must match.',-1)
00218 
00219     if (numprocs == 1) then
00220        res = dot_product(x1, x2)
00221        return
00222     end if
00223     
00224     if (sctls%input_type==DISTRIBUTION_TYPE_ASSEMBLED.or.&
00225          sctls%input_type==DISTRIBUTION_TYPE_STRUCTURED) then
00226       local_res=dot_product(x1(1:ninner),x2(1:ninner))
00227     else
00228       ! Interface freedoms defined in 'dot_intf_fmap' + local inner ones
00229       local_res = dot_product(x1(dot_intf_fmap), x2(dot_intf_fmap)) + &
00230            dot_product(x1(intf_end+1:), x2(intf_end+1:))
00231       !local_res=0.0_rk
00232       !do i=1,size(dot_intf_fmap)
00233       !  local_res=local_res+x1(dot_intf_fmap(i))*x2(dot_intf_fmap(i))
00234       !enddo
00235       !local_res = local_res + &
00236       !     dot_product(x1(intf_end+1:), x2(intf_end+1:))
00237     endif
00238     call MPI_ALLREDUCE(local_res, res, 1, MPI_fkind, MPI_SUM, &
00239          MPI_COMM_WORLD, ierr)
00240   end function Vect_dot_product
00241 
00242 
00243   !-------------------------------
00244   ! Permute vector
00245   !-------------------------------
00246   subroutine Vect_permute(x, perm)
00247     implicit none
00248 
00249     float(kind=rk), dimension(:), intent(in out) :: x
00250     integer,        dimension(:), intent(in)     :: perm
00251 
00252     integer                               :: n, i
00253     float(kind=rk), dimension(:), pointer :: xtmp
00254 
00255     n = size(x)
00256     if (n /= size(perm)) &
00257          call DOUG_abort('[Vect_permute] : SEVERE : size(x) /= size(perm)',-1)
00258 
00259     allocate(xtmp(n))
00260     do i = 1,n
00261        xtmp(i) = x(perm(i))
00262     end do
00263     x = xtmp
00264 
00265     deallocate(xtmp)
00266   end subroutine Vect_permute
00267 
00268   !--------------------------------------------
00272   subroutine Vect_Gather(xl, x, M, rank)
00273     implicit none
00274 
00275     float(kind=rk), dimension(:),     intent(in) :: xl ! local vector
00276     float(kind=rk), dimension(:), intent(in out) :: x  ! global vector
00277     type(Mesh),                       intent(in) :: M
00278     ! Rank of a process to gather vector on
00279     integer,                optional, intent(in) :: rank
00280 
00281     ! Array of number of local freedoms
00282     integer, dimension(:), pointer :: nlfs
00283     logical, dimension(:), pointer :: f_booked
00284     integer :: i, gf, h, p
00285 
00286     ! MPI
00287     ! Aux. receive buffers on master side
00288     integer,        dimension(:), pointer :: ibuf
00289     float(kind=rk), dimension(:), pointer :: fbuf
00290     integer, parameter :: D_TAG_NLF   = 222
00291     integer, parameter :: D_TAG_FMAP  = 333
00292     integer, parameter :: D_TAG_FDATA = 444
00293     integer :: bs, ierr, request, status(MPI_STATUS_SIZE)
00294     integer :: root
00295 
00296     if (numprocs==1) then
00297       x=xl
00298       return
00299     endif
00300     if (present(rank)) then
00301        root = rank
00302     else
00303        root = D_MASTER
00304     end if
00305 
00306     if (myrank == root) then
00307        ! Master simply copies local data
00308        x(M%lg_fmap(1:size(xl))) = xl
00309 
00310        ! Nothing else to do if I am alone
00311        if (M%nparts <= 1) &
00312             return
00313 
00314 !!!  >>> MASTER COULD DO IT BEFORE HAND EARLIER ...
00315        ! Number of freedoms local to processes
00316 !      allocate(nlfs(M%nparts),f_booked(M%nparts))
00317 !      nlfs = 0
00318 !      f_booked = .false.
00319 !      do gf = 1,M%ngf ! walk through global freedoms
00320 !         h = M%hashlook(int(gf/M%hscale)+1)
00321 !         do while (M%hash(h,1) > 0)
00322 !            if (M%hash(h,1) == gf) then ! found this freedom in hash table
00323 !               p = M%eptnmap(M%hash(h,2))
00324 !               if (.not.f_booked(p)) then
00325 !                  nlfs(p) = nlfs(p) + 1
00326 !                  f_booked(p) = .true.
00327 !               end if
00328 !            end if
00329 !            h = h + 1
00330 !         end do ! do while
00331 !         f_booked = .false.
00332 !      end do
00333 !!! <<< MASTER COULD DO IT BEFORE HAND EARLIER
00334        allocate(nlfs(M%nparts))
00335        do p = 1,M%nparts
00336          if (p-1 == root) then! Don't receive from myself
00337            nlfs(p)=0
00338            cycle
00339          endif
00340          call MPI_RECV(nlfs(p),1,MPI_INTEGER, &
00341               p-1, D_TAG_NLF, MPI_COMM_WORLD, status, ierr)
00342        enddo
00343        allocate(ibuf(maxval(nlfs))); ibuf = 0
00344        allocate(fbuf(maxval(nlfs))); fbuf = 0.0_rk
00345 
00346 ! TEMPORARY "static" solution : rewrite!
00347 ! -- possible workaround :
00348 !      allocate p-1 buffers
00349 !      post p-1 requests
00350 !      go into while cycle with MPI_WAITANY
00351 !         fill x vector with arrived data
00352 !      deallocate incoming buffers
00353        do p = 1,M%nparts
00354           if (p-1 == root) & ! Don't receive from myself
00355                cycle
00356           bs = nlfs(p) ! size of receive buffers
00357           ! Get freedoms' global indexes
00358           call MPI_RECV(ibuf, bs, MPI_INTEGER, &
00359                p-1, D_TAG_FMAP, MPI_COMM_WORLD, status, ierr)
00360           ! Get solution
00361           call MPI_RECV(fbuf, bs, MPI_fkind, &
00362                p-1, D_TAG_FDATA, MPI_COMM_WORLD, status, ierr)
00363           x(ibuf(1:bs)) = fbuf(1:bs)
00364        end do
00365        ! Destroy local dynamic data
00366 !       deallocate( &
00367 !            nlfs,    &
00368 !            f_booked,&
00369 !            ibuf,    &
00370 !            fbuf)
00371        deallocate(nlfs,ibuf,fbuf)
00372     ! Workers
00373     else
00374        ! nlf
00375        !!call MPI_ISEND(M%nlf, 1, MPI_INTEGER, &
00376        call MPI_ISEND(M%ninner, 1, MPI_INTEGER, &
00377             root, D_TAG_NLF, MPI_COMM_WORLD, request, ierr)
00378        ! Freedoms' global indexes
00379        !!call MPI_ISEND(M%lg_fmap, M%nlf, MPI_INTEGER, &
00380        call MPI_ISEND(M%lg_fmap, M%ninner, MPI_INTEGER, &
00381             root, D_TAG_FMAP, MPI_COMM_WORLD, request, ierr)
00382        ! Local vector to send
00383        !!call MPI_ISEND(xl, M%nlf, MPI_fkind, &
00384        call MPI_ISEND(xl, M%ninner, MPI_fkind, &
00385             root, D_TAG_FDATA, MPI_COMM_WORLD, request, ierr)
00386        call MPI_WAIT(request,status,ierr)
00387        !write(stream,*) 'WARNING: must wait for this send request to complete!'
00388     end if
00389   end subroutine Vect_Gather
00390 
00391   !--------------------------------------------
00395   subroutine Integer_Vect_Gather(il, ig, M,owner, rank)
00396     implicit none
00397     integer, dimension(:),     intent(in) :: il ! local vector
00398     integer, dimension(:), intent(in out) :: ig ! global vector
00399     type(Mesh),                intent(in) :: M
00400     integer,dimension(:),pointer,optional :: owner ! freedom owner rank 
00401     ! Rank of a process to gather vector on
00402     integer,                optional, intent(in) :: rank
00403     ! Array of number of local freedoms
00404     integer, dimension(:), pointer :: nlfs
00405     logical, dimension(:), pointer :: f_booked
00406     integer :: i, gf, h, p
00407     ! MPI
00408     ! Aux. receive buffers on master side
00409     integer, dimension(:), pointer :: ibuf
00410     integer, dimension(:), pointer :: iibuf
00411     integer, parameter :: D_TAG_NLF   = 222
00412     integer, parameter :: D_TAG_FMAP  = 333
00413     integer, parameter :: D_TAG_FDATA = 444
00414     integer :: bs, ierr, request, status(MPI_STATUS_SIZE)
00415     integer :: root
00416     if (numprocs==1) then
00417       ig=il
00418       if (present(owner)) then
00419         owner=0
00420       endif
00421       return
00422     endif
00423     if (present(rank)) then
00424        root = rank
00425     else
00426        root = D_MASTER
00427     end if
00428     if (myrank == root) then
00429        ! Master simply copies local data
00430        ig(M%lg_fmap(1:size(il))) = il
00431        if (present(owner)) then
00432          owner(M%lg_fmap(1:size(il)))=myrank+1
00433        endif
00434        ! Nothing else to do if I am alone
00435        if (M%nparts <= 1) &
00436             return
00437        allocate(nlfs(M%nparts))
00438        do p = 1,M%nparts
00439          if (p-1 == root) then! Don't receive from myself
00440            nlfs(p)=0
00441            cycle
00442          endif
00443          call MPI_RECV(nlfs(p),1,MPI_INTEGER, &
00444               p-1, D_TAG_NLF, MPI_COMM_WORLD, status, ierr)
00445        enddo
00446        allocate(ibuf(maxval(nlfs))); ibuf = 0
00447        allocate(iibuf(maxval(nlfs))); iibuf = 0.0_rk
00448        do p = 1,M%nparts
00449           if (p-1 == root) & ! Don't receive from myself
00450                cycle
00451           bs = nlfs(p) ! size of receive buffers
00452           ! Get freedoms' global indexes
00453           call MPI_RECV(ibuf, bs, MPI_INTEGER, &
00454                p-1, D_TAG_FMAP, MPI_COMM_WORLD, status, ierr)
00455           ! Get vector
00456           call MPI_RECV(iibuf, bs, MPI_INTEGER, &
00457                p-1, D_TAG_FDATA, MPI_COMM_WORLD, status, ierr)
00458           ig(ibuf(1:bs)) = iibuf(1:bs)
00459           if (present(owner)) then
00460             owner(ibuf(1:bs))=p
00461           endif
00462        end do
00463        deallocate(nlfs,ibuf,iibuf)
00464     ! Workers
00465     else
00466        ! nlf
00467        call MPI_ISEND(M%ninner, 1, MPI_INTEGER, &
00468             root, D_TAG_NLF, MPI_COMM_WORLD, request, ierr)
00469        ! Freedoms' global indexes
00470        call MPI_ISEND(M%lg_fmap, M%ninner, MPI_INTEGER, &
00471             root, D_TAG_FMAP, MPI_COMM_WORLD, request, ierr)
00472        ! Local vector to send
00473        call MPI_ISEND(il, M%ninner, MPI_INTEGER, &
00474             root, D_TAG_FDATA, MPI_COMM_WORLD, request, ierr)
00475        call MPI_WAIT(request,status,ierr)
00476     end if
00477   end subroutine Integer_Vect_Gather
00478 
00479   subroutine Print_Glob_Vect(x,M,text,rows,chk_endind)
00480     float(kind=rk),dimension(:),intent(in) :: x
00481     type(Mesh),intent(in) :: M ! Mesh
00482     character*(*),intent(in) :: text
00483     logical,intent(in),optional :: rows
00484     integer,intent(in),optional :: chk_endind
00485     logical :: rw
00486     float(kind=rk),dimension(:),pointer :: x_glob
00487     integer :: i,ierr,ei,fmap_size
00488     allocate(x_glob(M%ngf))
00489     call Vect_Gather(x,x_glob,M)
00490     if (present(rows).and.rows) then
00491       rw=.true.
00492     else
00493       rw=.false.
00494     endif
00495     if (present(chk_endind)) then
00496       ei=chk_endind
00497     else
00498       ei=M%ngf
00499     endif
00500     if (ismaster()) then
00501       if (rw) then
00502         write(stream,*) text
00503         do i=1,M%ngf
00504           write(stream,*)i,':',x_glob(i)
00505         enddo
00506       else
00507         write(stream,*) text,x_glob
00508       endif
00509     end if
00510     ! Perform also integrity check on the overlap:
00511     call MPI_BCAST(x_glob,M%ngf,MPI_fkind,D_MASTER,MPI_COMM_WORLD,ierr)
00512     fmap_size=0
00513     if (associated(M%gl_fmap)) fmap_size=size(M%gl_fmap)
00514     do i=1,fmap_size
00515       if (M%gl_fmap(i)/=0.and.M%gl_fmap(i)<=ei) then
00516         if (abs(x(M%gl_fmap(i))-x_glob(i))>1d-14) then
00517           write (*,*)'!!!!!### value mismatch on subd.',myrank,' globind=',i,&
00518           ' loc:',x(M%gl_fmap(i)),' glob:',x_glob(i)
00519           !call DOUG_abort('value mismatch'i,36)
00520         endif
00521       endif
00522     enddo
00523     deallocate(x_glob)
00524   end subroutine Print_Glob_Vect
00525 
00526   !--------------------------------------------------------------------
00527   ! Exchange RHS values
00528   !--------------------------------------------------------------------
00529   subroutine Vect_exchangeIntf(x, M)
00530     implicit none
00531 
00532     float(kind=rk), dimension(:), intent(in out) :: x ! RHS
00533     type(Mesh),                   intent(in)     :: M ! Mesh
00534 
00535     integer :: i, j, le, ge, lf, gf, lfp
00536     integer :: ierr, status(MPI_STATUS_SIZE)
00537     integer :: in_req_count, out_req_count
00538     integer :: n, count, pp
00539     integer, dimension(:), pointer :: in_reqs, out_reqs
00540 
00541     if (size(x) /= M%nlf) &
00542          call DOUG_abort('[Vect_assembleFromElem] : SEVERE : size(x)'//&
00543          ' /= M%nlf',-1)
00544 
00545     ! We've done if we are alone
00546     if (M%nparts <= 1) then
00547        return
00548     end if
00549 
00550     ! Initialize structures needed for communication
00551     call CommStructs_init(M)
00552 
00553     ! First, start non-blocking receives for RHSs from neighbours
00554     ! RHS - initialise receives
00555     allocate(in_reqs(M%nnghbrs), out_reqs(M%nnghbrs))
00556     in_req_count = 0
00557     do pp = 1,M%nparts
00558        if (M%nfreesend_map(pp) /= 0) then
00559           i = pid2indxV(pp)
00560           n = M%nfreesend_map(pp)
00561           in_req_count = in_req_count + 1
00562           call MPI_IRECV(inbufsV(i)%arr, n, MPI_fkind, &
00563                pp-1, D_TAG_FREE_INTERFFREE, MPI_COMM_WORLD, in_reqs(in_req_count), ierr)
00564        end if
00565     end do
00566 
00567     ! RHS - nonblockingly send
00568     out_req_count = 0
00569     do pp = 1,M%nparts
00570        if (M%nfreesend_map(pp) /= 0) then
00571           i = pid2indxV(pp)
00572           n = M%nfreesend_map(pp)
00573           outbufsV(i)%arr = x(fexchindxV(1:n,i))
00574           out_req_count = out_req_count + 1
00575           call MPI_ISEND(outbufsV(i)%arr, n, MPI_fkind, &
00576                pp-1, D_TAG_FREE_INTERFFREE, MPI_COMM_WORLD, out_reqs(out_req_count), ierr)
00577        end if
00578     end do
00579     ! x - wait for neighbours' interface freedoms
00580     do while (.true.)
00581        call MPI_WAITANY(in_req_count, in_reqs, i, status, ierr)
00582        if (i /= MPI_UNDEFINED) then
00583           count = M%nfreesend_map(M%nghbrs(i)+1)
00584           x(fexchindxV(1:count,i)) = &
00585                x(fexchindxV(1:count,i)) + inbufsV(i)%arr
00586        else
00587           exit
00588        end if
00589     end do
00590     call MPI_WAITALL(out_req_count, out_reqs, MPI_STATUSES_IGNORE, ierr)
00591 
00592     deallocate(in_reqs, out_reqs)
00593 
00594     ! Deallocate auxiliary data structures
00595     ! helped to assist with pmvm
00596     call CommStructs_destroy()
00597 
00598   end subroutine Vect_exchangeIntf
00599   
00600   !=======================================================
00603   subroutine Vect_readAndBroadcastRHS(b, Msh) 
00604         use globals
00605         implicit none
00606 
00607         type(Mesh),                intent(in) :: Msh !< Mesh
00608         float(kind=rk), dimension(:), pointer :: b !< local RHS
00609 
00610         float(kind=rk), dimension(:), pointer  :: x 
00611         integer ierr, i
00612 
00613         allocate(x(Msh%ngf))
00614 
00615         ! Read RHS data, if master
00616         if (ismaster()) &
00617           call Vect_ReadFromFile(x, mctls%assembled_rhs_file, mctls%assembled_rhs_format)
00618  
00619         print *, "BROADCAST", size(x)
00620         ! Broadcast the vector from master
00621         call MPI_BCAST(x, size(x), MPI_fkind, D_MASTER, MPI_COMM_WORLD, ierr)
00622 
00623         ! map global RHS to local RHS
00624         allocate( b (size(Msh%lg_fmap)) )
00625         do i = 1,size(Msh%lg_fmap)
00626           b(i) = x(Msh%lg_fmap(i))
00627         end do
00628         print *, "REMAP", size(b)
00629      
00630     deallocate(x)
00631     
00632   end subroutine Vect_readAndBroadcastRHS
00633   
00634   !=======================================================
00635   !
00636   ! Utility subroutines.
00637   !
00638   !-------------------------------------------------------
00639   ! Allocate and initialise data structures used to assist
00640   ! in vector assamblege at communications with neighbours
00641   !-------------------------------------------------------
00642   subroutine CommStructs_init(M)
00643     implicit none
00644 
00645     ! Mesh
00646     type(Mesh),              intent(in) :: M
00647 
00648     integer, dimension(:,:), pointer :: booked
00649     integer,   dimension(:), pointer :: counters
00650     integer :: p, j, h, lf, gf, ge, ptn, indx, n, f
00651 
00652     ! <<<
00653     ! Fill in indexes of freedoms to be
00654     ! exchanged between processors
00655     allocate(fexchindxV(maxval(M%nfreesend_map),M%nnghbrs))
00656     fexchindxV = 0
00657 
00658     ! Map from processes's ids to indeces in 'fexchindxV[:,pid2indxV(:)]'
00659     allocate(pid2indxV(M%nparts))
00660     pid2indxV = 0
00661     do p = 1,M%nparts
00662        do j = 1,M%nnghbrs
00663           if (M%nghbrs(j) == p-1) then
00664              pid2indxV(p) = j
00665              exit
00666           end if
00667        end do
00668     end do
00669 
00670     allocate(booked(M%nlf,M%nnghbrs))
00671     booked = 0
00672     allocate(counters(M%nnghbrs))
00673     counters = 0
00674     do lf = 1,M%nlf
00675        ! interface freedom
00676        if (M%inner_interf_fmask(lf) == D_FREEDOM_INTERF) then
00677           gf = M%lg_fmap(lf)
00678           h = M%hashlook(int(gf/M%hscale)+1)
00679           do while (M%hash(h,1) > 0)
00680              if (M%hash(h,1) == gf) then
00681                 ge = M%hash(h,2)
00682                 ptn = M%eptnmap(ge)
00683                 indx = pid2indxV(ptn)
00684                 if (indx /= 0) then
00685                    if (booked(lf,indx) /= 1) then ! book this freedom
00686                       booked(lf,indx) = 1
00687                       counters(indx) = counters(indx) + 1
00688                       fexchindxV(counters(indx),indx) = lf
00689                    end if
00690                 end if
00691              end if
00692              h = h + 1
00693           end do ! do while
00694        end if
00695     end do
00696 
00697     ! Bufers for incoming and outgoing messages
00698     allocate(inbufsV(M%nnghbrs), outbufsV(M%nnghbrs))
00699     do p = 1,M%nparts
00700        if (M%nfreesend_map(p) /= 0) then
00701           j = pid2indxV(p)
00702           n = M%nfreesend_map(p)
00703           allocate( inbufsV(j)%arr(n))
00704           allocate(outbufsV(j)%arr(n))
00705        end if
00706     end do
00707 
00708     ! Auxiliary arrays has been initialised
00709     D_COMMSTRUCTS_INITEDV = .true.
00710 
00711     deallocate(counters, booked)
00712 
00713   end subroutine CommStructs_init
00714 
00715 
00716   !------------------------------------
00717   ! Deallocate data structures used to
00718   ! assist with vector assemblage
00719   !------------------------------------
00720   subroutine CommStructs_destroy()
00721     implicit none
00722 
00723     integer :: i
00724 
00725     if (associated(fexchindxV)) deallocate(fexchindxV)
00726     if (associated(pid2indxV))  deallocate(pid2indxV)
00727 
00728     ! Destroy incoming buffers
00729     if (associated(inbufsV)) then
00730        do i = 1,size(inbufsV)
00731           if (associated(inbufsV(i)%arr)) deallocate(inbufsV(i)%arr)
00732        end do
00733        deallocate(inbufsV)
00734     end if
00735 
00736     ! Destroy outgoing buffers
00737     if (associated(outbufsV)) then
00738        do i = 1,size(outbufsV)
00739           if (associated(outbufsV(i)%arr)) deallocate(outbufsV(i)%arr)
00740        end do
00741        deallocate(outbufsV)
00742     end if
00743 
00744     D_COMMSTRUCTS_INITEDV = .false.
00745 
00746   end subroutine CommStructs_destroy
00747   !========================
00748   !
00749   ! I/O subroutines
00750   !
00751   !-----------------------------
00752   ! Print out double vector
00753   !-----------------------------
00754   subroutine DVect_Print(x, str)
00755     implicit none
00756     float(kind=rk), dimension(:), intent(in) :: x
00757     character*(*),      optional, intent(in) :: str
00758 
00759     integer :: i, n
00760 
00761     n = size(x)
00762     if (present(str)) then
00763        write(stream,'(/a,i6,a)') str//' :size [',n,']:'
00764     else
00765        write(stream,'(/a,i6,a)') 'vector :size [',n,']:'
00766     end if
00767     do i = 1,n
00768        write(stream, '(a,i6,a,e21.14)') ' [',i,']=',x(i)
00769     end do
00770     call flush(stream)
00771   end subroutine DVect_Print
00772 
00773 
00774   !-----------------------------
00775   ! Print out integer vector
00776   !-----------------------------
00777   subroutine IVect_Print(x, str)
00778     implicit none
00779     integer,   dimension(:), intent(in) :: x
00780     character*(*), optional, intent(in) :: str
00781 
00782     integer :: i, n
00783 
00784     n = size(x)
00785     if (present(str)) then
00786        write(stream,'(/a,i6,a)') str//' :size [',n,']:'
00787     else
00788        write(stream,'(/a,i6,a)') 'vector :size [',n,']:'
00789     end if
00790     do i = 1,n
00791        write(stream, '(a,i6,a,i10)') ' [',i,']=',x(i)
00792     end do
00793     call flush(stream)
00794   end subroutine IVect_Print
00795 
00796 
00797   !-----------------------------
00798   ! Print out integer vector
00799   !-----------------------------
00800   subroutine I1Vect_Print(x, str)
00801     implicit none
00802     integer(kind=1), dimension(:), intent(in) :: x
00803     character*(*),       optional, intent(in) :: str
00804 
00805     integer :: i, n
00806 
00807     n = size(x)
00808     if (present(str)) then
00809        write(stream,'(/a,i6,a)') str//' :size [',n,']:'
00810     else
00811        write(stream,'(/a,i6,a)') 'vector :size [',n,']:'
00812     end if
00813     do i = 1,n
00814        write(stream, '(a,i6,a,i10)') ' [',i,']=',x(i)
00815     end do
00816     call flush(stream)
00817   end subroutine I1Vect_Print
00818 
00819   !! Remap x to y
00820   subroutine Vect_remap(x,y,map,dozero)
00821      use RealKind
00822      
00823      float(kind=rk), intent(in) :: x(:)
00824      float(kind=rk), intent(out) :: y(:)
00825      integer, intent(in) :: map(:)
00826      logical, optional :: dozero
00827      
00828      integer :: i
00829 
00830      if (present(dozero)) y=0.0_rk
00831      
00832      do i=lbound(x,1),ubound(x,1)
00833         if (map(i)/=0) y(map(i))=x(i)
00834      enddo
00835   end subroutine Vect_remap
00836 
00837   !-----------------------------
00840   subroutine Vect_ReadFromFile(x, filename, format)
00841     implicit none
00842     
00843     float(kind=rk), dimension(:), pointer :: x !< the vector; before calling, the vector must be dimensioned with correct size
00844     character*(*), intent(in)             :: filename !< name of the file to read in
00845     integer, intent(in), optional         :: format !< In which format is the input data (default: D_RHS_BINARY)
00846 
00847     integer :: fmt
00848 
00849     ! Binary format is default
00850     if (.not.present(format)) then
00851       fmt = D_RHS_BINARY
00852     else
00853       fmt = format
00854     endif
00855 
00856     if (fmt == D_RHS_TEXT) then
00857       call Vect_ReadFromFile_Text(filename, x)
00858     elseif (fmt == D_RHS_BINARY) then
00859       call Vect_ReadFromFile_Binary(filename, x)
00860 #ifdef HAVE_LIBFXDR
00861     elseif (fmt == D_RHS_XDR) then
00862       call Vect_ReadFromFile_XDR(filename, x)
00863 #endif
00864     else
00865       call DOUG_abort('[Vect_ReadFromFile] Data format not recognized. 0=text, 1=binary, 2=XDR (if compiled in)', -1)
00866     endif
00867 
00868   end subroutine Vect_ReadFromFile
00869   
00870   !-----------------------------
00873   subroutine Vect_ReadFromFile_Text(filename, x)
00874     implicit none
00875     
00876     float(kind=rk), dimension(:), pointer :: x !< the vector; before calling, the vector must be dimensioned with correct size
00877     character*(*), intent(in)             :: filename !< name of the file to read in
00878     logical :: found
00879     integer :: iounit, n, i
00880  
00881     call FindFreeIOUnit(found, iounit)
00882     if (.NOT.found) &
00883       call DOUG_abort('[Vect_ReadFromFile_Text] : No free IO-Unit', -1)
00884  
00885     open(iounit,FILE=trim(filename),STATUS='OLD',FORM='FORMATTED', ERR=444)
00886     read(iounit, '(i9)', END=500) n
00887     if (n /= size(x)) &
00888       call DOUG_abort('[Vect_ReadFromFile_Text] : Number of vector elements in file is not as expected.', -1)
00889     
00890     do i=1,size(x)
00891       read(iounit, '(e23.16)', END=500) x(i)
00892     enddo
00893     return
00894     
00895 444 call DOUG_abort('[Vect_ReadFromFile_Text] : Unable to open vector file: '//trim(filename)//' ', -1)
00896 500 call DOUG_abort('[Vect_ReadFromFile_Text] : End of file reached to early.', -1)
00897   
00898   end subroutine Vect_ReadFromFile_Text
00899 
00900   !-----------------------------
00903   subroutine Vect_ReadFromFile_Binary(filename, x)
00904     implicit none
00905     
00906     float(kind=rk), dimension(:), pointer :: x !< the vector; before calling, the vector must be dimensioned with correct size
00907     character*(*), intent(in)             :: filename !< name of the file to read in
00908     logical :: found
00909     integer :: iounit, n, i
00910  
00911     call FindFreeIOUnit(found, iounit)
00912     if (.NOT.found) &
00913       call DOUG_abort('[Vect_ReadFromFile_Binary] : No free IO-Unit', -1)
00914  
00915     open(iounit, FILE=trim(filename), STATUS='OLD', FORM='UNFORMATTED', ERR=444)
00916     read (iounit, END=500) (x(i), i = 1,size(x))
00917     close(iounit)
00918     return
00919     
00920 444 call DOUG_abort('[Vect_ReadFromFile_Binary] : Unable to open vector file: '//trim(filename)//' ', -1)
00921 500 call DOUG_abort('[Vect_ReadFromFile_Binary] : End of file reached to early.', -1)
00922   
00923   end subroutine Vect_ReadFromFile_Binary
00924 
00925 #ifdef HAVE_LIBFXDR  
00926   !-----------------------------
00929   subroutine Vect_ReadFromFile_XDR(filename, x)
00930     implicit none
00931     include 'fxdr.inc' 
00932     
00933     float(kind=rk), dimension(:), pointer :: x !< the vector; before calling, the vector must be dimensioned with correct size
00934     character*(*), intent(in)             :: filename !< name of the file to read in
00935     integer :: fHandle, ierr, n, i
00936  
00937     fHandle = initxdr( trim(filename), 'r', .FALSE. )
00938     ierr  = ixdrint( fHandle, n )
00939 
00940     if (n /= size(x)) &
00941       call DOUG_abort('[Vect_ReadFromFile_XDR] : Number of vector elements in file is not as expected.', -1)
00942  
00943     do i=1,size(x)
00944       ierr = ixdrdouble( fHandle, x(i) )
00945     enddo
00946   
00947   end subroutine Vect_ReadFromFile_XDR
00948 
00949 #endif
00950 
00951 end module Vect_mod