|
DOUG 0.2
|
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
1.7.3-20110217