DOUG 0.2

SpMtx_op_AB.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 !!Operation A*B for Sparse Matrix
00024 !! A, B - sparse Matrix
00025 !!----------------------------------------------------------
00026 module SpMtx_op_AB
00027   use realkind
00028   use SpMtx_class
00029   use SpMtx_arrangement
00030   use SpMtx_util
00031   Implicit None
00032   Contains
00033 !----------------------------------------------------------
00034 !Operation A*B
00035 !    Arguments:
00036 !           A,B - sparse matrix (type SpMtx)
00037 !    Result: sparse matrix A*B=AB
00038 !----------------------------------------------------------
00039   function SpMtx_AB(A,B,AT,BT) result(AB)
00040     implicit none
00041     type(SpMtx), intent(inout) :: A, B
00042     logical,intent(in),optional :: AT,BT !A-transp, B-transp?
00043     type(SpMtx)                :: AB
00044     integer                    :: i,j,jj,k,kk,kkk,nnz,nind,Annz,Bnnz
00045     integer                    :: aj,bi,as,ae,bs,be,ac,bc
00046     integer                    :: counter,maxrowa,maxrowb,maxrow
00047     integer,dimension(:),allocatable :: indeces,coladder,colindeces,ABi,ABj,ABb
00048     real(kind=rk),dimension(:),allocatable :: indvals,ABv
00049     integer,dimension(:),pointer     :: Aindi,Aindj,Bindi,Bindj
00050     integer                          :: Anrows,Ancols,Bnrows,Bncols
00051     integer :: Asparsity,Bsparsity,nnzest
00052     real :: rl
00053     logical :: cont
00054     ! Timing:
00055     !real(kind=rk) :: t1, t2
00056 
00057     if (sctls%verbose>3) then
00058       write(stream,"(A,': ',I0,'x',I0,', ',I0,'x',I0,A,L,',',L)") &
00059            "Multiplying matrices with shapes", &
00060            A%nrows, A%ncols, B%nrows, B%ncols, &
00061            ", transposed:",AT,BT
00062     end if
00063     !t1 = MPI_WTIME()
00064     if (present(AT).and.AT) then
00065       Aindi=>A%indj
00066       Aindj=>A%indi
00067       !Anrows=A%Ncols
00068       !Ancols=A%Nrows
00069       Anrows=max(A%ncols, maxval(A%indj))
00070       Ancols=max(A%nrows, maxval(A%indi))
00071       call SpMtx_arrange(A,D_SpMtx_ARRNG_COLS,sort=.false.,nnz=max(A%nnz,A%ol0nnz), &
00072                          nrows=Ancols,ncols=Anrows)        
00073     else
00074       Aindi=>A%indi
00075       Aindj=>A%indj
00076       Anrows=max(A%nrows, maxval(Aindi))
00077       Ancols=max(A%ncols, maxval(Aindj))
00078       call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS,sort=.false.,nnz=max(A%nnz,A%ol0nnz), &
00079                          nrows=Anrows,ncols=Ancols)        
00080     endif
00081     if (present(BT).and.BT) then
00082       Bindi=>B%indj
00083       Bindj=>B%indi
00084       !Bnrows=B%Ncols
00085       !Bncols=B%Nrows
00086       Bnrows=max(B%ncols, maxval(B%indj))
00087       Bncols=max(B%nrows, maxval(B%indi))
00088       call SpMtx_arrange(B,D_SpMtx_ARRNG_COLS,sort=.false.,nnz=max(B%nnz,B%ol0nnz), &
00089                          nrows=Bncols,ncols=Bnrows)
00090     else
00091       Bindi=>B%indi
00092       Bindj=>B%indj
00093       !Bnrows=B%Nrows
00094       !Bncols=B%Ncols
00095       Bnrows=max(B%nrows, maxval(Bindi))
00096       Bncols=max(B%ncols, maxval(Bindj))
00097       call SpMtx_arrange(B,D_SpMtx_ARRNG_ROWS,sort=.false.,nnz=max(B%nnz,B%ol0nnz), &
00098                          nrows=Bnrows,ncols=Bncols)        
00099     endif
00100     if (Ancols /= Bnrows) then
00101       write(stream,*)"ERROR: A*B is impossible!!!"
00102       write(stream,*)'Ancols=',Ancols,'Bnrows=',Bnrows
00103       call DOUG_Abort("[SpMtx_AB] failed", -1)
00104     endif
00105 !print *,'Anrows,Ancols,Bnrows,Bncols:',Anrows,Ancols,Bnrows,Bncols
00106     rl=1.0*Anrows*Ancols
00107     Annz=max(A%nnz,A%ol0nnz)
00108     Bnnz=max(B%nnz,B%ol0nnz)
00109     Asparsity=int(rl/Annz)
00110 
00111     rl=1.0*Bnrows*Bncols
00112     Bsparsity=int(rl/Bnnz)
00113     rl=1.0*Anrows*Bncols
00114     !nnzest=int(rl/min(Asparsity,Bsparsity)*4.0)
00115     nnzest=int(rl/max(1,min(Asparsity,Bsparsity))*4.0)+max(Annz,Bnnz)
00116     maxrowa=maxval(A%M_bound(2:Anrows+1)-A%M_bound(1:Anrows))
00117     maxrowb=maxval(B%M_bound(2:Bnrows+1)-B%M_bound(1:Bnrows))
00118     maxrow=maxrowa*maxrowb
00119     allocate(indeces(Bncols)) 
00120     allocate(coladder(Bncols)) 
00121     allocate(colindeces(maxrow)) 
00122     allocate(indvals(maxrow)) !could be smaller with some estimate...
00123     allocate(ABb(Anrows+1))
00124     cont=.true.
00125     do while(cont)
00126       cont=.false.
00127       allocate(ABi(nnzest))
00128       allocate(ABj(nnzest))
00129       allocate(ABv(nnzest))
00130       counter=1
00131       nnz = 0
00132       indeces=0 
00133       coladder=0 
00134   try:do i=1,Anrows
00135         ABb(i)=counter
00136         nind=0
00137         as=A%M_bound(i);ae=A%M_bound(i+1)-1
00138         do k=as,ae
00139           kk=Aindj(k)
00140           bs=B%M_bound(kk);be=B%M_bound(kk+1)-1
00141           do j=bs,be
00142             jj=Bindj(j)
00143             if (coladder(jj)/=i) then
00144               coladder(jj)=i
00145               nind=nind+1
00146               indeces(jj)=nind ! haven't added yet here
00147               colindeces(nind)=jj
00148               indvals(nind)=A%val(k)*B%val(j)
00149               nnz=nnz+1
00150      !if (i==1.or.jj==2) write(*,'((i4) (i4) (f15.10))')i,jj,indvals(nind)
00151               if (nnz>nnzest) then
00152                 print *,'nnz estimate too small in SpMtx_op_AB.f90'
00153                 !stop
00154                 deallocate(ABv)
00155                 deallocate(ABj)
00156                 deallocate(ABi)
00157                 !nnzest=nnzest+A%nnz+B%nnz
00158                 nnzest=nnzest+max(Annz,Bnnz)
00159                 print *,'incresing nnzest to',nnzest
00160                 cont=.true.
00161                 exit try
00162               endif
00163             else !(coladder(jj)==i) This value has already been added to:
00164               kkk=indeces(jj)
00165      !if (i==1.or.jj==2) write(*,'((i4) (i4) (f15.10)"+"(f15.10) (i4) )') &
00166      !           i,jj,indvals(nind),A%val(k)*B%val(j),kkk
00167               indvals(kkk)=indvals(kkk)+A%val(k)*B%val(j)
00168             endif
00169           enddo
00170         enddo
00171         ABi(counter:nnz)=i
00172         ABj(counter:nnz)=colindeces(1:nind)
00173         ABv(counter:nnz)=indvals(1:nind)
00174         counter=counter+nind
00175       enddo try
00176     enddo
00177     !constructing new sparse matrix
00178     ABb(Anrows+1)=counter
00179     AB = SpMtx_newInit(nnz=nnz,             &
00180                 nblocks=1,                  &
00181                   nrows=Anrows,             &
00182                   ncols=Bncols,             &
00183              symmstruct=.false.,            &
00184             symmnumeric=.false.,            &
00185                    indi=ABi,                &
00186                    indj=ABj,                &
00187                     val=ABv,                &
00188            arrange_type=D_SpMtx_ARRNG_ROWS, &
00189                 M_bound=ABb)
00190 !write(stream,*)'AB:',AB%nrows,AB%ncols,AB%nnz
00191     deallocate(ABv)
00192     deallocate(ABj)
00193     deallocate(ABi)
00194     !print *,'nnz est. and actual in AB:',nnzest,nnz
00195     deallocate(ABb)
00196     deallocate(indvals) 
00197     deallocate(colindeces) 
00198     deallocate(coladder) 
00199     deallocate(indeces) 
00200   end function SpMtx_AB
00201 
00202   ! with memorising nonzeroe matches:
00203   function SpMtx_AB_nonopt(A,B,AT,BT) result(AB)
00204     implicit none
00205     type(SpMtx), intent(inout) :: A, B
00206     logical,intent(in),optional :: AT,BT !A-transp, B-transp?
00207     type(SpMtx)                :: AB
00208     integer                    :: i,j,k,nnz,counter
00209     integer                    :: aj,bi,as,ae,bs,be,ac,bc
00210     logical                    :: cont,nnz_not_added,nops_ok,enlarge_nops
00211     integer                    :: nops,done
00212     integer,dimension(:),allocatable :: imem,jmem,amem,bmem,zmem
00213     integer,dimension(:),pointer     :: Aindi,Aindj,Bindi,Bindj
00214     integer                          :: Anrows,Ancols,Bnrows,Bncols
00215     real :: ir,jr,nr
00216     ! Timing:
00217     real(kind=rk) :: t1, t2
00218 
00219     t1 = MPI_WTIME()
00220     if (present(AT).and.AT) then
00221       Aindi=>A%indj
00222       Aindj=>A%indi
00223       Anrows=A%Ncols
00224       Ancols=A%Nrows
00225       call SpMtx_arrange(A,D_SpMtx_ARRNG_COLS)        
00226     else
00227       Aindi=>A%indi
00228       Aindj=>A%indj
00229       Anrows=A%Nrows
00230       Ancols=A%Ncols
00231       call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS)        
00232     endif
00233     if (present(BT).and.BT) then
00234       Bindi=>B%indj
00235       Bindj=>B%indi
00236       Bnrows=B%Ncols
00237       Bncols=B%Nrows
00238       call SpMtx_arrange(B,D_SpMtx_ARRNG_ROWS)        
00239     else
00240       Bindi=>B%indi
00241       Bindj=>B%indj
00242       Bnrows=B%Nrows
00243       Bncols=B%Ncols
00244       call SpMtx_arrange(B,D_SpMtx_ARRNG_COLS)        
00245     endif
00246 
00247     if (Ancols /= Bnrows) then
00248       print*, "ERROR: A*B is impossible!!!"
00249       stop
00250     endif
00251     write(*,*) 'AB startup-time:',MPI_WTIME()-t1
00252     ! we are counting nonzero elements, but at the same time, we remember
00253     !   where some multiplications need to be done later...:
00254     ! need to estimate the number of multops
00255     !   find average number of nonzeroes in each A row and B column:
00256     i=sum(A%M_bound(2:Anrows)-A%M_bound(1:Anrows-1)) ; ir=1.0*i/Anrows
00257     !print *,'i,ir:',i,ir
00258     j=sum(B%M_bound(2:Bncols)-B%M_bound(1:Bncols-1)) ; jr=1.0*j/Bncols
00259     !print *,'j,jr:',j,jr
00260     ! probability, that an arbitrary index falls onto nonzero in matrix B col:
00261     jr=1.0*jr/Bnrows
00262     !print *,'jr=',jr
00263     ! => probable number of matching pairs for each row of A is:
00264     nr=ir*jr
00265     !print *,'nr=',nr
00266     ! ...and we have that many B columns per each row and that many A rows:
00267     nops=int(nr*Bncols*Anrows)
00268     !print *,'real:',nr*Bncols*Anrows,nr,Bncols*Anrows
00269     !print *,'int:',int(nr*Bncols*Anrows)
00270     !print *,'@@@@@@@@@',nops,nr,Bncols*Anrows,ir,i
00271     !Let's make it bigger, just in case:
00272     nops=int(nops*1.1)
00273     ! now, these were barely estimates, for 100% success, we allow resizing:
00274     nops_ok=.false.
00275     do while(.not.nops_ok)
00276       nnz = 0
00277       counter=0
00278       enlarge_nops=.false.
00279       allocate(imem(nops))
00280       allocate(jmem(nops))
00281       allocate(amem(nops))
00282       allocate(bmem(nops))
00283       allocate(zmem(nops))
00284 outer:do i=1,Anrows
00285         as=A%M_bound(i);ae=A%M_bound(i+1)-1
00286         do j=1,Bncols
00287           bs=B%M_bound(j);be=B%M_bound(j+1)-1
00288           cont=.true.
00289           nnz_not_added=.true.
00290           ac=as;bc=bs
00291           aj=Aindj(ac);bi=Bindi(bc)
00292           do while(cont)
00293             if (aj==bi) then
00294               if (nnz_not_added) then
00295                 nnz=nnz+1
00296                 nnz_not_added=.false.
00297               endif
00298               counter=counter+1
00299               if (counter>nops) then
00300                 enlarge_nops=.true.
00301                 cont=.false.
00302                 done=(i-1)*Anrows+j
00303                 !j=Bncols+1 ! (for ending the for-cycles...)
00304                 !i=Anrows
00305                 exit outer
00306               else
00307                 imem(counter)=i
00308                 jmem(counter)=j
00309                 amem(counter)=ac
00310                 bmem(counter)=bc
00311                 zmem(counter)=nnz
00312               endif
00313               if (ac<ae) then
00314                 ac=ac+1
00315                 aj=Aindj(ac)
00316               else
00317                 cont=.false.
00318               endif
00319               if (bc<be) then
00320                 bc=bc+1
00321                 bi=Bindi(bc)
00322               else
00323                 cont=.false.
00324               endif
00325             elseif (aj>bi) then ! advance in B:
00326               if (bc<be) then
00327                 bc=bc+1
00328                 bi=Bindi(bc)
00329               else
00330                 cont=.false.
00331               endif
00332             else ! (bi>aj) ! advance in A:
00333               if (ac<ae) then
00334                 ac=ac+1
00335                 aj=Aindj(ac)
00336               else
00337                 cont=.false.
00338               endif
00339             endif
00340           enddo
00341         enddo
00342       enddo outer
00343       if (enlarge_nops) then ! the estimation failed...
00344         deallocate(zmem)
00345         deallocate(bmem)
00346         deallocate(amem)
00347         deallocate(jmem)
00348         deallocate(imem)
00349         k=nops
00350         nops=nops+nops*(done+nops/10)/(Anrows*Bncols)
00351         print *,'Warning: estimated nops increased from ',k,' to', &
00352              nops,done,Anrows,Bncols,i,j,counter
00353         stop
00354       else
00355         nops_ok=.true.
00356       endif
00357     enddo
00358     print *,'Estimated nops was:',nops,' #ops is:',counter
00359     print *,'number of nonzeroes in AB is:',nnz
00360     !constructing new sparse matrix
00361     AB = SpMtx_newInit(nnz)
00362     ! Initialising AB%nrows and AB%ncols
00363     AB%nrows=Anrows; AB%ncols=Bncols
00364     AB%indi=0
00365     do i=1,counter ! now the real multiplication:
00366       k=zmem(i)
00367       if (AB%indi(k)==0) then
00368         AB%indi(k)=imem(i)
00369         AB%indj(k)=jmem(i)
00370         AB%val(k)=A%val(amem(i))*B%val(bmem(i))
00371       else
00372         AB%val(k)=AB%val(k)+A%val(amem(i))*B%val(bmem(i))
00373       endif
00374     enddo
00375     deallocate(zmem)
00376     deallocate(bmem)
00377     deallocate(amem)
00378     deallocate(jmem)
00379     deallocate(imem)
00380   end function SpMtx_AB_nonopt
00381 
00382   function SpMtx_AB2(A,B,AT,BT) result(AB)
00383     implicit none
00384     type(SpMtx), intent(inout) :: A, B
00385     logical,intent(in),optional :: AT,BT !A-transp, B-transp?
00386     type(SpMtx)                :: AB
00387     integer                    :: i,j,k,nnz,counter
00388     integer                    :: aj,bi,as,ae,bs,be,ac,bc
00389     logical                    :: cont,nnz_not_added,nops_ok,enlarge_nops
00390     integer                    :: nops,done
00391     integer,dimension(:),pointer     :: Aindi,Aindj,Bindi,Bindj
00392     integer                          :: Anrows,Ancols,Bnrows,Bncols
00393     real :: ir,jr,nr
00394     ! Timing:
00395     real(kind=rk) :: t1, t2
00396 
00397     t1 = MPI_WTIME()
00398     if (present(AT).and.AT) then
00399       Aindi=>A%indj
00400       Aindj=>A%indi
00401       Anrows=A%Ncols
00402       Ancols=A%Nrows
00403       call SpMtx_arrange(A,D_SpMtx_ARRNG_COLS,.true.)        
00404     else
00405       Aindi=>A%indi
00406       Aindj=>A%indj
00407       Anrows=A%Nrows
00408       Ancols=A%Ncols
00409       call SpMtx_arrange(A,D_SpMtx_ARRNG_ROWS,.true.)        
00410     endif
00411     if (present(BT).and.BT) then
00412       Bindi=>B%indj
00413       Bindj=>B%indi
00414       Bnrows=B%Ncols
00415       Bncols=B%Nrows
00416       call SpMtx_arrange(B,D_SpMtx_ARRNG_ROWS,.true.)        
00417     else
00418       Bindi=>B%indi
00419       Bindj=>B%indj
00420       Bnrows=B%Nrows
00421       Bncols=B%Ncols
00422       call SpMtx_arrange(B,D_SpMtx_ARRNG_COLS,.true.)        
00423     endif
00424 
00425     if (Ancols /= Bnrows) then
00426       print*, "ERROR: A*B is impossible!!!"
00427       stop
00428     endif
00429     write(*,*) 'AB startup-time:',MPI_WTIME()-t1
00430     ! we are counting nonzero elements    do while(.not.nops_ok)
00431     nnz = 0
00432     do i=1,Anrows
00433       as=A%M_bound(i);ae=A%M_bound(i+1)-1
00434       if(ae<as) cycle
00435       do j=1,Bncols
00436         bs=B%M_bound(j);be=B%M_bound(j+1)-1
00437         if (be<bs) cycle
00438         nnz_not_added=.true.
00439         ac=as;bc=bs
00440         aj=Aindj(ac);bi=Bindi(bc)
00441         cont=.TRUE.
00442         do while(cont)
00443           if (aj==bi) then
00444             if (nnz_not_added) then
00445               nnz=nnz+1
00446               nnz_not_added=.false.
00447             endif
00448             if (ac<ae) then
00449               ac=ac+1
00450               aj=Aindj(ac)
00451             else
00452               cont=.false.
00453             endif
00454             if (bc<be) then
00455               bc=bc+1
00456               bi=Bindi(bc)
00457             else
00458               cont=.false.
00459             endif
00460           elseif (aj>bi) then ! advance in B:
00461             if (bc<be) then
00462               bc=bc+1
00463               bi=Bindi(bc)
00464             else
00465               cont=.false.
00466             endif
00467           else ! (bi>aj) ! advance in A:
00468             if (ac<ae) then
00469               ac=ac+1
00470               aj=Aindj(ac)
00471             else
00472               cont=.false.
00473             endif
00474           endif
00475         enddo
00476       enddo
00477     enddo
00478     print *,'number of nonzeroes in AB is:',nnz
00479     !constructing new sparse matrix
00480     AB = SpMtx_newInit(nnz)
00481     ! Initialising AB%nrows and AB%ncols
00482     AB%nrows=Anrows; AB%ncols=Bncols
00483     !2nd pass starts here:
00484     nnz = 0
00485     do i=1,Anrows
00486       as=A%M_bound(i);ae=A%M_bound(i+1)-1
00487       if(ae<as) cycle
00488       do j=1,Bncols
00489         bs=B%M_bound(j);be=B%M_bound(j+1)-1
00490         if (be<bs) cycle
00491         nnz_not_added=.true.
00492         ac=as;bc=bs
00493         aj=Aindj(ac);bi=Bindi(bc)
00494         cont=.TRUE.
00495         do while(cont)
00496           if (aj==bi) then
00497             if (nnz_not_added) then
00498               nnz=nnz+1
00499               nnz_not_added=.false.
00500               AB%indi(nnz)=i
00501               AB%indj(nnz)=j
00502               AB%val(nnz)=A%val(ac)*B%val(bc)
00503             else
00504               AB%val(nnz)=AB%val(nnz)+A%val(ac)*B%val(bc)
00505             endif
00506             if (ac<ae) then
00507               ac=ac+1
00508               aj=Aindj(ac)
00509             else
00510               cont=.false.
00511             endif
00512             if (bc<be) then
00513               bc=bc+1
00514               bi=Bindi(bc)
00515             else
00516               cont=.false.
00517             endif
00518           elseif (aj>bi) then ! advance in B:
00519             if (bc<be) then
00520               bc=bc+1
00521               bi=Bindi(bc)
00522             else
00523               cont=.false.
00524             endif
00525           else ! (bi>aj) ! advance in A:
00526             if (ac<ae) then
00527               ac=ac+1
00528               aj=Aindj(ac)
00529             else
00530               cont=.false.
00531             endif
00532           endif
00533         enddo
00534       enddo
00535     enddo
00536   end function SpMtx_AB2
00537 
00538 end module SpMtx_op_AB
00539