DOUG 0.2

SpMtx_operation.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 ! Include some useful operations:
00024 !     Dense to Sparse
00025 !     Sparse to Dense
00026 !     Full AB operation
00027 !----------------------------------------------------------
00028 Module SpMtx_operation
00029 
00030   use SpMtx_class
00031   use SpMtx_permutation
00032   use SpMtx_util
00033   use SpMtx_op_Ax
00034   use Vect_mod
00035   use RealKind
00036   use globals
00037 
00038   Implicit None
00039 
00040 #include<doug_config.h>
00041 
00042 #ifdef D_COMPLEX
00043 #define float complex
00044 #else
00045 #define float real
00046 #endif
00047 
00048 !!$  type iarray
00049 !!$     integer, dimension(:), pointer :: arr
00050 !!$  end type iarray
00051 
00052   type farray
00053      float(kind=rk), dimension(:), pointer :: arr
00054   end type farray
00055 
00056   ! Important to know for comm. after preconditioning:
00057   logical :: D_PMVM_USING_AGHOST = .false.
00058 
00059   type OperationCache
00061     logical :: D_PMVM_AUXARRS_INITED
00062 
00063     type(farray), dimension(:), pointer :: inbufs
00064 
00065     type(farray), dimension(:), pointer :: outbufs
00066 
00068     integer, dimension(:,:), pointer :: fexchindx
00069 
00070     integer,   dimension(:), pointer :: pid2indx
00071   end type OperationCache
00072 
00073 CONTAINS
00074 !---------------------------------------------------
00075 ! Convert Dense matrix to Sparse Matrix type
00076 !         M - Regular Matrix
00077 !         A - Sparse Matrix Type
00078 !---------------------------------------------------
00079   Function SpMtx_DenseToSparse(M) result(A)
00080     Implicit None
00081     Type(SpMtx)                          :: A   !Sparse matrix
00082     Real(kind=rk), intent(in), dimension(:,:):: M   !dense matrix
00083     Integer                                  :: i, j, nnz, n
00084     !- - - - - - - - - - - - - - - - - - - - - - - - -
00085     nnz=0                            !count number of non-zero elements
00086     do i=1,size(M,dim=1)
00087       do j=1,size(M, dim=2)
00088         if (M(i,j) /= 0.) nnz=nnz+1
00089       end do
00090     end do
00091     A = SpMtx_newInit(nnz)
00092     A%nrows=size(M,dim=1)
00093     A%ncols=size(M,dim=2)
00094     n=0
00095     do i=1,A%nrows                     !Adding values to sparse matrix
00096       do j=1,A%ncols
00097         if (M(i,j) /= 0.) then
00098           n=n+1
00099           A%indi(n)=i
00100           A%indj(n)=j
00101           A%val(n)=M(i,j)
00102         end if
00103       end do
00104     end do
00105     if (n /= nnz) &
00106          call DOUG_abort('[SpMtx_DenseToSparse] : SEVERE : Dimensions'//&
00107          ' conflict.', -1)
00108   End Function SpMtx_DenseToSparse
00109 !---------------------------------------------------
00110 ! Convert Sparse Matrix to Dense one
00111 !         M - Regular Matrix
00112 !         A - Sparse Matrix Type
00113 !   M must be allocated before.
00114 !---------------------------------------------------
00115   Function SpMtx_SparseToDense(A) result(M)
00116     Implicit None
00117     Type(SpMtx),                   intent(in) :: A !Sparse matrix
00118     Real(kind=rk), dimension(A%nrows,A%ncols) :: M !dense matrix
00119     Integer                                   :: i !counter
00120     !- - - - - - - - - - - - - - - - - - - - -
00121     M=0.0_rk
00122     do i=1, A%nnz
00123       M(A%indi(i), A%indj(i))=A%val(i)
00124     end do
00125   End Function SpMtx_SparseToDense
00126 !------------------------------------------------
00127 ! Full operation AB
00128 !  A, B, AB : sparse Matrix
00129 !------------------------------------------------
00130   Function SpMtx_Full_AB(A,B) result(AB)
00131     Implicit None
00132     Type(SpMtx), intent(in):: A, B !matrix (in)
00133     Type(SpMtx)            :: AB   !result matrix (out)
00134     Integer                    :: counter, i, j
00135     !- - - - - - - - - - - - - - - - - -
00136     counter = 0
00137     do i = 1, A%nrows
00138       do j = 1, B%ncols
00139         counter=counter+1
00140         AB%indi(counter) = i
00141         AB%indj(counter) = j
00142         AB%val(counter) = get_product(A,B,i,j)
00143       enddo
00144     enddo
00145   End Function SpMtx_Full_AB
00146 !----------------------------------------------------------
00147 ! Useful function for Full_AB
00148 !----------------------------------------------------------
00149   function get_product(A,B,ci,cj) result(res)
00150     implicit none
00151     type(SpMtx), intent(in) :: A, B   !sparse matrix
00152     integer, intent(in)        :: ci, cj !matrix element indexes (i,j)
00153     real(kind=rk)              :: res    !result
00154     integer                    :: i, j   !counters
00155     !- - - - - - - - - - - - - - - - - - -
00156     res = 0
00157     do i = 1, A%nnz
00158       if(A%indi(i) == ci) then
00159         do j = 1, B%nnz
00160           if(B%indj(j) == cj) then
00161             if(B%indi(j) == A%indj(i)) then
00162               res = res + B%val(j) * A%val(i)
00163             end if
00164           end if
00165         end do
00166       end if
00167     end do
00168   end function get_product
00169 
00170   subroutine SpMtx_pmvm_elemental(y,A,x,M,C)
00171     implicit none
00172     type(SpMtx),                      intent(in) :: A ! System matrix
00173     float(kind=rk),dimension(:),pointer          :: x ! Vector
00174     type(Mesh),                       intent(in) :: M ! Mesh
00175     float(kind=rk),dimension(:),pointer          :: y ! Result
00176     type(OperationCache) :: C
00177     integer :: i, n, p
00178     ! MPI
00179     integer,dimension(:),pointer :: in_reqs
00180     integer                      :: ierr,out_req,status(MPI_STATUS_SIZE)
00181     integer,parameter            :: D_TAG_FREE_INTERFFREE=676
00182     float(kind=rk),dimension(:),pointer :: yp,xp
00183 
00184     y=0.0_rk
00185 
00186     ! Innerface freedoms matrix-vector multiplication
00187     !  ---------- -----------
00188     ! |^ interf. |  interf./ |
00189     ! |         v|  inner    |
00190     !  ----------+-----------
00191     ! | inner/   |           |
00192     ! | interf.  |   inner   |
00193     ! |          |           |
00194     !  ---------- ------------
00195     call SpMtx_mvm(A,x,y, &
00196          A%mtx_bbs(1,1),A%mtx_bbe(1,1))
00197     !call Vect_Print(y, 'after 1:1 mult:')
00198 
00199     ! look -- flipping in file datatypes/ElemMtxs_assemble.f90
00200     ! todo: is the flipping still needed, or correct at all?
00201     ! Inner/iterface freedoms matrix-vector multiplication
00202     !  --------- ----------
00203     ! | interf. | interf./ |
00204     ! |         | inner    |
00205     !  ---------+----------
00206     ! |^ inner/ |          |
00207     ! |  interf.| inner    |
00208     ! |        v|          |
00209     !  ---------- ----------
00210     call SpMtx_mvm(A,x,y, &
00211          A%mtx_bbs(A%nblocks+1,1),A%mtx_bbe(A%nblocks+1,1))
00212     ! y - initialise receives
00213     allocate(in_reqs(M%nnghbrs))
00214     do p=1,M%nparts
00215       if (M%nfreesend_map(p)/=0) then
00216         i=C%pid2indx(p)
00217         n=M%nfreesend_map(p)
00218         call MPI_IRECV(C%inbufs(i)%arr,n,MPI_fkind, &
00219              p-1,D_TAG_FREE_INTERFFREE,MPI_COMM_WORLD,in_reqs(i),ierr)
00220       endif
00221     enddo
00222     ! y - nonblockingly send
00223     do p=1,M%nparts
00224       if (M%nfreesend_map(p)/=0) then
00225         i=C%pid2indx(p)
00226         n=M%nfreesend_map(p)
00227         C%outbufs(i)%arr(1:n)=y(C%fexchindx(1:n,i))
00228         call MPI_ISEND(C%outbufs(i)%arr,n,MPI_fkind, &
00229              p-1,D_TAG_FREE_INTERFFREE,MPI_COMM_WORLD,out_req,ierr)
00230       endif
00231     enddo
00232     ! Innerface/inner freedoms matrix-vector multiplication
00233     !  --------- -----------
00234     ! | interf. |^ interf./ |
00235     ! |         |  inner   v|
00236     !  ---------+-----------
00237     ! | inner/  |           |
00238     ! | interf. |   inner   |
00239     ! |         |           |
00240     !  --------- ------------
00241     call SpMtx_mvm(A,x,y, &
00242          A%mtx_bbs(1,A%nblocks+1),A%mtx_bbe(1,A%nblocks+1))
00243 
00244     ! Inner freedoms matrix-vector multiplication
00245     !  ---------- ----------
00246     ! |  interf. | interf./ |
00247     ! |          | inner    |
00248     !  ----------+----------
00249     ! |  inner/  |^         |
00250     ! |  interf. | inner    |
00251     ! |          |         v|
00252     !  ---------- ----------
00253     call SpMtx_mvm(A,x,y, &
00254          A%mtx_bbs(A%nblocks+1,A%nblocks+1),A%mtx_bbe(A%nblocks+1,A%nblocks+1))
00255     !call Vect_Print(y, 'after 4:4 mult:')
00256 
00257     ! y - wait for neighbours' interface freedoms
00258     do while (.true.)
00259       call MPI_WAITANY(M%nnghbrs,in_reqs,i,status,ierr)
00260       if (i/=MPI_UNDEFINED) then
00261         n=M%nfreesend_map(M%nghbrs(i)+1)
00262         y(C%fexchindx(1:n,i))=&
00263              y(C%fexchindx(1:n,i))+C%inbufs(i)%arr(1:n)
00264         !write(stream,*) 'got from:', M%nghbrs(i),'buf:', inbufs(i)%arr
00265         !write(stream,*) 'y=',y
00266       else
00267         exit
00268       endif
00269     enddo
00270     deallocate(in_reqs)
00271   end subroutine SpMtx_pmvm_elemental
00272 
00273   subroutine SpMtx_pmvm_assembled(y,A,x,M,C)
00274     implicit none
00275 
00276     type(SpMtx),                      intent(in) :: A ! System matrix
00277     float(kind=rk),dimension(:),pointer          :: x ! Vector
00278     type(Mesh),                       intent(in) :: M ! Mesh
00279     float(kind=rk),dimension(:),pointer          :: y ! Result
00280     type(OperationCache) :: C
00281     integer :: i, n, p
00282     ! MPI
00283     integer,dimension(:),pointer :: in_reqs
00284     integer                      :: ierr,out_req,status(MPI_STATUS_SIZE)
00285     integer,parameter            :: D_TAG_FREE_INTERFFREE=676
00286     float(kind=rk),dimension(:),pointer :: yp,xp
00287 
00288     y=0.0_rk
00289 
00290     ! Innerface freedoms matrix-vector multiplication
00291     !  ---------- -----------
00292     ! |^ interf. |  interf./ |
00293     ! |         v|  inner    |
00294     !  ----------+-----------
00295     ! | inner/   |           |
00296     ! | interf.  |   inner   |
00297     ! |          |           |
00298     !  ---------- ------------
00299     call SpMtx_mvm(A,x,y, &
00300          A%mtx_bbs(1,1),A%mtx_bbe(1,1))
00301     !call vect_print(y, 'after 1:1 mult:')
00302     call SpMtx_mvm(A,x,y, &
00303        A%mtx_bbs(1,2),A%mtx_bbe(1,2))
00304 
00305     ! y - initialise receives
00306     allocate(in_reqs(M%nnghbrs))
00307     do i=1,M%nnghbrs
00308       n=M%ax_recvidx(i)%ninds
00309       p=M%nghbrs(i)
00310 !write(stream,*) '**** starting non-blocking recv from ',p
00311       call MPI_IRECV(C%inbufs(i)%arr,n,MPI_fkind, &
00312                p,D_TAG_FREE_INTERFFREE,MPI_COMM_WORLD,in_reqs(i),ierr)
00313     enddo
00314     ! y - nonblockingly send
00315     do i=1,M%nnghbrs
00316       n=M%ax_sendidx(i)%ninds
00317       p=M%nghbrs(i)
00318       C%outbufs(i)%arr(1:M%ax_sendidx(i)%ninds)=y(M%ax_sendidx(i)%inds)
00319       call MPI_ISEND(C%outbufs(i)%arr,n,MPI_fkind, &
00320                p,D_TAG_FREE_INTERFFREE,MPI_COMM_WORLD,out_req,ierr)
00321 !write(stream,*) '**** sending to ',p,outbufs(i)%arr(1:M%ax_sendidx(i)%ninds)
00322     enddo
00323     call SpMtx_mvm(A,x,y, &
00324          A%mtx_bbs(2,1),A%mtx_bbe(2,1))
00325 
00326     ! Inner freedoms matrix-vector multiplication
00327     !  ---------- ----------
00328     ! |  interf. | interf./ |
00329     ! |          | inner    |
00330     !  ----------+----------
00331     ! |  inner/  |^         |
00332     ! |  interf. | inner    |
00333     ! |          |         v|
00334     !  ---------- ----------
00335     call SpMtx_mvm(A,x,y,A%mtx_bbs(2,2),A%mtx_bbe(2,2))
00336     !call Vect_Print(y, 'after 4:4 mult:')
00337 
00338 !!!write(stream,*)'----- q before comm:------ ',y
00339     ! y - wait for neighbours' interface freedoms
00340     do while (.true.)
00341       call MPI_WAITANY(M%nnghbrs,in_reqs,i,status,ierr)
00342       if (i/=MPI_UNDEFINED) then
00343 !write(stream,*)'**** received from ',M%nghbrs(i),inbufs(i)%arr(1:M%ax_recvidx(i)%ninds)
00344         y(M%ax_recvidx(i)%inds)=C%inbufs(i)%arr(1:M%ax_recvidx(i)%ninds)
00345       else
00346         exit
00347       endif
00348     enddo
00349     deallocate(in_reqs)
00350   end subroutine SpMtx_pmvm_assembled
00351 
00352   subroutine SpMtx_pmvm_assembled_ol0(y,A,x,M,C)
00353     implicit none
00354 
00355     type(SpMtx),                      intent(in) :: A ! System matrix
00356     float(kind=rk),dimension(:),pointer          :: x ! Vector
00357     type(Mesh),                       intent(in) :: M ! Mesh
00358     float(kind=rk),dimension(:),pointer          :: y ! Result
00359     type(OperationCache) :: C
00360     integer :: i, n, p
00361     ! MPI
00362     integer,dimension(:),pointer :: in_reqs
00363     integer                      :: ierr,out_req,status(MPI_STATUS_SIZE)
00364     integer,parameter            :: D_TAG_FREE_INTERFFREE=676
00365     float(kind=rk),dimension(:),pointer :: yp,xp
00366 
00367     y=0.0_rk
00368     ! y - initialise receives
00369     allocate(in_reqs(M%nnghbrs))
00370     do i=1,M%nnghbrs
00371       n=M%ax_recvidx(i)%ninds
00372       p=M%nghbrs(i)
00373 !write(stream,*) '**** starting non-blocking recv from ',p
00374       call MPI_IRECV(C%inbufs(i)%arr,n,MPI_fkind, &
00375                p,D_TAG_FREE_INTERFFREE,MPI_COMM_WORLD,in_reqs(i),ierr)
00376     enddo
00377     ! x - nonblockingly send
00378     do i=1,M%nnghbrs
00379       n=M%ax_sendidx(i)%ninds
00380       p=M%nghbrs(i)
00381       C%outbufs(i)%arr(1:M%ax_sendidx(i)%ninds)=x(M%ax_sendidx(i)%inds)
00382       call MPI_ISEND(C%outbufs(i)%arr,n,MPI_fkind, &
00383                p,D_TAG_FREE_INTERFFREE,MPI_COMM_WORLD,out_req,ierr)
00384 !write(stream,*) '**** sending to ',p,outbufs(i)%arr
00385     enddo
00386     ! perform all the calculations on the subdomain
00387     call SpMtx_mvm(A,x,y, &
00388          A%mtx_bbs(1,1),A%mtx_bbe(1,1))
00389     call SpMtx_mvm(A,x,y, &
00390          A%mtx_bbs(2,1),A%mtx_bbe(2,1))
00391     call SpMtx_mvm(A,x,y, &
00392          A%mtx_bbs(1,2),A%mtx_bbe(1,2))
00393     call SpMtx_mvm(A,x,y, &
00394          A%mtx_bbs(2,2),A%mtx_bbe(2,2))
00395 !!!write(stream,*)'----- q before comm:------ ',y
00396     ! x - wait for neighbours' interface freedoms
00397     do while (.true.)
00398       call MPI_WAITANY(M%nnghbrs,in_reqs,i,status,ierr)
00399       if (i/=MPI_UNDEFINED) then
00400 !write(stream,*)'**** received from ',M%nghbrs(i),inbufs(i)%arr
00401         x(M%ax_recvidx(i)%inds)=C%inbufs(i)%arr(1:M%ax_recvidx(i)%ninds)
00402       else
00403         exit
00404       endif
00405     enddo
00406     deallocate(in_reqs)
00407     ! In the case of OL=0 the part after comm. needs to be done:
00408 !write(stream,*)'y before A(end)x:',y
00409 !write(stream,*)'A end:'
00410 !call SpMtx_printRaw(A=A,startnz=A%mtx_bbe(2,2)+1,endnz=A%nnz)
00411     call SpMtx_mvm(A,x,y,A%mtx_bbe(2,2)+1,A%nnz)
00412 !write(stream,*)'y after A(end)x:',y
00413   end subroutine SpMtx_pmvm_assembled_ol0
00414 
00415   subroutine exch_aggr_nums(aggr,M)
00416     implicit none
00417 
00418     integer,dimension(:),pointer          :: aggr ! Vector
00419     type(Mesh),                       intent(in) :: M ! Mesh
00420     integer :: i, n, p
00421     ! MPI
00422     integer,dimension(:),pointer :: in_reqs,out_reqs
00423     integer                      :: ierr,status(MPI_STATUS_SIZE)
00424     integer,parameter            :: D_TAG_FREE_INTERFFREE=777
00425     ! Input/output bufers for send/receiv freedoms
00426     type(indlist),dimension(:),pointer :: iinbufs
00427     type(indlist),dimension(:),pointer :: ioutbufs
00428 
00429     ! initialise receives
00430     allocate(out_reqs(M%nnghbrs))
00431     allocate(in_reqs(M%nnghbrs))
00432     allocate(iinbufs(M%nnghbrs))
00433     allocate(ioutbufs(M%nnghbrs))
00434     do i=1,M%nnghbrs
00435       n=M%ol_outer(i)%ninds
00436       allocate(iinbufs(i)%inds(n))
00437       p=M%nghbrs(i)
00438 !write(stream,*) '**** starting non-blocking recv from ',p
00439       call MPI_IRECV(iinbufs(i)%inds,n,MPI_INTEGER, &
00440                p,D_TAG_FREE_INTERFFREE,MPI_COMM_WORLD,in_reqs(i),ierr)
00441     enddo
00442     ! x - nonblockingly send
00443     do i=1,M%nnghbrs
00444       n=M%ol_inner(i)%ninds
00445       allocate(ioutbufs(i)%inds(n))
00446       p=M%nghbrs(i)
00447       ioutbufs(i)%inds(1:M%ol_inner(i)%ninds)=aggr(M%ol_inner(i)%inds)
00448       call MPI_ISEND(ioutbufs(i)%inds,n,MPI_INTEGER, &
00449                p,D_TAG_FREE_INTERFFREE,MPI_COMM_WORLD,out_reqs(i),ierr)
00450 !write(stream,*) '**** sending to ',p,outbufs(i)%arr
00451     enddo
00452     ! could actually perform some calculations here... (TODO)
00453     ! wait for neighbours' interface freedoms
00454     do while (.true.)
00455       call MPI_WAITANY(M%nnghbrs,in_reqs,i,status,ierr)
00456       if (i/=MPI_UNDEFINED) then
00457 !write(stream,*)'**** received from ',M%nghbrs(i),inbufs(i)%arr
00458         aggr(M%ol_outer(i)%inds)=iinbufs(i)%inds(1:M%ol_outer(i)%ninds)
00459         deallocate(iinbufs(i)%inds)
00460       else
00461         exit
00462       endif
00463     enddo
00464     deallocate(in_reqs)
00465     do while (.true.)
00466       call MPI_WAITANY(M%nnghbrs,out_reqs,i,status,ierr)
00467       if (i/=MPI_UNDEFINED) then
00468         deallocate(ioutbufs(i)%inds)
00469       else
00470         exit
00471       endif
00472     enddo
00473     deallocate(out_reqs)
00474     deallocate(ioutbufs)
00475     deallocate(iinbufs)
00476     ! do some post-calculations here if needed... todo
00477   end subroutine exch_aggr_nums
00478 
00479   subroutine exch_aggr_nums_ol0(aggr,M)
00480     implicit none
00481 
00482     integer,dimension(:),pointer          :: aggr ! Vector
00483     type(Mesh),                       intent(in) :: M ! Mesh
00484     integer :: i, n, p
00485     ! MPI
00486     integer,dimension(:),pointer :: in_reqs,out_reqs
00487     integer                      :: ierr,status(MPI_STATUS_SIZE)
00488     integer,parameter            :: D_TAG_FREE_INTERFFREE=777
00489     ! Input/output bufers for send/receiv freedoms
00490     type(indlist),dimension(:),pointer :: iinbufs
00491     type(indlist),dimension(:),pointer :: ioutbufs
00492 
00493     ! initialise receives
00494     allocate(out_reqs(M%nnghbrs))
00495     allocate(in_reqs(M%nnghbrs))
00496     allocate(iinbufs(M%nnghbrs))
00497     allocate(ioutbufs(M%nnghbrs))
00498     do i=1,M%nnghbrs
00499       n=M%ax_recvidx(i)%ninds
00500       allocate(iinbufs(i)%inds(n))
00501       p=M%nghbrs(i)
00502 !write(stream,*) '**** starting non-blocking recv from ',p
00503       call MPI_IRECV(iinbufs(i)%inds,n,MPI_INTEGER, &
00504                p,D_TAG_FREE_INTERFFREE,MPI_COMM_WORLD,in_reqs(i),ierr)
00505     enddo
00506     ! x - nonblockingly send
00507     do i=1,M%nnghbrs
00508       n=M%ax_sendidx(i)%ninds
00509       allocate(ioutbufs(i)%inds(n))
00510       p=M%nghbrs(i)
00511       ioutbufs(i)%inds(1:M%ax_sendidx(i)%ninds)=aggr(M%ax_sendidx(i)%inds)
00512       call MPI_ISEND(ioutbufs(i)%inds,n,MPI_INTEGER, &
00513                p,D_TAG_FREE_INTERFFREE,MPI_COMM_WORLD,out_reqs(i),ierr)
00514 !write(stream,*) '**** sending to ',p,outbufs(i)%arr
00515     enddo
00516     ! could actually perform some calculations here... (TODO)
00517     ! wait for neighbours' interface freedoms
00518     do while (.true.)
00519       call MPI_WAITANY(M%nnghbrs,in_reqs,i,status,ierr)
00520       if (i/=MPI_UNDEFINED) then
00521 !write(stream,*)'**** received from ',M%nghbrs(i),inbufs(i)%arr
00522         aggr(M%ax_recvidx(i)%inds)=iinbufs(i)%inds(1:M%ax_recvidx(i)%ninds)
00523         deallocate(iinbufs(i)%inds)
00524       else
00525         exit
00526       endif
00527     enddo
00528     deallocate(in_reqs)
00529     do while (.true.)
00530       call MPI_WAITANY(M%nnghbrs,out_reqs,i,status,ierr)
00531       if (i/=MPI_UNDEFINED) then
00532         deallocate(ioutbufs(i)%inds)
00533       else
00534         exit
00535       endif
00536     enddo
00537     deallocate(out_reqs)
00538     deallocate(ioutbufs)
00539     deallocate(iinbufs)
00540     ! do some post-calculations here if needed... todo
00541   end subroutine exch_aggr_nums_ol0
00542 
00543 !!$  !--------------------------------------
00544 !!$  ! Parallel matrix-vector multiplication
00545 !!$  !--------------------------------------
00546 !!$  subroutine SpMtx_pmvm(A,x,y, M)
00547 !!$    implicit none
00548 !!$
00549 !!$    type(SpMtx),                      intent(in) :: A ! System matrix
00550 !!$    float(kind=rk), dimension(:), intent(in out) :: x ! Vector
00551 !!$    float(kind=rk), dimension(:), intent(in out) :: y ! Result
00552 !!$
00553 !  ==> NEED TO GET RID OF IT HERE
00554 !!$    type(Mesh),     intent(in) :: M ! Mesh
00555 !!$
00556 !!$    integer :: i, n, p, count
00557 !!$
00558 !!$    ! MPI
00559 !!$    integer, dimension(:), pointer :: in_reqs
00560 !!$    integer                        :: ierr, out_req, status(MPI_STATUS_SIZE)
00561 !!$    integer, parameter             :: D_TAG_FREE_INTERFFREE = 666
00562 !!$
00563 !!$
00564 !!$    if (size(x) /= size(y)) &
00565 !!$         call DOUG_abort('[SpMtx_mvm] : size(x) /= size(y)',-1)
00566 !!$
00567 !!$    ! Initialise auxiliary data structures
00568 !!$    ! to assist with pmvm
00569 !!$    if (.not.D_PMVM_AUXARRS_INITED) &
00570 !!$         call pmvmCommStructs_init(A, M)
00571 !!$
00572 !!$
00573 !!$    ! Innerface freedoms matrix-vector multiplication
00574 !!$    !  ---------- -----------
00575 !!$    ! |^ interf. |  interf./ |
00576 !!$    ! |         v|  inner    |
00577 !!$    !  ----------+-----------
00578 !!$    ! | inner/   |           |
00579 !!$    ! | interf.  |   inner   |
00580 !!$    ! |          |           |
00581 !!$    !  ---------- ------------
00582 !!$    call SpMtx_mvm(A, x, y, &
00583 !!$         A%mtx_bbs(1,1), A%mtx_bbe(A%nblocks,A%nblocks))
00584 !!$
00585 !    write(stream,*) 'after 1:1 mult:'
00586 !    do i = 1,M%nlf
00587 !       write(stream,'(a,i2,a,f9.5)') 'y(',i,') =', y(i)
00588 !    end do
00589 !!$
00590 !!$    ! Innerface/inner freedoms matrix-vector multiplication
00591 !!$    !  --------- -----------
00592 !!$    ! | interf. |^ interf./ |
00593 !!$    ! |         |  inner   v|
00594 !!$    !  ---------+-----------
00595 !!$    ! | inner/  |           |
00596 !!$    ! | interf. |   inner   |
00597 !!$    ! |         |           |
00598 !!$    !  --------- ------------
00599 !!$    call SpMtx_mvm(A, x, y, &
00600 !!$         A%mtx_bbs(1,A%nblocks+1), A%mtx_bbe(A%nblocks,2*A%nblocks))
00601 !!$
00602 !    write(stream,*) 'after 3:3 mult:'
00603 !    do i = 1,M%nlf
00604 !       write(stream,'(a,i2,a,f9.5)') 'y(',i,') =', y(i)
00605 !    end do
00606 !!$
00607 !!$
00608 !!$    ! y - initialise receives
00609 !!$    allocate(in_reqs(M%nnghbrs))
00610 !!$    do p = 1,M%nparts
00611 !!$       if (M%nfreesend_map(p) /= 0) then
00612 !!$          i = pid2indx(p)
00613 !!$          n = M%nfreesend_map(p)
00614 !!$          call MPI_IRECV(inbufs(i)%arr, n, MPI_fkind, &
00615 !!$               p-1, D_TAG_FREE_INTERFFREE, MPI_COMM_WORLD, in_reqs(i), ierr)
00616 !!$       end if
00617 !!$    end do
00618 !!$
00619 !!$    ! y - nonblockingly send
00620 !!$    do p = 1,M%nparts
00621 !!$       if (M%nfreesend_map(p) /= 0) then
00622 !!$          i = pid2indx(p)
00623 !!$          n = M%nfreesend_map(p)
00624 !!$          outbufs(i)%arr = y(fexchindx(1:n,i))
00625 !!$          call MPI_ISEND(outbufs(i)%arr, n, MPI_fkind, &
00626 !!$               p-1, D_TAG_FREE_INTERFFREE, MPI_COMM_WORLD, out_req, ierr)
00627 !!$       end if
00628 !!$    end do
00629 !!$
00630 !!$
00631 !!$    ! Inner/iterface freedoms matrix-vector multiplication
00632 !!$    !  --------- ----------
00633 !!$    ! | interf. | interf./ |
00634 !!$    ! |         | inner    |
00635 !!$    !  ---------+----------
00636 !!$    ! |^ inner/ |          |
00637 !!$    ! |  interf.| inner    |
00638 !!$    ! |        v|          |
00639 !!$    !  ---------- ----------
00640 !!$    call SpMtx_mvm(A, x, y, &
00641 !!$         A%mtx_bbs(A%nblocks+1,1), A%mtx_bbe(2*A%nblocks,A%nblocks))
00642 !!$
00643 !    write(stream,*) 'after 2:2 mult:'
00644 !    do i = 1,M%nlf
00645 !       write(stream,'(a,i2,a,f9.5)') 'y(',i,') =', y(i)
00646 !    end do
00647 !!$
00648 !!$    ! Inner freedoms matrix-vector multiplication
00649 !!$    !  ---------- ----------
00650 !!$    ! |  interf. | interf./ |
00651 !!$    ! |          | inner    |
00652 !!$    !  ----------+----------
00653 !!$    ! |  inner/  |^         |
00654 !!$    ! |  interf. | inner    |
00655 !!$    ! |          |         v|
00656 !!$    !  ---------- ----------
00657 !!$    call SpMtx_mvm(A, x, y, &
00658 !!$         A%mtx_bbs(A%nblocks+1,A%nblocks+1),A%mtx_bbe(2*A%nblocks,2*A%nblocks))
00659 !!$
00660 !    write(stream,*) 'after 4:4 mult:'
00661 !    do i = 1,M%nlf
00662 !       write(stream,'(a,i2,a,f9.5)') 'y(',i,') =', y(i)
00663 !    end do
00664 !!$
00665 !!$
00666 !!$    ! y - wait for neighbours' interface freedoms
00667 !!$    do while (1)
00668 !!$       call MPI_WAITANY(M%nnghbrs, in_reqs, i, status, ierr)
00669 !!$       if (i /= MPI_UNDEFINED) then
00670 !!$          count = M%nfreesend_map(M%nghbrs(i)+1)
00671 !!$          y(fexchindx(1:count,i)) = &
00672 !!$               y(fexchindx(1:count,i)) + inbufs(i)%arr
00673 !!$
00674 !          write(stream,*) 'got from:', M%nghbrs(i),'buf:', inbufs(i)%arr
00675 !          write(stream,*) 'y=',y
00676 !!$
00677 !!$       else
00678 !!$          exit
00679 !!$       end if
00680 !!$    end do
00681 !!$
00682 !!$    deallocate(in_reqs)
00683 !!$  end subroutine SpMtx_pmvm
00684 
00685 
00686   !----------------------------------------------
00687   ! General parallel matrix-vector multiplication
00688   !----------------------------------------------
00689 !!$  subroutine SpMtx_gmvm(A,alpha,x,beta,y) ! y = beta*y + alpha*A*x
00690 !!$
00691 !!$  end subroutine SpMtx_gmvm
00692 
00693 
00694   !--------------------------------------
00695   ! Matrix-vector multiplication
00696   !--------------------------------------
00697   subroutine SpMtx_mvm(A, x, y, bbs, bbe)
00698     implicit none
00699 
00700     type(SpMtx),intent(in)              :: A ! System matrix
00701     float(kind=rk),dimension(:),pointer :: x
00702     float(kind=rk),dimension(:),pointer :: y ! resulting vector
00703     integer,intent(in),optional         :: bbs, bbe
00704 
00705     integer :: i, row, col
00706     integer :: s, e
00707 
00708     s = 1
00709     e = A%nnz
00710 
00711     if (present(bbs).and.(.not.present(bbe))) &
00712          call DOUG_abort('[SpMtx_mvm] : SEVERE : "bbe" must be present'//&
00713          ' along with "bbs".',-1)
00714     if (present(bbs)) then
00715        s = bbs
00716        e = bbe
00717     end if
00718 
00719     do i = s,e
00720        row = A%indi(i)
00721        col = A%indj(i)
00722        y(row) = y(row) + A%val(i)*x(col)
00723     end do
00724 
00725   end subroutine SpMtx_mvm
00726 
00728   function SpMtx_addAll(As, koefs) result(C)
00729     type(SpMtx), intent(in) :: As(:)
00730     real(kind=rk), intent(in) :: koefs(:)
00731     type(SpMtx) :: C
00732 
00733     integer :: i, nnz, from, to
00734     integer, dimension(:), pointer :: indi, indj, val
00735 
00736     nnz = sum(As%nnz)
00737 
00738     C = SpMtx_newInit(nnz)
00739     !write(stream,*) nnz, As%nnz, size(As(1)%indi), C%nnz, size(C%indi), size(C%indj)
00740     from = 1
00741     do i=1,size(As)
00742       if (As(i)%nnz==0) cycle
00743        to = from + As(i)%nnz - 1
00744        C%indi(from:to) = As(i)%indi(1:As(i)%nnz)
00745        C%indj(from:to) = As(i)%indj(1:As(i)%nnz)
00746        C%val(from:to) = koefs(i)*As(i)%val(1:As(i)%nnz)
00747        from = to+1
00748     end do
00749 
00750     C%nrows = maxval(As%nrows)
00751     C%ncols = maxval(As%ncols)
00752     
00753     ! add up duplicates
00754     call SpMtx_consolidate(C, .TRUE.)
00755 
00756     ! debug
00757     !call SpMtx_printRaw(C)
00758 
00759   end function SpMtx_addAll
00760 
00762   function SpMtx_add(A,B,ka,kb) result(C)
00763     type(SpMtx), intent(in) :: A,B
00764     real(kind=rk), intent(in) :: ka, kb
00765     type(SpMtx) :: C
00766     
00767     C = SpMtx_addAll((/A,B/),(/ka,kb/))
00768   end function SpMtx_add
00769 
00770 End Module SpMtx_operation
00771 !----------------------------------------------------------
00772 !$Log: SpMtx_operation.f90,v $
00773 !Revision 1.2  2004/04/30 08:47:17  elmo
00774 !Formatting improved
00775 !
00776 !Revision 1.1  2004/04/06 07:56:38  elmo
00777 !Finishing AMG cycle. Makefile update.
00778 !
00779 !----------------------------------------------------------