|
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 Module SpMtx_arrangement 00023 !!-------------------------------------------------------- 00024 !!Arrange elements in sparse matrix 00025 !!-------------------------------------------------------- 00026 00027 use RealKind 00028 use SpMtx_class 00029 use SpMtx_util 00030 use Mesh_class 00031 use globals 00032 use Aggregate_mod 00033 00034 Implicit None 00035 00036 #include<doug_config.h> 00037 00038 ! "on-the-fly" real/complex picking 00039 #ifdef D_COMPLEX 00040 #define float complex 00041 #else 00042 #define float real 00043 #endif 00044 00045 logical, parameter :: arr_by_j = .true., arr_by_i = .false. 00046 00047 CONTAINS 00048 00049 include 'SpMtx_arrange_qs.F90' 00050 00051 !---------------------------------------------------------- 00052 ! arrange for ROWs or COLUMNs 00053 ! Arguments: 00054 ! M - sparse matrix 00055 ! optional: 00056 ! arrange_type -- D_SpMtx_ARRNG_ROWS or D_SpMtx_ARRNG_COLS 00057 !---------------------------------------------------------- 00058 subroutine SpMtx_arrange(M,arrange_type,sort,nnz,nrows,ncols) 00059 Implicit None 00060 Type(SpMtx), intent(in out) :: M !Initial matrix 00061 integer, optional, intent(in) :: arrange_type ! how to arrange? 00062 logical, optional, intent(in) :: sort ! wheather entries ascending? 00063 integer, optional, intent(in) :: nnz ! if the M%nnz to be overridden 00064 integer, optional, intent(in) :: nrows ! if the M%nrows to be overridden 00065 integer, optional, intent(in) :: ncols ! if the M%ncols to be overridden 00066 logical :: columnwise,dosort,ok 00067 integer :: i,j,k,kk,ind_beg,ind,Mnnz,Mnrows,Mncols 00068 integer, dimension(:), allocatable :: el !elements vector 00069 integer, dimension(:), allocatable :: indi,indj !helper vectors 00070 float(kind=rk),dimension(:),allocatable :: val 00071 integer, dimension(:), allocatable :: sortref !helper for sorting 00072 integer :: at 00073 !- - - - - - - - - - - - - - - - - - - - - - - - 00074 ok=.true. 00075 if (present(arrange_type)) then 00076 at=arrange_type 00077 else 00078 at=D_SpMtx_ARRNG_ROWS 00079 endif 00080 if (M%arrange_type/=at) then 00081 ok=.false. 00082 endif 00083 if (present(nnz)) then 00084 Mnnz=nnz 00085 else 00086 Mnnz=M%nnz 00087 endif 00088 if (present(nrows)) then 00089 Mnrows=nrows 00090 else 00091 Mnrows=M%nrows 00092 endif 00093 if (present(ncols)) then 00094 Mncols=ncols 00095 else 00096 Mncols=M%ncols 00097 endif 00098 dosort=.false. 00099 if (present(sort)) then 00100 if (sort) then 00101 dosort=sort 00102 endif 00103 endif 00104 if (at==D_SpMtx_ARRNG_ROWS) then 00105 columnwise =.false. 00106 if (ok) then 00107 if (size(M%M_bound)/=Mnrows+1) then 00108 ok=.false. 00109 endif 00110 endif 00111 elseif (at==D_SpMtx_ARRNG_COLS) then 00112 columnwise = .true. 00113 if (ok) then 00114 if (size(M%M_bound)/=Mncols+1) then 00115 ok=.false. 00116 endif 00117 endif 00118 else 00119 write(stream,*)'WARNING: SpMtx_arrange to ',at,' not done' 00120 return 00121 endif 00122 if (ok) then 00123 write(stream,*)'WARNING: SpMtx_arrange to ',at,' not needed...' 00124 return 00125 endif 00126 M%arrange_type=at 00127 00128 !===== allocate memory and control arrange_type 00129 if (columnwise) then !!!columns 00130 allocate(el(Mncols)) 00131 allocate(M%M_bound(Mncols+1)) 00132 else !!!rows 00133 allocate(el(Mnrows)) 00134 allocate(M%M_bound(Mnrows+1)) 00135 end if 00136 00137 00138 !===== 1.find how many elements are every row/col 00139 if (.not.present(nnz).and.associated(M%mtx_bbe)) then 00140 if (M%mtx_bbe(2,2)>0) then 00141 Mnnz=M%mtx_bbe(2,2) 00142 endif 00143 endif 00144 el = 0 00145 do i = 1, Mnnz 00146 if (columnwise) then 00147 k=M%indj(i) 00148 else 00149 k=M%indi(i) 00150 end if 00151 el(k)=el(k)+1 00152 end do 00153 !===== generate M_bound vector 00154 if (size(M%M_bound)>0) M%M_bound(1) = 1 00155 do i = 1, size(M%M_bound) - 1 00156 M%M_bound(i+1) = M%M_bound(i) + el(i) 00157 end do 00158 allocate(indi(Mnnz),indj(Mnnz),val(Mnnz)) 00159 if (dosort) then 00160 ! 2.find the order: 00161 allocate(sortref(Mnnz)) 00162 el = 0 00163 if (columnwise) then 00164 do i = 1,Mnnz 00165 k=M%indj(i) 00166 ind_beg = M%M_bound(k) 00167 if (el(k)==0) then ! the first element 00168 sortref(ind_beg)=i 00169 el(k)=1 00170 else ! rank and insert into the sorted list: 00171 ind = M%M_bound(k)+el(k) 00172 j=ind_beg 00173 whilc:do while (.true.) 00174 if (j==ind) then ! adding to the end: 00175 sortref(j)=i 00176 el(k)=el(k)+1 00177 exit whilc 00178 elseif (M%indi(i)<M%indi(sortref(j))) then ! add betw.: 00179 do kk=ind,j+1,-1 ! advance the rest of the sequence by 1 step: 00180 sortref(kk)=sortref(kk-1) 00181 enddo 00182 sortref(j)=i 00183 el(k)=el(k)+1 00184 exit whilc 00185 else !advance to the next sorted element: 00186 j=j+1 00187 endif 00188 enddo whilc 00189 endif 00190 enddo 00191 else ! rowwise: 00192 do i = 1,Mnnz 00193 k=M%indi(i) 00194 ind_beg = M%M_bound(k) 00195 if (el(k)==0) then ! the first element 00196 sortref(ind_beg)=i 00197 el(k)=1 00198 else ! rank and insert into the sorted list: 00199 ind = M%M_bound(k)+el(k) 00200 j=ind_beg 00201 whilr:do while (.true.) 00202 if (j==ind) then ! adding to the end: 00203 sortref(j)=i 00204 el(k)=el(k)+1 00205 exit whilr 00206 elseif (M%indj(i)<M%indj(sortref(j))) then ! add betw.: 00207 do kk=ind,j+1,-1 ! advance the rest of the sequence by 1 step: 00208 sortref(kk)=sortref(kk-1) 00209 enddo 00210 sortref(j)=i 00211 el(k)=el(k)+1 00212 exit whilr 00213 else !advance to the next sorted element: 00214 j=j+1 00215 endif 00216 enddo whilr 00217 endif 00218 enddo 00219 endif 00220 ! 3.rearrange arrays: 00221 do i = 1,Mnnz 00222 ind=sortref(i) ! where to get the values for this pos 00223 indi(i) = M%indi(ind) 00224 indj(i) = M%indj(ind) 00225 val(i) = M%val(ind) 00226 enddo 00227 deallocate(sortref) 00228 else 00229 el = 0 00230 do i = 1,Mnnz 00231 if (columnwise) then 00232 k=M%indj(i) 00233 else 00234 k=M%indi(i) 00235 end if 00236 ind = M%M_bound(k)+el(k) 00237 el(k)=el(k)+1 00238 indi(ind) = M%indi(i) 00239 indj(ind) = M%indj(i) 00240 val(ind) = M%val(i) 00241 end do 00242 endif 00243 !===== finishing ... 00244 M%indi(1:Mnnz) = indi 00245 M%indj(1:Mnnz) = indj 00246 M%val(1:Mnnz) = val 00247 deallocate(val,indj,indi) 00248 deallocate(el) 00249 00250 ! destroy strong 00251 if (associated(M%strong)) deallocate(M%strong) 00252 end subroutine SpMtx_arrange 00253 00254 !------------------------------------------------------ 00255 ! Matrix consolidation: 00256 ! find duplicate elements and add them together 00257 !------------------------------------------------------ 00258 subroutine SpMtx_consolidate(M,add) 00259 use RealKind 00260 use SpMtx_class 00261 use Mesh_class 00262 Implicit None 00263 Type(SpMtx), intent(inout) :: M !Initial matrix 00264 logical :: add 00265 integer :: i, k 00266 00267 ! Sort the matrix by indices (I pray it works for duplicates too) 00268 call SpMtx_arrange(M,sort=.true.) 00269 ! Consolidate it in one pass 00270 k=1; 00271 do i=2,M%nnz 00272 if (M%indi(i)==M%indi(k) .and. M%indj(i)==M%indj(k)) then 00273 if (add) then 00274 M%val(k)=M%val(k)+M%val(i); 00275 endif 00276 else 00277 k=k+1; 00278 if (k/=i) then 00279 M%indi(k)=M%indi(i) 00280 M%indj(k)=M%indj(i) 00281 M%val(k)=M%val(i) 00282 endif 00283 endif 00284 enddo 00285 00286 ! Remove arrangement indications (they could remain valid with extra work) 00287 M%arrange_type=D_SpMtx_ARRNG_NO 00288 deallocate(M%M_bound) 00289 00290 ! And resize it 00291 call SpMtx_resize(M,k) 00292 end subroutine SpMtx_consolidate 00293 00294 !------------------------------------------------------ 00295 ! Diagonal scaling of matrix: 00296 ! diagonal value is stored in diag 00297 !------------------------------------------------------ 00298 subroutine SpMtx_scale(M,A_ghost) 00299 Implicit None 00300 Type(SpMtx), intent(in out) :: M 00301 Type(SpMtx),intent(in out),optional :: A_ghost 00302 integer :: i,j,ndiags 00303 float(kind=rk), dimension(:), pointer :: scalerval 00304 if (M%scaling==D_SpMtx_SCALE_NO.or.M%scaling==D_SpMtx_SCALE_UNDEF) then 00305 !ndiags=min(M%nrows,M%ncols) 00306 ndiags=max(M%nrows,M%ncols) 00307 if (.not.associated(M%diag)) then 00308 allocate(M%diag(ndiags)) 00309 endif 00310 allocate(scalerval(ndiags)) 00311 do i=1,M%mtx_bbe(2,2) 00312 if (M%indi(i)==M%indj(i)) then 00313 M%diag(M%indi(i))=M%val(i) 00314 endif 00315 enddo 00316 if (present(A_ghost).and.associated(A_ghost%indi)) then 00317 do i=1,A_ghost%nnz 00318 if (A_ghost%indi(i)==A_ghost%indj(i)) then 00319 M%diag(A_ghost%indi(i))=A_ghost%val(i) 00320 endif 00321 enddo 00322 endif 00323 if (M%symmstruct.and.M%symmnumeric) then ! the symmetric case: 00324 do i=1,ndiags 00325 scalerval(i)=dsqrt(dabs(M%diag(i))) 00326 enddo 00327 do i=1,M%mtx_bbe(2,2) 00328 M%val(i)=M%val(i)/scalerval(M%indi(i)) 00329 M%val(i)=M%val(i)/scalerval(M%indj(i)) 00330 enddo 00331 else ! unsymmetric case: 00332 do i=1,ndiags 00333 scalerval(i)=dabs(M%diag(i)) 00334 enddo 00335 do i=1,M%nnz 00336 M%val(i)=M%val(i)/scalerval(M%indi(i)) 00337 enddo 00338 endif 00339 deallocate(scalerval) 00340 M%scaling=D_SpMtx_SCALE_DIAG 00341 else 00342 if (D_MSGLVL>1) then 00343 write(stream,*) 'WARNING: matrix already in scaled format:',M%scaling 00344 endif 00345 endif 00346 end subroutine SpMtx_scale 00347 !------------------------------------------------------ 00348 ! Unscaling of matrix with that stored in diag 00349 !------------------------------------------------------ 00350 subroutine SpMtx_unscale(M) 00351 Implicit None 00352 Type(SpMtx), intent(in out) :: M 00353 integer :: i,ndiags,nnz 00354 float(kind=rk), dimension(:), pointer :: scalerval 00355 if (M%mtx_bbe(2,2)>0) then 00356 nnz=M%mtx_bbe(2,2) 00357 else 00358 nnz=M%nnz 00359 endif 00360 if (M%scaling==D_SpMtx_SCALE_DIAG.or. & 00361 M%scaling==D_SpMtx_SCALE_DIAG_FILTERED) then 00362 ndiags=max(M%nrows,M%ncols) 00363 allocate(scalerval(ndiags)) 00364 if (M%symmstruct.and.M%symmnumeric) then ! the symmetric case: 00365 do i=1,ndiags 00366 scalerval(i)=dsqrt(dabs(M%diag(i))) 00367 enddo 00368 do i=1,nnz 00369 M%val(i)=M%val(i)*scalerval(M%indi(i)) 00370 M%val(i)=M%val(i)*scalerval(M%indj(i)) 00371 enddo 00372 else ! unsymmetric case: 00373 do i=1,ndiags 00374 scalerval(i)=dabs(M%diag(i)) 00375 enddo 00376 do i=1,nnz 00377 M%val(i)=M%val(i)*scalerval(M%indi(i)) 00378 enddo 00379 endif 00380 deallocate(scalerval) 00381 deallocate(M%diag) 00382 M%scaling=D_SpMtx_SCALE_NO 00383 else 00384 if (D_MSGLVL>1) then 00385 write(stream,*) 'WARNING: matrix already in unscaled format:',M%scaling 00386 endif 00387 endif 00388 end subroutine SpMtx_unscale 00389 00390 !---------------------------------------------------------- 00391 ! build matrix reference arrays: 00392 ! for quick strong row reference: 00393 ! rowstart(1:n+1), colnrs(1:noffdels) 00394 ! colstart(1:m+1), rownrs(1:noffdels) 00395 !---------------------------------------------------------- 00396 subroutine SpMtx_build_refs(M,noffdels,rowstart,colnrs,colstart,rownrs) 00397 Implicit None 00398 Type(SpMtx), intent(in out) :: M 00399 integer, intent(out) :: noffdels 00400 integer, dimension(:), pointer :: rowstart,colnrs,colstart,rownrs 00401 integer :: i,nn,mm,nnz,ind,ii,jj 00402 integer, dimension(:), allocatable :: nnz_in_row ! #nonzeroes in each row 00403 integer, dimension(:), allocatable :: nnz_in_col ! #nonzeroes in each col 00404 !- - - - - - - - - - - - - - - - - - - - - - - - 00405 nn=M%nrows 00406 mm=M%ncols 00407 nnz=M%nnz 00408 !===== allocate memory 00409 allocate(nnz_in_row(nn)) 00410 allocate(nnz_in_col(mm)) 00411 !===== find how many elements are every row/col 00412 nnz_in_row=0 00413 nnz_in_col=0 00414 do i = 1,nnz 00415 if (M%strong(i)) then 00416 ii=M%indi(i) 00417 jj=M%indj(i) 00418 if (ii/=jj) then 00419 nnz_in_row(ii) = nnz_in_row(ii) + 1 00420 nnz_in_col(jj) = nnz_in_col(jj) + 1 00421 endif 00422 endif 00423 end do 00424 allocate(rowstart(nn+1)) ! do not forget to deallocate somewhere! 00425 allocate(colstart(mm+1)) 00426 rowstart(1)=1 00427 do i=1,nn 00428 rowstart(i+1)=rowstart(i)+nnz_in_row(i) 00429 enddo 00430 colstart(1)=1 00431 do i=1,mm 00432 colstart(i+1)=colstart(i)+nnz_in_col(i) 00433 enddo 00434 noffdels=colstart(mm+1)-1 00435 allocate(colnrs(noffdels)) 00436 allocate(rownrs(noffdels)) 00437 nnz_in_row=0 00438 nnz_in_col=0 00439 do i = 1,nnz 00440 if (M%strong(i)) then 00441 ii=M%indi(i) 00442 jj=M%indj(i) 00443 if (ii/=jj) then 00444 colnrs(rowstart(ii)+nnz_in_row(ii))=jj 00445 nnz_in_row(ii) = nnz_in_row(ii) + 1 00446 rownrs(colstart(jj)+nnz_in_col(jj))=ii 00447 nnz_in_col(jj) = nnz_in_col(jj) + 1 00448 endif 00449 endif 00450 end do 00451 !===== finising ... 00452 deallocate(nnz_in_row) 00453 deallocate(nnz_in_col) 00454 end subroutine SpMtx_build_refs 00455 00456 !---------------------------------------------------------- 00461 subroutine SpMtx_build_refs_symm(A,noffdels,rowstart,colnrs,& 00462 sortdown,diagstrong) 00463 Implicit None 00464 Type(SpMtx), intent(in out) :: A 00465 integer, intent(out) :: noffdels 00466 integer, dimension(:), pointer :: rowstart,colnrs 00467 logical,intent(in),optional :: sortdown ! wheather dominant connections first? 00468 logical, dimension(:), pointer, optional :: diagstrong 00469 integer :: i,j,k,nn,nnz,ind,ind_beg,ii,jj 00470 integer, dimension(:), allocatable :: nnz_in_row ! #nonzeroes in each row 00471 integer, dimension(:), allocatable :: sortref ! for getting the sorted list 00472 logical :: diagonals 00473 !- - - - - - - - - - - - - - - - - - - - - - - - 00474 nn=A%nrows 00475 if (A%mtx_bbe(2,2)>0) then 00476 nnz=A%mtx_bbe(2,2) 00477 else 00478 nnz=A%nnz 00479 endif 00480 if (present(diagstrong)) then 00481 diagonals=.true. 00482 diagstrong=.false. 00483 else 00484 diagonals=.false. 00485 endif 00486 !===== allocate memory 00487 allocate(nnz_in_row(nn)) 00488 !===== find how many elements are every row/col 00489 nnz_in_row=0 00490 do i = 1,nnz 00491 if (A%strong(i)) then 00492 ii=A%indi(i) 00493 jj=A%indj(i) 00494 if (ii/=jj) then 00495 if (jj<=nn) then 00496 nnz_in_row(ii) = nnz_in_row(ii) + 1 00497 endif 00498 elseif(diagonals) then 00499 diagstrong(ii)=.true. 00500 endif 00501 endif 00502 end do 00503 allocate(rowstart(nn+1)) ! do not forget to deallocate somewhere! 00504 rowstart(1)=1 00505 do i=1,nn 00506 rowstart(i+1)=rowstart(i)+nnz_in_row(i) 00507 enddo 00508 noffdels=rowstart(nn+1)-1 00509 allocate(colnrs(noffdels)) 00510 nnz_in_row=0 00511 if (present(sortdown).and.sortdown) then 00512 allocate(sortref(noffdels)) 00513 do i = 1,nnz 00514 if (A%strong(i)) then 00515 ii=A%indi(i) 00516 jj=A%indj(i) 00517 if (ii/=jj) then 00518 if (jj<=nn) then 00519 ind_beg=rowstart(ii) 00520 if (nnz_in_row(ii)==0) then ! the first element 00521 sortref(ind_beg)=i 00522 nnz_in_row(ii)=1 00523 else ! rank and insert into the sorted list: 00524 ind = rowstart(ii)+nnz_in_row(ii)-1 ! last existing rowentry 00525 j=ind_beg 00526 whilr:do while (.true.) 00527 if (j==ind+1) then ! adding to the end: 00528 sortref(j)=i 00529 nnz_in_row(ii)=nnz_in_row(ii)+1 00530 exit whilr 00531 elseif (abs(A%val(i))>abs(A%val(sortref(j)))) then ! add betw.: 00532 do k=ind+1,j+1,-1 ! advance the rest of the sequence by 1 step: 00533 sortref(k)=sortref(k-1) 00534 enddo 00535 sortref(j)=i 00536 nnz_in_row(ii)=nnz_in_row(ii)+1 00537 exit whilr 00538 else !advance to the next sorted element: 00539 j=j+1 00540 endif 00541 enddo whilr 00542 endif 00543 endif 00544 endif 00545 endif 00546 end do 00547 ! set colnrs in sorted order: 00548 do i=1,noffdels 00549 colnrs(i)=A%indj(sortref(i)) 00550 enddo 00551 deallocate(sortref) 00552 else 00553 do i = 1,nnz 00554 if (A%strong(i)) then 00555 ii=A%indi(i) 00556 jj=A%indj(i) 00557 if (ii/=jj) then 00558 if (jj<=nn) then 00559 colnrs(rowstart(ii)+nnz_in_row(ii))=jj 00560 nnz_in_row(ii) = nnz_in_row(ii) + 1 00561 endif 00562 endif 00563 endif 00564 end do 00565 endif 00566 !===== finishing ... 00567 deallocate(nnz_in_row) 00568 end subroutine SpMtx_build_refs_symm 00569 00570 00571 !---------------------------------------------------------- 00574 subroutine SpMtx_roughly_aggregate(A,aggr,neighood,maxaggrsize,alpha) 00575 use globals 00576 use Mesh_class 00577 Implicit None 00578 Type(SpMtx),intent(in out) :: A ! our matrix 00579 type(AggrInfo),intent(out) :: aggr !< aggregates 00580 integer,intent(in) :: neighood ! 1-neighood,2-neighood or r-neighood... 00581 integer,intent(in) :: maxaggrsize 00582 float(kind=rk),intent(in) :: alpha 00583 integer :: wavelen 00584 integer :: i,j,j1,j2,jj,k,kk,n 00585 integer,dimension(:),pointer :: ii,stat 00586 integer,dimension(neighood) :: seeds 00587 integer,dimension(maxaggrsize,neighood) :: layer 00588 integer,dimension(neighood) :: layersize,aggrsize 00589 integer,dimension(2*maxaggrsize) :: newlayer 00590 integer :: newlayersize 00591 integer,dimension(neighood) :: actanum ! active aggregate number 00592 integer :: nact ! #active 00593 integer :: counter,rs,re,layernode,cn,noffdels,weakend 00594 integer :: nagrs ! number of aggregates 00595 logical :: overwrite=.false. 00596 logical :: rounding=.true. 00597 logical, dimension(:), pointer :: diagstrong 00598 00599 aggr = AggrInfo_New() 00600 00601 n=A%nrows 00602 allocate(ii(n)) 00603 allocate(aggr%inner%num(n)) 00604 aggr%inner%num(1:n)=0 00605 ii=(/ (i,i=1,n) /) !built-in array initialisation to 1:n 00606 call random_permutation(n,ii) 00607 !write(stream,*)'permutation is:',ii 00608 !write(stream,*)i,ii(i) 00609 00610 ! Start aggregating taking seeds from ii 00611 ! if node j aggregated set ii(j)=-abs(ii(j)) 00612 allocate(stat(n)) 00613 allocate(diagstrong(n)) 00614 stat=0 00615 layersize=0 00616 aggr%inner%nagr=0 00617 nact=0 00618 counter=0 00619 call SpMtx_build_refs_symm(A,noffdels, & 00620 A%strong_rowstart,A%strong_colnrs, & 00621 sortdown=.true.,diagstrong=diagstrong) 00622 aggrsize=0 00623 do i=1,n 00624 ! find next unaggregated seed in permutation order: 00625 if (stat(ii(i))==0) then 00626 !if weak, put exchange it with a strong one from the end 00627 if (.not.diagstrong(ii(i))) then 00628 if (weakend>i) then 00629 do while(.true.) 00630 weakend=weakend-1 00631 if (weakend<=i) then 00632 exit 00633 endif 00634 if (stat(ii(weakend))/=0) then 00635 cycle 00636 endif 00637 if (diagstrong(ii(weakend))) then !exchange 00638 j=ii(i) 00639 ii(i)=ii(weakend) 00640 ii(weakend)=j 00641 exit 00642 else 00643 cycle 00644 endif 00645 enddo 00646 endif 00647 endif 00648 !Add an aggregate 00649 aggr%inner%nagr=aggr%inner%nagr+1 00650 aggr%inner%num(ii(i))=aggr%inner%nagr 00651 if (nact<neighood) then 00652 nact=nact+1 00653 endif 00654 counter=counter+1 00655 if (counter>neighood) then 00656 counter=1 00657 endif 00658 actanum(counter)=aggr%inner%nagr 00659 aggrsize(counter)=1 00660 layersize(counter)=1 00661 layer(1,counter)=ii(i) 00662 stat(ii(i))=1 ! mark it 00663 if (nact>1) then !{ 00664 ! build another layer for still-to-grow aggregates: 00665 do j=1,nact-1 00666 k=counter-j 00667 if (k<1) then 00668 k=neighood 00669 endif 00670 newlayersize=0 00671 aa:do kk=1,layersize(k) 00672 layernode=layer(kk,k) 00673 rs=A%strong_rowstart(layernode) 00674 re=A%strong_rowstart(layernode+1)-1 00675 do jj=rs,re ! put the neighs into the list if needed 00676 cn=A%strong_colnrs(jj) 00677 if (stat(cn)==0) then !node not in any aggregate 00678 newlayersize=newlayersize+1 00679 newlayer(newlayersize)=cn 00680 aggr%inner%num(cn)=actanum(k) 00681 stat(cn)=1 ! mark it 00682 aggrsize(k)=aggrsize(k)+1 00683 if (aggrsize(k)>=maxaggrsize) then 00684 ! close this aggregate 00685 newlayersize=0 00686 exit aa 00687 endif 00688 endif 00689 enddo 00690 enddo aa 00691 ! keep only outermost layer always replacing it with new in the end 00692 layersize(k)=newlayersize 00693 layer(1:newlayersize,k)=newlayer(1:newlayersize) 00694 ! Now continue with a rounding step: 00695 if (rounding.and.aggrsize(k)<maxaggrsize) then 00696 newlayersize=0 00697 do kk=1,layersize(k) 00698 layernode=layer(kk,k) 00699 rs=A%strong_rowstart(layernode) 00700 re=A%strong_rowstart(layernode+1)-1 00701 do jj=rs,re ! put the neighs into the list if needed 00702 cn=A%strong_colnrs(jj) 00703 if (stat(cn)<=0) then !node not in any aggregate 00704 newlayersize=newlayersize+1 00705 newlayer(newlayersize)=cn 00706 stat(cn)=stat(cn)-1 ! count it 00707 endif 00708 enddo 00709 enddo 00710 ! add only rounding nodes 00711 bb:do kk=1,newlayersize 00712 if (stat(newlayer(kk))<-1) then 00713 layersize(k)=layersize(k)+1 00714 layer(layersize(k),k)=newlayer(kk) 00715 aggr%inner%num(newlayer(kk))=actanum(k) 00716 aggrsize(k)=aggrsize(k)+1 00717 if (aggrsize(k)>=maxaggrsize) then 00718 ! close this aggregate 00719 stat(newlayer(kk+1:newlayersize))=0 00720 layersize(k)=0 00721 exit bb 00722 endif 00723 else 00724 stat(newlayer(kk))=0 00725 endif 00726 enddo bb 00727 endif ! rounding 00728 enddo !j 00729 endif !} nact 00730 endif 00731 enddo 00732 if (sctls%plotting==1.or.sctls%plotting==3) then 00733 write(stream,*)'Rough Aggregates are:' 00734 call color_print_aggrs(n=n,aggrnum=aggr%inner%num,overwrite=overwrite) 00735 endif 00736 deallocate(A%strong_colnrs) 00737 deallocate(A%strong_rowstart) 00738 deallocate(diagstrong) 00739 deallocate(stat) 00740 deallocate(ii) 00741 00742 end subroutine SpMtx_roughly_aggregate 00743 00744 subroutine SpMtx_Build_lggl(A,A_ghost,M) 00745 implicit none 00746 type(SpMtx), intent(in out) :: A !Initial matrix 00747 type(SpMtx), intent(in out) :: A_ghost !matrix on ghost nodes 00748 type(Mesh) :: M !Mesh object 00749 00750 integer :: i,j,k,ntobsent,ninonol,ninner,indepoutol,nlf 00751 ! we are organising local freedoms as follows: 00752 !1,2,...,M%ntobsent,...,M%ninonol,...,M%ninner,...,M%indepoutol,...,M%nlf| 00753 !<-feedoms4send -> |<-rest inol->|<-independ.>|<-indep.onoutol>|<receivd>| 00754 !<- inner overlap -> |<-freedoms->|<- outer overlap ->| 00755 ! the goal: 00756 ! -1 -2 -3 2 1 | 00757 !-- the algorithm for marking the freedoms:--start with -3, then: -------+ 00758 ! M%ol_inner (-2) |############| M%ol_outer (2) | 00759 ! M%ax_sendidx | !M%recvidx| 00760 ! set to -1 | | set to 1| 00761 !------------------------------------------------------------------------+ 00762 M%ngf=max(A%nrows,A_ghost%nrows) 00763 allocate(M%gl_fmap(M%ngf)) 00764 M%gl_fmap=0 00765 do i=1,A%nnz 00766 M%gl_fmap(A%indi(i))=-3 00767 M%gl_fmap(A%indj(i))=-3 00768 enddo 00769 !do k=1,M%nnghbrs 00770 ! do i=1,M%ol_solve(k)%ninds 00771 ! M%gl_fmap(M%ol_solve(k)%inds(i))=2 00772 ! enddo 00773 !enddo 00774 do k=1,M%nnghbrs 00775 do i=1,M%ol_inner(k)%ninds 00776 M%gl_fmap(M%ol_inner(k)%inds(i))=-2 00777 enddo 00778 do i=1,M%ol_outer(k)%ninds 00779 M%gl_fmap(M%ol_outer(k)%inds(i))=2 00780 enddo 00781 enddo 00782 do k=1,M%nnghbrs 00783 do i=1,M%ax_sendidx(k)%ninds 00784 j=M%ax_sendidx(k)%inds(i) 00785 M%gl_fmap(j)=-1 00786 enddo 00787 do i=1,M%ax_recvidx(k)%ninds 00788 j=M%ax_recvidx(k)%inds(i) 00789 M%gl_fmap(j)=1 00790 enddo 00791 enddo 00792 !write(stream,*)'global freedom organisation:' 00793 !do i=1,M%ngf 00794 ! write(stream,*)i,':::',M%gl_fmap(i) 00795 !enddo 00796 00797 ! count how many corresponding freedoms there are: 00798 ntobsent=0 ! -1 00799 ninonol=0 ! -2 00800 ninner=0 ! -3 00801 indepoutol=0 ! 2 00802 nlf=0 ! 1 00803 do i=1,M%ngf 00804 if (M%gl_fmap(i)==-1) then 00805 ntobsent=ntobsent+1 00806 elseif (M%gl_fmap(i)==-2) then 00807 ninonol=ninonol+1 00808 elseif (M%gl_fmap(i)==-3) then 00809 ninner=ninner+1 00810 elseif (M%gl_fmap(i)==2) then 00811 indepoutol=indepoutol+1 00812 elseif (M%gl_fmap(i)==1) then 00813 nlf=nlf+1 00814 else 00815 M%gl_fmap(i)=0 00816 endif 00817 enddo 00818 !write(stream,*)'ntobsent=',ntobsent 00819 !write(stream,*)'ninonol=',ninonol 00820 !write(stream,*)'ninner=',ninner 00821 !write(stream,*)'indepoutol=',indepoutol 00822 !write(stream,*)'nlf=',nlf 00823 j=ntobsent+ninonol+ninner+indepoutol+nlf 00824 allocate(M%lg_fmap(ntobsent+ninonol+ninner+indepoutol+nlf)) 00825 ! now put also the rest of the freedom numbers in: 00826 ! [re]start the counters: 00827 nlf=ntobsent+ninonol+ninner+indepoutol 00828 indepoutol=ntobsent+ninonol+ninner 00829 ninner=ntobsent+ninonol 00830 ninonol=ntobsent 00831 ntobsent=0 00832 do i=1,M%ngf 00833 if (M%gl_fmap(i)==-1) then 00834 ntobsent=ntobsent+1 00835 M%gl_fmap(i)=ntobsent 00836 M%lg_fmap(ntobsent)=i 00837 elseif (M%gl_fmap(i)==-2) then 00838 ninonol=ninonol+1 00839 M%gl_fmap(i)=ninonol 00840 M%lg_fmap(ninonol)=i 00841 elseif (M%gl_fmap(i)==-3) then 00842 ninner=ninner+1 00843 M%gl_fmap(i)=ninner 00844 M%lg_fmap(ninner)=i 00845 elseif (M%gl_fmap(i)==2) then 00846 indepoutol=indepoutol+1 00847 M%gl_fmap(i)=indepoutol 00848 M%lg_fmap(indepoutol)=i 00849 elseif (M%gl_fmap(i)==1) then 00850 nlf=nlf+1 00851 M%gl_fmap(i)=nlf 00852 M%lg_fmap(nlf)=i 00853 endif 00854 enddo 00855 M%ntobsent=ntobsent 00856 M%ninonol=ninonol 00857 M%ninner=ninner 00858 M%indepoutol=indepoutol 00859 M%nlf=nlf 00860 if (nlf/=j) then 00861 call DOUG_abort('SpMtx_Build_lggl: nlf not right...',756) 00862 endif 00863 if (sctls%verbose>3) then 00864 write(stream,*)'inner freedoms:',M%lg_fmap(1:M%ninner) 00865 endif 00866 end subroutine SpMtx_Build_lggl 00867 00870 recursive subroutine SpMtx_addFront(sendmask, A, M, i, p, ol) 00871 implicit none 00872 integer(kind=1), dimension(:), pointer :: sendmask !< distribution mask 00873 type(SpMtx), intent(in) :: A !< original matrix 00874 type(Mesh), intent(in) :: M !< mesh corresponding to A 00875 integer, intent(in) :: i !< freedom 00876 integer, intent(in) :: p !< processor id (1..) 00877 integer, intent(in) :: ol !< overlap, non-negative 00878 !------------------------------------------- 00879 integer :: j, k 00880 00881 sendmask(p+i*numprocs)=max(ol+1, sendmask(p+i*numprocs)) 00882 if (ol>0) then 00883 do k=A%M_bound(i),A%M_bound(i+1)-1 00884 j=A%indj(k) 00885 if (p/=M%eptnmap(j) .and. sendmask(p+j*numprocs)<ol) then 00886 call SpMtx_addFront(sendmask, A, M, j, p, ol-1) 00887 end if 00888 end do 00889 end if 00890 end subroutine SpMtx_addFront 00891 00895 subroutine SpMtx_distributeWithOverlap(A, b, M, ol) 00896 implicit none 00897 type(SpMtx), intent(inout) :: A !< original matrix in case of master; in case of slave matrix data is ignored but structure (dimensions, symmetry, etc) should be initialized 00898 float(kind=rk), dimension(:), pointer :: b !< original RHS vector in case of master; in case of slave, b should be large enough to contain whole RHS (contents are ignored) 00899 type(Mesh), intent(in) :: M !< mesh corresponding to A 00900 integer, intent(in) :: ol !< overlap 00901 !------------------------------------ 00902 integer(kind=1), dimension(:), pointer :: sendmask 00903 integer, dimension(:), pointer :: nnz_p 00904 float(kind=rk), dimension(:), pointer :: val 00905 integer, dimension(:), pointer :: indi, indj 00906 integer :: i, j, k, p, ierr, status(MPI_STATUS_SIZE) 00907 00908 ! Number of nodes to distribute to each processor 00909 allocate(nnz_p(numprocs)) 00910 00911 if (ismaster()) then 00912 call SpMtx_arrange(A,arrange_type=D_SpMtx_ARRNG_ROWS,sort=.true.) 00913 00914 ! Create send vector for each node 00915 allocate(sendmask((A%nrows+1)*numprocs)) 00916 sendmask=0 00917 do k=1,A%nnz 00918 i=A%indi(k) 00919 call SpMtx_addFront(sendmask, A, M, i, M%eptnmap(i), max(1,ol*2)) 00920 end do 00921 00922 ! Calculate total number of nodes for each process 00923 nnz_p=0 00924 do p=1,numprocs 00925 do i=1,A%nrows 00926 if (sendmask(p+i*numprocs)>0) then 00927 do k=A%M_bound(i),A%M_bound(i+1)-1 00928 if (sendmask(p+A%indj(k)*numprocs)>0) then 00929 nnz_p(p)=nnz_p(p)+1 00930 end if 00931 end do 00932 end if 00933 end do 00934 end do 00935 00936 ! Debugging 00937 if (sctls%verbose>3) then 00938 write(stream,*) 'nnz_p:',nnz_p(1:numprocs) 00939 end if 00940 end if 00941 00942 ! Broadcast number of nodes for each processor 00943 call MPI_BCAST(nnz_p, numprocs, MPI_INTEGER, D_MASTER, MPI_COMM_WORLD, ierr) 00944 00945 ! Send/receive matrix nodes 00946 if (ismaster()) then 00947 do p=numprocs,1,-1 00948 allocate(val(nnz_p(p)), indi(nnz_p(p)), indj(nnz_p(p))) 00949 j=0 00950 do i=1,A%nrows 00951 if (sendmask(p+i*numprocs)>0) then 00952 do k=A%M_bound(i),A%M_bound(i+1)-1 00953 if (sendmask(p+A%indj(k)*numprocs)>0) then 00954 j=j+1 00955 val(j)=A%val(k) 00956 indi(j)=A%indi(k) 00957 indj(j)=A%indj(k) 00958 end if 00959 end do 00960 end if 00961 end do 00962 if (j/=nnz_p(p)) call DOUG_abort('[SpMtx_distributeWithOverlap] : j/=nnz_p(p)') 00963 if (p/=1) then 00964 call MPI_SEND(val, nnz_p(p), MPI_fkind, p-1, & 00965 D_TAG_ASSEMBLED_VALS, MPI_COMM_WORLD, status, ierr) 00966 call MPI_SEND(indi, nnz_p(p), MPI_INTEGER, p-1, & 00967 D_TAG_ASSEMBLED_IDXS_I, MPI_COMM_WORLD, status, ierr) 00968 call MPI_SEND(indj, nnz_p(p), MPI_INTEGER, p-1, & 00969 D_TAG_ASSEMBLED_IDXS_J, MPI_COMM_WORLD, status, ierr) 00970 deallocate(val, indi, indj) 00971 else 00972 deallocate(A%val, A%indi, A%indj) 00973 A%nnz=nnz_p(p) 00974 A%val=>val 00975 A%indi=>indi 00976 A%indj=>indj 00977 end if 00978 end do 00979 deallocate(sendmask) 00980 else 00981 p=myrank+1 00982 allocate(val(nnz_p(p)), indi(nnz_p(p)), indj(nnz_p(p))) 00983 call MPI_RECV(val, nnz_p(p), MPI_fkind, D_MASTER, & 00984 D_TAG_ASSEMBLED_VALS, MPI_COMM_WORLD, status, ierr) 00985 call MPI_RECV(indi, nnz_p(p), MPI_INTEGER, D_MASTER, & 00986 D_TAG_ASSEMBLED_IDXS_I, MPI_COMM_WORLD, status, ierr) 00987 call MPI_RECV(indj, nnz_p(p), MPI_INTEGER, D_MASTER, & 00988 D_TAG_ASSEMBLED_IDXS_J, MPI_COMM_WORLD, status, ierr) 00989 A%nnz=nnz_p(p) 00990 A%val=>val 00991 A%indi=>indi 00992 A%indj=>indj 00993 end if 00994 00995 deallocate(nnz_p) 00996 00997 ! rebuild A%M_bound 00998 if (associated(A%M_bound)) deallocate(A%M_bound) 00999 if (A%nnz>0) then 01000 allocate(A%M_bound(A%indi(A%nnz)+1)) 01001 j=1 01002 A%M_bound(j)=1 01003 do i=1,A%nnz 01004 do while (j<A%indi(i)) 01005 j=j+1 01006 A%M_bound(j)=i 01007 end do 01008 end do 01009 A%M_bound(j+1)=A%nnz+1 01010 if (j/=A%indi(A%nnz)) call DOUG_abort('[SpMtx_distributeWithOverlap] : j/=A%indi(A%nnz)') 01011 else 01012 allocate(A%M_bound(1)) 01013 A%M_bound(1)=1 01014 end if 01015 01016 A%arrange_type=D_SpMtx_ARRNG_ROWS 01017 01018 ! TODO: distribute only needed freedoms 01019 call MPI_BCAST(b, A%nrows, MPI_fkind, D_MASTER, MPI_COMM_WORLD, ierr) 01020 01021 end subroutine SpMtx_distributeWithOverlap 01022 01028 subroutine SpMtx_build_ghost(clr,ol,A,A_ghost,M,clrorder,clrstarts) 01029 implicit none 01030 !use SpMtx_class, only: indlist 01031 integer,intent(in) :: clr !the color # we are keeping 01032 integer,intent(in) :: ol !overlap size 01033 type(SpMtx), intent(in out) :: A !Initial matrix 01034 type(SpMtx), intent(in out) :: A_ghost !matrix on ghost nodes 01035 type(Mesh) :: M !Mesh object 01036 integer,dimension(:),pointer :: clrorder 01037 !order of matrix rows (columns) so that color i is found in rows (columns): 01038 ! clrorder(clrstarts(i):clrstarts(i+1)-1) 01039 integer,dimension(:),pointer :: clrstarts !(allocated earlier) 01040 !local: 01041 integer :: i,j,jj,k,clrnode,clrneigh,nfront,layer,lastlayer,neigh,node,nnz 01042 integer :: maxleadind,sendcnt,nfront1,sendnodecnt,recvcnt,neighnum,ol0nnz 01043 integer :: hl,offset,klayer 01044 integer,dimension(:),pointer :: neighmap,front 01045 integer,dimension(:),pointer :: onfront 01046 ! to each neighbour (Ax op) 01047 integer(kind=1),dimension(:),pointer :: sendnodes 01048 integer,dimension(:),pointer :: sendnodeidx ! to mark fred.s that 01049 ! will be communicated from my subdomain wherever (Ax op) 01050 integer,dimension(:),pointer :: frontstart,frontend 01051 integer :: a_ghostsz,a_gsz,ol0connfrom_outside,ol0connfrom_inside 01052 integer :: ol0cfo,ol0cfi,nol_on_neigh,nol_off_neigh 01053 integer,dimension(:),pointer :: itmp,jtmp,btmp 01054 float(kind=rk),dimension(:),pointer :: rtmp 01055 integer,dimension(:), pointer :: nnodesonclrol,ccount 01056 integer,dimension(2,2) :: bbe 01057 integer,dimension(:),pointer :: ol_outer 01058 01059 ! we assume this now: 01060 !! if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then 01061 !!elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then 01062 !! !TODO? 01063 !!else 01064 !! call DOUG_abort('SpMtx_keep_subd_wol: Matrix arrangment not done!',19) 01065 !!endif 01066 01067 allocate(neighmap(numprocs)) 01068 neighmap=0 01069 allocate(front(A%nrows)) ! for keeping track of the front 01070 allocate(onfront(A%nrows)) ! for keeping track on the front 01071 onfront=0 01072 lastlayer=max(ol,1) 01073 allocate(frontstart(0:2*lastlayer)) 01074 allocate(frontend(0:2*lastlayer)) 01075 allocate(nnodesonclrol(numprocs)) 01076 nnodesonclrol=0 01077 nfront=0 01078 ! mark clr nodes as the very first step. 01079 frontstart(0)=1 01080 do i=clrstarts(clr),clrstarts(clr+1)-1 01081 node=clrorder(i) 01082 if (onfront(node)==0) then 01083 nfront=nfront+1 01084 front(nfront)=node 01085 onfront(node)=-1 01086 endif 01087 enddo 01088 frontend(0)=nfront 01089 ! first, add 2*ol layers to the subdomain 01090 frontstart(1)=nfront+1 01091 do i=clrstarts(clr),clrstarts(clr+1)-1 01092 node=clrorder(i) 01093 do j=A%M_bound(node),A%M_bound(node+1)-1 01094 neigh=A%indj(j) 01095 clrneigh=M%eptnmap(neigh) 01096 if (clrneigh/=clr.and.onfront(neigh)==0) then 01097 onfront(neigh)=1 01098 nnodesonclrol(clrneigh)=nnodesonclrol(clrneigh)+1 01099 nfront=nfront+1 01100 front(nfront)=neigh 01101 endif 01102 enddo 01103 enddo 01104 frontend(1)=nfront 01105 if (ol>0) then 01106 !Now, let's go on with next layers...: 01107 do layer=2,2*lastlayer 01108 frontstart(layer)=nfront+1 01109 do i=frontstart(layer-1),frontend(layer-1) 01110 node=front(i) 01111 do j=A%M_bound(node),A%M_bound(node+1)-1 01112 neigh=A%indj(j) 01113 clrneigh=M%eptnmap(neigh) 01114 if (clrneigh/=clr.and.onfront(neigh)==0) then 01115 onfront(neigh)=layer 01116 nnodesonclrol(clrneigh)=nnodesonclrol(clrneigh)+1 01117 nfront=nfront+1 01118 front(nfront)=neigh 01119 endif 01120 enddo 01121 enddo 01122 frontend(layer)=nfront 01123 enddo 01124 else 01125 frontend(2)=frontend(1) 01126 endif 01127 ! now we are ready to decrease A the first time: 01128 nnz=0 01129 maxleadind=0 01130 do i=1,A%nnz 01131 if (onfront(A%indi(i))/=0.and. & 01132 onfront(A%indj(i))/=0) then 01133 nnz=nnz+1 01134 if (A%indi(i)>maxleadind) then 01135 maxleadind=A%indi(i) 01136 endif 01137 endif 01138 enddo 01139 allocate(itmp(nnz)) 01140 allocate(jtmp(nnz)) 01141 allocate(rtmp(nnz)) 01142 allocate(btmp(maxleadind+1)) 01143 btmp=0 01144 nnz=0 01145 do i=1,A%nnz 01146 if (onfront(A%indi(i))/=0.and. & 01147 onfront(A%indj(i))/=0) then 01148 nnz=nnz+1 01149 node=A%indi(i) 01150 itmp(nnz)=node 01151 btmp(node+1)=btmp(node+1)+1 01152 jtmp(nnz)=A%indj(i) 01153 rtmp(nnz)=A%val(i) 01154 endif 01155 enddo 01156 A%nnz=nnz 01157 deallocate(A%indi) 01158 allocate(A%indi(1:A%nnz)) 01159 A%indi(1:A%nnz)=itmp(1:nnz) 01160 deallocate(itmp) 01161 ! indj 01162 deallocate(A%indj) 01163 allocate(A%indj(1:A%nnz)) 01164 A%indj(1:A%nnz)=jtmp(1:nnz) 01165 deallocate(jtmp) 01166 !M_bound 01167 btmp(1)=1 01168 do i=2,maxleadind+1 01169 btmp(i)=btmp(i-1)+btmp(i) 01170 enddo 01171 deallocate(A%M_bound) 01172 allocate(A%M_bound(maxleadind+1)) 01173 A%M_bound=btmp 01174 deallocate(btmp) 01175 ! val 01176 deallocate(A%val) 01177 allocate(A%val(1:A%nnz)) 01178 A%val(1:A%nnz)=rtmp(1:nnz) 01179 deallocate(rtmp) 01180 ! take a look what are the neighbours: 01181 j=0 01182 do i=1,numprocs 01183 if (nnodesonclrol(i)>0) then 01184 j=j+1 01185 endif 01186 enddo 01187 M%nnghbrs=j 01188 allocate(M%nghbrs(M%nnghbrs+1)) 01189 j=0 01190 do i=1,numprocs 01191 if (nnodesonclrol(i)>0) then 01192 j=j+1 01193 M%nghbrs(j)=i-1 01194 neighmap(i)=j !shows now, where the subdomain is in M%nghbrs 01195 ! and is 0 if the subdomain is not neighbour 01196 ! (it is actually pid2indx in SpMtx_operation) 01197 elseif (i==clr) then 01198 M%nghbrs(M%nnghbrs+1)=i-1 01199 neighmap(i)=M%nnghbrs+1 01200 endif 01201 enddo 01202 ! now we are able to determine recvnodes already... i.e. nodes 01203 ! where onfront==lastlayer or actually the nodes at front( 01204 ! frontstart(lastlayer):frontend(lastlayer)) 01205 allocate(M%ax_recvidx(M%nnghbrs)) 01206 M%ax_recvidx%ninds = 0 01207 do i=frontstart(lastlayer),frontend(lastlayer) 01208 node=front(i) 01209 neighnum=neighmap(M%eptnmap(node)) 01210 M%ax_recvidx(neighnum)%ninds=M%ax_recvidx(neighnum)%ninds+1 01211 enddo 01212 do neighnum=1,M%nnghbrs 01213 allocate(M%ax_recvidx(neighnum)%inds(M%ax_recvidx(neighnum)%ninds)) 01214 enddo 01215 M%ax_recvidx%ninds = 0 01216 do i=frontstart(lastlayer),frontend(lastlayer) 01217 node=front(i) 01218 neighnum=neighmap(M%eptnmap(node)) 01219 M%ax_recvidx(neighnum)%ninds=M%ax_recvidx(neighnum)%ninds+1 01220 M%ax_recvidx(neighnum)%inds(M%ax_recvidx(neighnum)%ninds)=node 01221 enddo 01222 ! but now we need to rework clrorder aswell: 01223 allocate(ccount(M%nnghbrs+1)) 01224 ccount=0 01225 do i=frontstart(0),frontend(2*lastlayer) 01226 j=front(i) 01227 ccount(neighmap(M%eptnmap(j)))=ccount(neighmap(M%eptnmap(j)))+1 01228 enddo 01229 allocate(clrstarts(M%nnghbrs+2)) 01230 clrstarts(1)=1 01231 do i=1,M%nnghbrs+1 01232 clrstarts(i+1)=clrstarts(i)+ccount(i) 01233 end do 01234 allocate(clrorder(maxleadind)) 01235 ccount(1:M%nnghbrs+1)=clrstarts(1:M%nnghbrs+1)-1 01236 ! todo siin viga: also clr nodes need to be added! 01237 !rite(stream,*)'maxleadind=',maxleadind 01238 !rite(stream,*)'fff front=',front 01239 !rite(stream,*)'fff part(front)=',M%eptnmap(front) 01240 !rite(stream,*)'fff clrstarts=',clrstarts 01241 !rite(stream,*)'fff neighmap=',neighmap 01242 !rite(stream,*)'fff M%nghbrs=',M%nghbrs 01243 do i=frontstart(0),frontend(2*lastlayer) 01244 j=front(i) 01245 ccount(neighmap(M%eptnmap(j)))=ccount(neighmap(M%eptnmap(j)))+1 01246 clrorder(ccount(neighmap(M%eptnmap(j))))=j 01247 enddo 01248 if (sctls%verbose>3.and.A%nrows<200) then 01249 do i=1,M%nnghbrs+1 ! 01250 write(stream,*)'partition ',M%nghbrs(i),' is in:', & ! 01251 clrorder(clrstarts(i):clrstarts(i+1)-1) ! 01252 enddo ! 01253 endif 01254 deallocate(ccount) 01255 ! 01256 ! First, we must repair the onfront back to 0 from where -1n 01257 do i=frontstart(0),frontend(0) 01258 onfront(front(i))=0 01259 enddo 01260 ! Now we may start out from each neighbour individually concerning 01261 ! ol-outer, ol_inner, ol_solve and ax_sendidx 01262 ! mark clr nodes as the very first step. 01263 allocate(M%ax_sendidx(M%nnghbrs)) 01264 M%ax_sendidx%ninds = 0 01265 allocate(sendnodes(A%nrows)) ! indicator of nodes that will be sent to... 01266 sendnodes=0 ! whichever neighbour 01267 allocate(M%ol_inner(M%nnghbrs)) 01268 M%ol_inner%ninds = 0 01269 allocate(M%ol_outer(M%nnghbrs)) 01270 M%ol_outer%ninds = 0 01271 allocate(M%ol_solve(M%nnghbrs)) 01272 hl=2*lastlayer+1 ! the highest layer#+1 for offset def 01273 do k=1,M%nnghbrs 01274 nfront=0 01275 offset=k*hl 01276 clrnode=M%nghbrs(k)+1 01277 frontstart(0)=1 01278 !rite(stream,*)'starting with: onfront=',onfront 01279 do i=clrstarts(k),clrstarts(k+1)-1 01280 node=clrorder(i) 01281 !rite(stream,*)'node is:',node 01282 layer=modulo(onfront(node),hl) 01283 !rite(stream,*)'layer is:',layer 01284 if (layer<=lastlayer) then ! the node on ol_outer 01285 M%ol_outer(k)%ninds=M%ol_outer(k)%ninds+1 01286 endif 01287 onfront(node)=offset+layer 01288 nfront=nfront+1 01289 front(nfront)=node 01290 enddo 01291 frontend(0)=nfront 01292 do klayer=1,lastlayer 01293 frontstart(klayer)=nfront+1 01294 do i=frontstart(klayer-1),frontend(klayer-1) 01295 node=front(i) 01296 do j=A%M_bound(node),A%M_bound(node+1)-1 01297 neigh=A%indj(j) 01298 if (onfront(neigh)<offset) then ! the neigh not included yet 01299 layer=modulo(onfront(neigh),hl) 01300 if (layer==0) then ! the node inside the clr 01301 if (klayer==lastlayer) then ! the node for ax_sendidx 01302 M%ax_sendidx(k)%ninds=M%ax_sendidx(k)%ninds+1 01303 endif 01304 M%ol_inner(k)%ninds=M%ol_inner(k)%ninds+1 01305 elseif (layer<=lastlayer) then ! the node on ol_outer 01306 M%ol_outer(k)%ninds=M%ol_outer(k)%ninds+1 01307 endif 01308 onfront(neigh)=offset+layer 01309 nfront=nfront+1 01310 front(nfront)=neigh 01311 endif 01312 enddo 01313 enddo 01314 frontend(klayer)=nfront 01315 !rite(stream,*)'neigh:',k,' is:',M%nghbrs(k),' onfront=',onfront 01316 !rite(stream,*)'lastlayer:',lastlayer 01317 enddo 01318 allocate(M%ax_sendidx(k)%inds(M%ax_sendidx(k)%ninds)) 01319 allocate(M%ol_inner(k)%inds(M%ol_inner(k)%ninds)) 01320 allocate(ol_outer(M%ol_outer(k)%ninds)) 01321 M%ax_sendidx(k)%ninds = 0 01322 M%ol_inner(k)%ninds = 0 01323 M%ol_outer(k)%ninds = 0 01324 !do i=frontstart(1),frontend(lastlayer-1) 01325 !rite(stream,*)'the front in ',k,' is:',front(frontstart(0):frontend(lastlayer-1)) 01326 !rite(stream,*)'the front out ',k,' is:',front(frontstart(lastlayer):frontend(lastlayer)) 01327 do i=frontstart(0),frontend(lastlayer-1) 01328 node=front(i) 01329 layer=modulo(onfront(node),hl) 01330 if (layer==0) then ! the node inside the clr 01331 !rite(stream,*)'adding to inner: node=',node,' layer=',layer 01332 M%ol_inner(k)%ninds=M%ol_inner(k)%ninds+1 01333 M%ol_inner(k)%inds(M%ol_inner(k)%ninds)=node 01334 elseif (layer<=lastlayer) then ! the node on ol_outer 01335 !rite(stream,*)'adding to 0 outer: node=',node,' layer=',layer 01336 M%ol_outer(k)%ninds=M%ol_outer(k)%ninds+1 01337 ol_outer(M%ol_outer(k)%ninds)=node 01338 endif 01339 enddo 01340 do i=frontstart(lastlayer),frontend(lastlayer) 01341 node=front(i) 01342 layer=modulo(onfront(node),hl) 01343 if (layer==0) then ! the node inside the clr 01344 M%ax_sendidx(k)%ninds=M%ax_sendidx(k)%ninds+1 01345 M%ax_sendidx(k)%inds(M%ax_sendidx(k)%ninds)=node 01346 sendnodes(node)=1 ! the node added to the send pool 01347 !rite(stream,*)'adding to inner: node=',node,' layer=',layer 01348 M%ol_inner(k)%ninds=M%ol_inner(k)%ninds+1 01349 M%ol_inner(k)%inds(M%ol_inner(k)%ninds)=node 01350 elseif (layer<=lastlayer) then ! the node on ol_outer 01351 !rite(stream,*)'adding to outer: node=',node,' layer=',layer 01352 M%ol_outer(k)%ninds=M%ol_outer(k)%ninds+1 01353 ol_outer(M%ol_outer(k)%ninds)=node 01354 endif 01355 enddo 01356 call quicksort(M%ol_inner(k)%ninds,M%ol_inner(k)%inds) 01357 call quicksort(M%ax_sendidx(k)%ninds,M%ax_sendidx(k)%inds) 01358 call quicksort(M%ax_recvidx(k)%ninds,M%ax_recvidx(k)%inds) 01359 M%ol_solve(k)%ninds=M%ol_outer(k)%ninds+& 01360 M%ol_inner(k)%ninds 01361 allocate(M%ol_solve(k)%inds(M%ol_solve(k)%ninds)) 01362 j=M%ol_outer(k)%ninds 01363 M%ol_solve(k)%inds(1:j)=ol_outer(:) 01364 jj=j+M%ol_inner(k)%ninds 01365 M%ol_solve(k)%inds(j+1:jj)=M%ol_inner(k)%inds(:) 01366 call quicksort(M%ol_solve(k)%ninds,M%ol_solve(k)%inds) 01367 ! Now take out all foreign nodes from ol_outer: 01368 j=0 01369 do i=1,M%ol_outer(k)%ninds 01370 if (M%eptnmap(ol_outer(i))==clrnode) then 01371 j=j+1 01372 if (j<i) then 01373 ol_outer(j)=ol_outer(i) 01374 endif 01375 endif 01376 enddo 01377 M%ol_outer(k)%ninds=j 01378 allocate(M%ol_outer(k)%inds(M%ol_outer(k)%ninds)) 01379 M%ol_outer(k)%inds=ol_outer(1:j) 01380 deallocate(ol_outer) 01381 call quicksort(M%ol_outer(k)%ninds,M%ol_outer(k)%inds) 01382 if (sctls%verbose>3.and.A%nrows<300) then 01383 write(stream,*)k,':ol_solve:::',M%ol_solve(k)%inds 01384 write(stream,*)k,':ol_inner:::',M%ol_inner(k)%inds 01385 write(stream,*)k,':ax_sendidx:::',M%ax_sendidx(k)%inds 01386 write(stream,*)k,':ax_recvidx:::',M%ax_recvidx(k)%inds 01387 write(stream,*)k,':ol_outer:::',M%ol_outer(k)%inds 01388 endif 01389 enddo ! k 01390 ! note: actually, interf/inner may contain also interf/receive_nodes 01391 ! connections 01392 ! --------- ---------- 01393 ! | interf. | interf./ | 01394 ! | 11| inner 12| 01395 ! ---------+---------- 01396 ! |^ inner/ | | 01397 ! | interf.| inner | 01398 ! | 21| 22| 01399 ! ---------- ---------- 01400 ! 01401 ! onfront>lastlayer => M%ax_recvidx 01402 ! 01403 ! Count A arrays size 01404 A%mtx_bbe(1,1)=0 01405 A%mtx_bbe(1,2)=0 01406 A%mtx_bbe(2,1)=0 01407 A%mtx_bbe(2,2)=0 01408 a_ghostsz=0 01409 maxleadind=0 01410 ol0connfrom_outside=0 01411 ! start with the inner ones... 01412 !do i=clrstarts(clr),clrstarts(clr+1)-1 01413 do i=clrstarts(M%nnghbrs+1),clrstarts(M%nnghbrs+2)-1 01414 node=clrorder(i) 01415 if (sendnodes(node)==1) then ! node which will be sent, treat it as an 01416 ! interface node 01417 do j=A%M_bound(node),A%M_bound(node+1)-1 01418 neigh=A%indj(j) 01419 if (sendnodes(neigh)==1) then 01420 if (node>maxleadind) then 01421 maxleadind=node 01422 endif 01423 A%mtx_bbe(1,1)=A%mtx_bbe(1,1)+1 01424 elseif (M%eptnmap(neigh)==clr) then 01425 if (node>maxleadind) then 01426 maxleadind=node 01427 endif 01428 A%mtx_bbe(1,2)=A%mtx_bbe(1,2)+1 01429 elseif (ol==0) then ! this must be connection from outside 01430 if (node>maxleadind) then 01431 maxleadind=node 01432 endif 01433 ol0connfrom_outside=ol0connfrom_outside+1 01434 else 01435 if (node>maxleadind) then 01436 maxleadind=node 01437 endif 01438 A%mtx_bbe(1,2)=A%mtx_bbe(1,2)+1 01439 endif 01440 enddo 01441 else 01442 do j=A%M_bound(node),A%M_bound(node+1)-1 01443 neigh=A%indj(j) 01444 if (sendnodes(neigh)==1) then 01445 if (node>maxleadind) then 01446 maxleadind=node 01447 endif 01448 A%mtx_bbe(2,1)=A%mtx_bbe(2,1)+1 01449 else 01450 if (node>maxleadind) then 01451 maxleadind=node 01452 endif 01453 A%mtx_bbe(2,2)=A%mtx_bbe(2,2)+1 01454 endif 01455 enddo 01456 endif 01457 enddo 01458 ! now go outside: 01459 ol0connfrom_inside=0 01460 do k=1,M%nnghbrs 01461 !clrnode=M%nghbrs(k)+1 01462 do i=clrstarts(k),clrstarts(k+1)-1 01463 node=clrorder(i) 01464 layer=modulo(onfront(node),hl) 01465 if (layer==lastlayer) then ! gets value from comm. 01466 ! but we need to take it to the ghost matrix (if it is 01467 ! within the domain with overlap)! 01468 do j=A%M_bound(node),A%M_bound(node+1)-1 01469 neigh=A%indj(j) 01470 if (ol>0) then 01471 if (modulo(onfront(neigh),hl)<=lastlayer) then 01472 a_ghostsz=a_ghostsz+1 01473 endif 01474 elseif (sendnodes(neigh)==1) then !conn.from clr 01475 if (node>maxleadind) then 01476 maxleadind=node 01477 endif 01478 ol0connfrom_inside=ol0connfrom_inside+1 01479 endif 01480 enddo 01481 elseif (layer<lastlayer) then 01482 ! may still have conn. from internal node to be sent out 01483 do j=A%M_bound(node),A%M_bound(node+1)-1 01484 neigh=A%indj(j) 01485 if (sendnodes(neigh)==1) then 01486 if (node>maxleadind) then 01487 maxleadind=node 01488 endif 01489 A%mtx_bbe(2,1)=A%mtx_bbe(2,1)+1 01490 elseif (modulo(onfront(neigh),hl)<=lastlayer) then 01491 if (node>maxleadind) then 01492 maxleadind=node 01493 endif 01494 A%mtx_bbe(2,2)=A%mtx_bbe(2,2)+1 01495 endif! 01496 enddo ! 01497 endif 01498 enddo 01499 enddo 01500 A%mtx_bbs(1,1)=1!; A%mtx_bbe(1,1) remains as it is 01501 A%mtx_bbs(1,2)=A%mtx_bbe(1,1)+1 01502 A%mtx_bbe(1,2)=A%mtx_bbs(1,2)+A%mtx_bbe(1,2)-1 01503 A%mtx_bbs(2,1)=A%mtx_bbe(1,2)+1 01504 A%mtx_bbe(2,1)=A%mtx_bbs(2,1)+A%mtx_bbe(2,1)-1 01505 A%mtx_bbs(2,2)=A%mtx_bbe(2,1)+1 01506 A%mtx_bbe(2,2)=A%mtx_bbs(2,2)+A%mtx_bbe(2,2)-1 01507 if (ol==0) then ! We put the additional part of the matrix that does not 01508 ! participate in solves but still in the Ax-operation between 01509 ! A%mtx_bbe(2,2) and A%nnz 01510 nnz=A%mtx_bbe(2,2)+ol0connfrom_outside 01511 ol0nnz=nnz+ol0connfrom_inside 01512 !write(stream,*)'ol0connfrom_inside=',ol0connfrom_inside 01513 ol0cfo=A%mtx_bbe(2,2) 01514 ol0cfi=nnz 01515 else 01516 nnz=A%mtx_bbe(2,2) 01517 a_gsz=0 01518 endif 01519 bbe(1,1)=A%mtx_bbs(1,1)-1 01520 bbe(1,2)=A%mtx_bbs(1,2)-1 01521 bbe(2,1)=A%mtx_bbs(2,1)-1 01522 bbe(2,2)=A%mtx_bbs(2,2)-1 01523 A%nnz=nnz 01524 if (ol==0) then 01525 A%ol0nnz=ol0nnz 01526 else 01527 A%ol0nnz=nnz 01528 ol0nnz=nnz 01529 endif 01530 allocate(itmp(ol0nnz)) 01531 allocate(jtmp(ol0nnz)) 01532 allocate(rtmp(ol0nnz)) 01533 !write(stream,*)'maxleadind:',maxleadind 01534 allocate(btmp(maxleadind+1)) 01535 if (ol>0) then 01536 A_ghost=SpMtx_newInit(a_ghostsz) 01537 endif 01538 01539 ! start with the inner ones... 01540 !do i=clrstarts(clr),clrstarts(clr+1)-1 01541 do i=clrstarts(M%nnghbrs+1),clrstarts(M%nnghbrs+2)-1 01542 node=clrorder(i) 01543 if (sendnodes(node)==1) then ! node which will be sent, treat it as an 01544 ! interface node 01545 do j=A%M_bound(node),A%M_bound(node+1)-1 01546 neigh=A%indj(j) 01547 if (sendnodes(neigh)==1) then 01548 bbe(1,1)=bbe(1,1)+1 01549 btmp(node+1)=btmp(node+1)+1 01550 itmp(bbe(1,1))=node 01551 jtmp(bbe(1,1))=neigh 01552 rtmp(bbe(1,1))=A%val(j) 01553 elseif (M%eptnmap(neigh)==clr) then !!!! siin 01554 bbe(1,2)=bbe(1,2)+1 01555 btmp(node+1)=btmp(node+1)+1 01556 itmp(bbe(1,2))=node 01557 jtmp(bbe(1,2))=neigh 01558 rtmp(bbe(1,2))=A%val(j) 01559 elseif (ol==0) then ! this must be connection from outside 01560 ol0cfo=ol0cfo+1 01561 btmp(node+1)=btmp(node+1)+1 01562 itmp(ol0cfo)=node 01563 jtmp(ol0cfo)=neigh 01564 rtmp(ol0cfo)=A%val(j) 01565 else 01566 bbe(1,2)=bbe(1,2)+1 01567 btmp(node+1)=btmp(node+1)+1 01568 itmp(bbe(1,2))=node 01569 jtmp(bbe(1,2))=neigh 01570 rtmp(bbe(1,2))=A%val(j) 01571 endif 01572 enddo 01573 else 01574 do j=A%M_bound(node),A%M_bound(node+1)-1 01575 neigh=A%indj(j) 01576 if (sendnodes(neigh)==1) then 01577 bbe(2,1)=bbe(2,1)+1 01578 btmp(node+1)=btmp(node+1)+1 01579 itmp(bbe(2,1))=node 01580 jtmp(bbe(2,1))=neigh 01581 rtmp(bbe(2,1))=A%val(j) 01582 else 01583 bbe(2,2)=bbe(2,2)+1 01584 btmp(node+1)=btmp(node+1)+1 01585 itmp(bbe(2,2))=node 01586 jtmp(bbe(2,2))=neigh 01587 rtmp(bbe(2,2))=A%val(j) 01588 endif 01589 enddo 01590 endif 01591 enddo 01592 ! now go outside: 01593 do k=1,M%nnghbrs 01594 !clrnode=M%nghbrs(k)+1 01595 !do i=clrstarts(clrnode),clrstarts(clrnode+1)-1 01596 do i=clrstarts(k),clrstarts(k+1)-1 01597 node=clrorder(i) 01598 layer=modulo(onfront(node),hl) 01599 if (layer==lastlayer) then ! gets value from comm. 01600 ! but we need to take it to the ghost matrix (if it is 01601 ! within the domain with overlap)! 01602 do j=A%M_bound(node),A%M_bound(node+1)-1 01603 neigh=A%indj(j) 01604 if (ol>0) then 01605 if (modulo(onfront(neigh),hl)<=lastlayer) then 01606 a_gsz=a_gsz+1 01607 A_ghost%indi(a_gsz)=node 01608 A_ghost%indj(a_gsz)=neigh 01609 A_ghost%val(a_gsz)=A%val(j) 01610 endif 01611 elseif (sendnodes(neigh)==1) then !conn.from clr 01612 ol0cfi=ol0cfi+1 01613 btmp(node+1)=btmp(node+1)+1 01614 itmp(ol0cfi)=node 01615 jtmp(ol0cfi)=neigh 01616 rtmp(ol0cfi)=A%val(j) 01617 endif 01618 enddo 01619 elseif (layer<lastlayer) then 01620 ! may still have conn. from internal node to be sent out 01621 do j=A%M_bound(node),A%M_bound(node+1)-1 01622 neigh=A%indj(j) 01623 if (sendnodes(neigh)==1) then 01624 bbe(2,1)=bbe(2,1)+1 01625 btmp(node+1)=btmp(node+1)+1 01626 itmp(bbe(2,1))=node 01627 jtmp(bbe(2,1))=neigh 01628 rtmp(bbe(2,1))=A%val(j) 01629 elseif (modulo(onfront(neigh),hl)<=lastlayer) then 01630 bbe(2,2)=bbe(2,2)+1 01631 btmp(node+1)=btmp(node+1)+1 01632 itmp(bbe(2,2))=node 01633 jtmp(bbe(2,2))=neigh 01634 rtmp(bbe(2,2))=A%val(j) 01635 endif! 01636 enddo ! 01637 endif 01638 enddo 01639 enddo 01640 ! write(stream,*)'bbe(1,1),A%mtx_bbe(1,1):',bbe(1,1),A%mtx_bbe(1,1) 01641 ! write(stream,*)'A%bbs(1,2):',A%mtx_bbs(1,2) 01642 ! write(stream,*)'bbe(1,2),A%mtx_bbe(1,2):',bbe(1,2),A%mtx_bbe(1,2) 01643 ! write(stream,*)'A%bbs(2,1):',A%mtx_bbs(2,1) 01644 ! write(stream,*)'bbe(2,1),A%mtx_bbe(2,1):',bbe(2,1),A%mtx_bbe(2,1) 01645 ! write(stream,*)'A%bbs(2,2):',A%mtx_bbs(2,2) 01646 ! write(stream,*)'bbe(2,2),A%mtx_bbe(2,2):',bbe(2,2),A%mtx_bbe(2,2) 01647 ! write(stream,*)'a_gsz,a_ghostsz:',a_gsz,a_ghostsz 01648 if (bbe(1,1)/=A%mtx_bbe(1,1)) then 01649 write(stream,*)'bbe(1,1),A%mtx_bbe(1,1):',bbe(1,1),A%mtx_bbe(1,1) 01650 call DOUG_abort('SpMtx_build_ghost -- bbe(1,1) wrong!',67) 01651 endif 01652 if (bbe(1,2)/=A%mtx_bbe(1,2)) then 01653 write(stream,*)'bbe(1,2),A%mtx_bbe(1,2):',bbe(1,2),A%mtx_bbe(1,2) 01654 call DOUG_abort('SpMtx_build_ghost -- bbe(1,2) wrong!',67) 01655 endif 01656 if (bbe(2,1)/=A%mtx_bbe(2,1)) then 01657 write(stream,*)'bbe(2,1),A%mtx_bbe(2,1):',bbe(2,1),A%mtx_bbe(2,1) 01658 call DOUG_abort('SpMtx_build_ghost -- bbe(2,1) wrong!',67) 01659 endif 01660 if (bbe(2,2)/=A%mtx_bbe(2,2)) then 01661 write(stream,*)'bbe(2,2),A%mtx_bbe(2,2):',bbe(2,2),A%mtx_bbe(2,2) 01662 call DOUG_abort('SpMtx_build_ghost -- bbe(2,2) wrong!',67) 01663 endif 01664 if (ol==0) then 01665 if (ol0cfo/=ol0connfrom_outside+A%mtx_bbe(2,2)) then 01666 call DOUG_abort('SpMtx_build_ghost -- ol0cfo!',67) 01667 endif 01668 if (ol0cfi/=ol0connfrom_outside+A%mtx_bbe(2,2)+ol0connfrom_inside) then 01669 call DOUG_abort('SpMtx_build_ghost -- ol0cfi!',67) 01670 endif 01671 else 01672 if (a_gsz/=a_ghostsz) then 01673 write(stream,*)'a_gsz,a_ghostsz:',a_gsz,a_ghostsz 01674 call DOUG_abort('SpMtx_build_ghost -- a_gsz wrong!',67) 01675 endif 01676 endif 01677 01678 ! Resize actually the matrix: 01679 ! indi 01680 deallocate(A%indi) 01681 allocate(A%indi(1:A%ol0nnz)) 01682 A%indi(1:A%ol0nnz)=itmp(1:ol0nnz) 01683 deallocate(itmp) 01684 ! indj 01685 deallocate(A%indj) 01686 allocate(A%indj(1:A%ol0nnz)) 01687 A%indj(1:A%ol0nnz)=jtmp(1:ol0nnz) 01688 deallocate(jtmp) 01689 !M_bound 01690 btmp(1)=1 01691 do i=2,maxleadind+1 01692 btmp(i)=btmp(i-1)+btmp(i) 01693 enddo 01694 deallocate(A%M_bound) 01695 allocate(A%M_bound(maxleadind+1)) 01696 A%M_bound=btmp 01697 deallocate(btmp) 01698 ! val 01699 deallocate(A%val) 01700 allocate(A%val(1:A%ol0nnz)) 01701 A%val(1:A%ol0nnz)=rtmp(1:ol0nnz) 01702 deallocate(rtmp) 01703 01704 deallocate(frontend) 01705 deallocate(frontstart) 01706 deallocate(onfront) 01707 deallocate(front) 01708 deallocate(neighmap) 01709 end subroutine SpMtx_build_ghost 01710 01711 subroutine SpMtx_build_ghost_v01(clr,ol,A,A_ghost,M,clrorder,clrstarts) 01712 implicit none 01713 !use SpMtx_class, only: indlist 01714 integer,intent(in) :: clr !the color # we are keeping 01715 integer,intent(in) :: ol !overlap size 01716 Type(SpMtx), intent(in out) :: A !Initial matrix 01717 Type(SpMtx), intent(in out) :: A_ghost !matrix on ghost nodes 01718 type(Mesh) :: M !Mesh object 01719 integer,dimension(:),pointer :: clrorder 01720 !order of matrix rows (columns) so that color i is found in rows (columns): 01721 ! clrorder(clrstarts(i):clrstarts(i+1)-1) 01722 integer,dimension(:),pointer :: clrstarts !(allocated earlier) 01723 !local: 01724 integer :: i,j,jj,k,clrnode,clrneigh,nfront,layer,lastlayer,neigh,node,nnz 01725 integer :: maxleadind,sendcnt,nfront1,sendnodecnt,recvcnt 01726 integer,dimension(:),pointer :: neighmap,front 01727 integer,dimension(:),pointer :: onfront 01728 ! to each neighbour (Ax op) 01729 integer(kind=1),dimension(:),pointer :: sendnodes 01730 integer,dimension(:),pointer :: sendnodeidx ! to mark fred.s that 01731 ! will be communicated from my subdomain wherever (Ax op) 01732 integer,dimension(:),pointer :: frontstart,frontend 01733 integer,dimension(:,:),pointer :: neigfstart,neigfend 01734 integer :: a_ghostsz,a_gsz,ol0connfrom_outside 01735 integer :: ol0cfo,nol_on_neigh,nol_off_neigh 01736 integer,dimension(:),pointer :: itmp,jtmp,btmp 01737 float(kind=rk),dimension(:),pointer :: rtmp 01738 integer,dimension(:), pointer :: ndirectneigs 01739 integer,dimension(:,:), pointer :: directneigs 01740 integer,dimension(2,2) :: bbe 01741 type(indlist),dimension(:),pointer :: ol_outer_on !(off-neigh) 01742 01743 allocate(neighmap(numprocs)) 01744 neighmap=0 01745 allocate(front(A%nrows)) ! for keeping track of the front 01746 allocate(onfront(A%nrows)) ! for keeping track on the front 01747 onfront=0 01748 lastlayer=max(ol,1) 01749 allocate(frontstart(-lastlayer:0)) 01750 allocate(frontend(-lastlayer:0)) 01751 allocate(ndirectneigs(numprocs)) 01752 ndirectneigs=0 01753 nfront=0 01754 frontstart(0)=1 01755 ! first, count direct neighbours per each proc.: 01756 if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then 01757 do i=clrstarts(clr),clrstarts(clr+1)-1 01758 node=clrorder(i) 01759 do j=A%M_bound(node),A%M_bound(node+1)-1 01760 neigh=A%indj(j) 01761 clrneigh=M%eptnmap(neigh) 01762 if (clrneigh/=clr.and.onfront(neigh)/=-111) then 01763 onfront(neigh)=-111 01764 ndirectneigs(clrneigh)=ndirectneigs(clrneigh)+1 01765 if (onfront(node)==0) then 01766 nfront=nfront+1 01767 front(nfront)=node 01768 onfront(node)=1 01769 endif 01770 endif 01771 enddo 01772 enddo 01773 elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then 01774 !TODO? 01775 else 01776 call DOUG_abort('SpMtx_keep_subd_wol: Matrix arrangment not done!',19) 01777 endif 01778 frontend(0)=nfront 01779 !write(stream,*)'I have ',sum(neighmap),' neighbours!' 01780 ! counts # nodes on overlap with every other proc.: 01781 allocate(M%nfreesend_map(numprocs)) 01782 M%nfreesend_map=ndirectneigs 01783 j=0 01784 do i=1,numprocs 01785 if (ndirectneigs(i)>0) then 01786 j=j+1 01787 endif 01788 enddo 01789 if (j==0) then 01790 write(*,*) 'I am process', myrank,' but I have got no own freedoms!!!?' 01791 !write(*,*) ' clrstarts(clr),clrstarts(clr+1)-1:',clrstarts(clr),clrstarts(clr+1)-1 01792 call DOUG_abort('empty set of freedoms om process!',3433) 01793 endif 01794 M%nnghbrs=j 01795 allocate(M%nghbrs(M%nnghbrs)) 01796 allocate(neigfstart(1:lastlayer,M%nnghbrs)) 01797 allocate(neigfend(1:lastlayer,M%nnghbrs)) 01798 j=0 01799 do i=1,numprocs 01800 if (ndirectneigs(i)>0) then 01801 j=j+1 01802 M%nghbrs(j)=i-1 01803 neighmap(i)=j !shows now, where the subdomain is in M%nghbrs 01804 ! and is 0 if the subdomain is not neighbour 01805 ! (it is actually pid2indx in SpMtx_operation) 01806 endif 01807 enddo 01808 allocate(directneigs(maxval(ndirectneigs),M%nnghbrs)) 01809 allocate(M%ax_recvidx(M%nnghbrs)) 01810 M%ax_recvidx%ninds = 0 01811 allocate(M%ol_inner(M%nnghbrs)) 01812 allocate(M%ol_outer(M%nnghbrs)) 01813 allocate(M%ol_solve(M%nnghbrs)) 01814 allocate(ol_outer_on(M%nnghbrs)) 01815 M%ol_inner%ninds = 0 01816 M%ol_outer%ninds = 0 01817 M%ol_solve%ninds = 0 01818 ol_outer_on%ninds = 0 01819 M%nfreesend_map=0 01820 !write(stream,*)'My neighbours are: ',(M%nghbrs(i),i=1,M%nnghbrs) 01821 01822 !I overlap is 0, then only nodes on the fron are used as 01823 ! ghost values in Ax operation 01824 !if ol>0, then the subdomain expands in ol layers (for solves) 01825 01826 ! also communication structures are built here. 01827 ! 1. communication for Ax operation: 01828 ! if ol>0: 01829 ! -- only the outermost layer of ghost nodes get received 01830 ! ie. the innermost layer of my nodes to each neighbour 01831 ! need to be sent first, call them INTERFACE NODES (aswell, 01832 ! although they are actually inside ones) 01833 ! -- node of neighbour in M%ax_recvidx if: 01834 ! a) on outermost layer, OR 01835 ! b) the node has a connection to outside 01836 ! (might add some nodes in vain but this makes the 01837 ! algorithm deterministic from both sides...) 01838 ! ol==0: 01839 ! -- as bebore, the communication is done in the beginning! 01840 ! 01841 ! 2. communication after applying the subdomain solves in Prec. 01842 ! -- all the ghost values + inner overlap layers (ie. ghostnodes 01843 ! to neighbours) get communicated 01844 ! but: seems that now a few datastructures are in SpMtx_operation 01845 ! (fexchindx,pid2indx) and we cannot do it right away... 01846 ! ax_recvidx,ax_sendidx,ol_inner,ol_outer introduced instead 01847 01848 !Now, let's go on with next layers...: 01849 !if ol==0 we still need 1st layer for Ax operation... 01850 ndirectneigs=0 01851 if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then 01852 ! starting with layer=1 globally: 01853 ! look through the neighs of layer_0 nodes: 01854 do i=frontstart(0),frontend(0) 01855 node=front(i) 01856 do j=A%M_bound(node),A%M_bound(node+1)-1 01857 01858 neigh=A%indj(j) 01859 clrneigh=M%eptnmap(neigh) 01860 if (clrneigh/=clr) then 01861 if (onfront(neigh)==-111) then 01862 ndirectneigs(clrneigh)=ndirectneigs(clrneigh)+1 01863 directneigs(ndirectneigs(clrneigh),neighmap(clrneigh))=neigh 01864 nfront=nfront+1 01865 front(nfront)=neigh 01866 onfront(neigh)=1 01867 endif 01868 endif 01869 enddo 01870 enddo 01871 ! now go on with each neighbour individually: 01872 do k=1,M%nnghbrs 01873 clrnode=M%nghbrs(k)+1 01874 ! layer 1: 01875 if (k==1) then 01876 neigfstart(1,k)=1 01877 neigfend(1,k)=ndirectneigs(clrnode) 01878 nfront=neigfend(1,k) 01879 else 01880 neigfstart(1,k)=nfront+1 01881 neigfend(1,k)=nfront+ndirectneigs(clrnode) 01882 nfront=neigfend(1,k) 01883 endif 01884 if (ol<2) then 01885 front(neigfstart(1,k):neigfend(1,k))=& 01886 directneigs(1:ndirectneigs(clrnode),neighmap(clrnode)) 01887 onfront(front(neigfstart(1,k):neigfend(1,k)))=& 01888 numprocs+clrnode ! mark the nodes for recv 01889 recvcnt=ndirectneigs(clrnode) 01890 else 01891 front(neigfstart(1,k):neigfend(1,k))=& 01892 directneigs(1:ndirectneigs(clrnode),neighmap(clrnode)) 01893 onfront(front(neigfstart(1,k):neigfend(1,k)))=clrnode ! mark the nodes 01894 ! successive layers: 01895 recvcnt=0 01896 do layer=2,lastlayer-1 01897 neigfstart(layer,k)=nfront+1 01898 do i=neigfstart(layer-1,k),neigfend(layer-1,k) 01899 node=front(i) 01900 do j=A%M_bound(node),A%M_bound(node+1)-1 01901 neigh=A%indj(j) 01902 clrneigh=M%eptnmap(neigh) 01903 if (clrneigh==clrnode) then ! the same colour 01904 if (onfront(neigh)/=clrnode.and.& 01905 onfront(neigh)/=numprocs+clrnode) then 01906 onfront(neigh)=clrnode 01907 nfront=nfront+1 01908 front(nfront)=neigh 01909 endif 01910 elseif (clrneigh/=clr.and.onfront(node)/=numprocs+clrnode) then 01911 ! node to be included to recv part 01912 onfront(node)=numprocs+clrnode 01913 recvcnt=recvcnt+1 01914 endif 01915 enddo 01916 enddo 01917 neigfend(layer,k)=nfront 01918 enddo ! layer 01919 ! last layer: TODO siin: if ol==1 then we do not have neigfstart(0)!!! 01920 !!neigfstart(lastlayer,k)=nfront+1 01921 do i=neigfstart(lastlayer-1,k),neigfend(lastlayer-1,k) 01922 node=front(i) 01923 do j=A%M_bound(node),A%M_bound(node+1)-1 01924 neigh=A%indj(j) 01925 clrneigh=M%eptnmap(neigh) 01926 if (clrneigh==clrnode) then 01927 if (onfront(neigh)/=clrnode.and.& 01928 onfront(neigh)/=numprocs+clrnode) then 01929 ! neigh to be included to recv part 01930 onfront(neigh)=numprocs+clrnode 01931 nfront=nfront+1 01932 front(nfront)=neigh 01933 recvcnt=recvcnt+1 01934 endif 01935 elseif (clrneigh/=clr.and.onfront(node)/=numprocs+clrnode) then 01936 ! node to be included to recv part 01937 onfront(node)=numprocs+clrnode 01938 recvcnt=recvcnt+1 01939 endif 01940 enddo 01941 enddo 01942 endif ! ol 01943 neigfend(lastlayer,k)=nfront 01944 M%ax_recvidx(k)%ninds=recvcnt 01945 allocate(M%ax_recvidx(k)%inds(recvcnt)) 01946 recvcnt=0 01947 do i=neigfstart(1,k),neigfend(lastlayer,k) 01948 node=front(i) 01949 if (onfront(node)==numprocs+clrnode) then 01950 recvcnt=recvcnt+1 01951 !if (recvcnt>M%ax_recvidx(k)%ninds) then 01952 !write (stream,*)'so far aboutto receive:',M%ax_recvidx(k)%inds 01953 !write (stream,*)'wanna add:',node 01954 !write (stream,*)'neigfstart(1,k),neigfend(lastlayer,k):',neigfstart(1,k),neigfend(lastlayer,k) 01955 !write (stream,*)'front(neigfstart(1,k):neigfend(lastlayer,k)):',& 01956 ! front(neigfstart(1,k):neigfend(lastlayer,k)) 01957 !write (stream,*)'ndirectneigs(clrnode)=',ndirectneigs(clrnode) 01958 !endif 01959 M%ax_recvidx(k)%inds(recvcnt)=node 01960 endif 01961 enddo 01962 call quicksort(M%ax_recvidx(k)%ninds,M%ax_recvidx(k)%inds) 01963 if (sctls%verbose>3.and.A%nrows<200) then 01964 write(stream,*)myrank,'*** Ax:Recving from ',M%nghbrs(k),' nodes:',& 01965 M%ax_recvidx(k)%inds(1:M%ax_recvidx(k)%ninds) 01966 endif 01967 nol_on_neigh=neigfend(lastlayer,k)-neigfstart(1,k)+1 01968 M%ol_outer(k)%ninds=nol_on_neigh 01969 allocate(M%ol_outer(k)%inds(M%ol_outer(k)%ninds)) 01970 M%ol_outer(k)%inds(1:nol_on_neigh)=front(neigfstart(1,k):neigfend(lastlayer,k)) 01971 call quicksort(M%ol_outer(k)%ninds,M%ol_outer(k)%inds) 01972 if (sctls%verbose>3.and.A%nrows<200) then 01973 write(stream,*)myrank,'*** outer OL with: ',M%nghbrs(k),' is:',& 01974 M%ol_outer(k)%inds(1:M%ol_outer(k)%ninds) 01975 !write(stream,*)'ooo2 onfront=',onfront 01976 endif 01977 enddo 01978 if (ol>0) then ! expand the ol_outer: 01979 nfront1=nfront 01980 do k=1,M%nnghbrs 01981 !write(stream,*)'ooo onfront=',onfront 01982 !TODO: ol_outer to get expanded to the ol_outer of the neigh. as well 01983 ! for this we go through *from all nodes* of the neighbour with a 01984 ! depth-first algorithm upto ol expansion 01985 ! to mark the nodes that belong to some ol_outer of the current clr 01986 ! -- for storage of such nodes we (re)use still free part of front 01987 ! -- everywhere on clr ol_outer - onfront>0 01988 clrnode=M%nghbrs(k)+1 01989 do i=clrstarts(clrnode),clrstarts(clrnode+1)-1 01990 node=clrorder(i) 01991 do j=A%M_bound(node),A%M_bound(node+1)-1 01992 neigh=A%indj(j) 01993 clrneigh=M%eptnmap(neigh) 01994 if (clrneigh/=clr.and.clrneigh/=clrnode) then 01995 if (onfront(neigh)>0) then ! node from the clr outer overlap 01996 if (onfront(neigh)/=2*numprocs+clrnode) then ! this node 01997 ! has not been included yet... 01998 nfront=nfront+1 01999 front(nfront)=neigh 02000 onfront(neigh)=2*numprocs+clrnode 02001 !write(stream,*) 'adding node neigh=',neigh,front(nfront),nfront 02002 !write(stream,*)'ooo3 onfront=',onfront 02003 endif 02004 if (ol>1) then 02005 write(stream,*)'warning: ol>1 not completed yet.' 02006 ! go on with successive neighbours recursively keeping to 02007 ! the colour clrneigh 02008 endif 02009 endif 02010 endif 02011 enddo 02012 enddo 02013 !nol_on_neigh=neigfend(lastlayer,k)-neigfstart(1,k)+1 02014 nol_off_neigh=nfront-nfront1 02015 ol_outer_on(k)%ninds=nol_off_neigh 02016 allocate(ol_outer_on(k)%inds(ol_outer_on(k)%ninds)) 02017 ol_outer_on(k)%inds(1:ol_outer_on(k)%ninds)=front(nfront1+1:nfront) 02018 !call quicksort(ol_outer_on(k)%ninds,ol_outer_on(k)%inds) 02019 if (sctls%verbose>3.and.A%nrows<200) then 02020 write(stream,*)myrank,'*** outer OL off-neigh: ',M%nghbrs(k),' is:',& 02021 ol_outer_on(k)%inds(1:ol_outer_on(k)%ninds) 02022 !write(stream,*)'ooo2 onfront=',onfront 02023 endif 02024 nfront=nfront1 02025 enddo ! loop over neighbours 02026 !call MPI_BARRIER(MPI_COMM_WORLD,i) 02027 !call DOUG_abort('[SpMtx_arrangement] : testing ol_outer', -1) 02028 endif 02029 elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then 02030 !TODO 02031 endif 02032 02033 ! Now do the opposite part (in the inside) 02034 allocate(M%ax_sendidx(M%nnghbrs)) 02035 M%ax_sendidx%ninds = 0 02036 allocate(sendnodes(A%nrows)) ! indicator of nodes that will be sent to... 02037 allocate(sendnodeidx(A%nrows)) ! indeces of nodes that will be sent to 02038 sendnodeidx=0 02039 sendnodes=0 ! whichever neighbour 02040 sendnodecnt=0 02041 nfront1=nfront 02042 if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then 02043 do k=1,M%nnghbrs 02044 sendcnt=0 02045 clrnode=M%nghbrs(k)+1 02046 nfront=nfront1 02047 frontstart(-1)=nfront+1 02048 if (ol>1) then 02049 ! Layer 1 first: 02050 do i=1,ndirectneigs(M%nghbrs(k)+1) 02051 node=directneigs(i,k) 02052 do j=A%M_bound(node),A%M_bound(node+1)-1 02053 neigh=A%indj(j) 02054 clrneigh=M%eptnmap(neigh) 02055 if (clrneigh==clr) then ! my node 02056 if (onfront(neigh)/=-clrnode) then 02057 nfront=nfront+1 02058 front(nfront)=neigh 02059 onfront(neigh)=-clrnode 02060 endif 02061 endif 02062 enddo 02063 enddo 02064 frontend(-1)=nfront 02065 ! successive layers: 02066 do layer=2,lastlayer-1 02067 frontstart(-layer)=nfront+1 02068 do i=frontstart(-layer+1),frontend(-layer+1) 02069 node=front(i) 02070 do j=A%M_bound(node),A%M_bound(node+1)-1 02071 neigh=A%indj(j) 02072 clrneigh=M%eptnmap(neigh) 02073 if (clrneigh==clr) then ! my node 02074 if (onfront(neigh)/=-clrnode.and.& 02075 onfront(neigh)/=-numprocs-clrnode) then 02076 onfront(neigh)=-clrnode 02077 nfront=nfront+1 02078 front(nfront)=neigh 02079 endif 02080 elseif (clrneigh/=clrnode.and.onfront(node)/=-numprocs-clrnode) then 02081 ! node to be included to send part 02082 onfront(node)=-numprocs-clrnode 02083 sendcnt=sendcnt+1 02084 if (sendnodes(node)==0) then 02085 sendnodes(node)=1 02086 sendnodecnt=sendnodecnt+1 02087 sendnodeidx(sendnodecnt)=node 02088 endif 02089 endif 02090 enddo 02091 enddo 02092 frontend(-layer)=nfront 02093 enddo ! layer 02094 frontstart(-lastlayer)=nfront+1 02095 ! last layer: 02096 do i=frontstart(-lastlayer+1),frontend(-lastlayer+1) 02097 node=front(i) 02098 do j=A%M_bound(node),A%M_bound(node+1)-1 02099 neigh=A%indj(j) 02100 clrneigh=M%eptnmap(neigh) 02101 if (clrneigh==clr) then 02102 if (onfront(neigh)/=-clrnode.and.& 02103 onfront(neigh)/=-numprocs-clrnode) then 02104 ! neigh to be included to send part 02105 onfront(neigh)=-numprocs-clrnode 02106 nfront=nfront+1 02107 front(nfront)=neigh 02108 sendcnt=sendcnt+1 02109 if (sendnodes(neigh)==0) then 02110 sendnodes(neigh)=1 02111 sendnodecnt=sendnodecnt+1 02112 sendnodeidx(sendnodecnt)=neigh 02113 endif 02114 endif 02115 elseif (clrneigh/=clrnode.and.onfront(node)/=-numprocs-clrnode) then 02116 ! node to be included to send part 02117 onfront(node)=-numprocs-clrnode 02118 sendcnt=sendcnt+1 02119 if (sendnodes(node)==0) then 02120 sendnodes(node)=1 02121 sendnodecnt=sendnodecnt+1 02122 sendnodeidx(sendnodecnt)=node 02123 endif 02124 endif 02125 enddo 02126 enddo 02127 frontend(-lastlayer)=nfront 02128 else !(ol<2) 02129 ! Layer 1 only: 02130 do i=1,ndirectneigs(M%nghbrs(k)+1) 02131 node=directneigs(i,k) 02132 do j=A%M_bound(node),A%M_bound(node+1)-1 02133 neigh=A%indj(j) 02134 clrneigh=M%eptnmap(neigh) 02135 if (clrneigh==clr) then 02136 if (onfront(neigh)/=-numprocs-clrnode) then 02137 ! neigh to be included to send part 02138 onfront(neigh)=-numprocs-clrnode 02139 nfront=nfront+1 02140 front(nfront)=neigh 02141 sendcnt=sendcnt+1 02142 if (sendnodes(neigh)==0) then 02143 sendnodes(neigh)=1 02144 sendnodecnt=sendnodecnt+1 02145 sendnodeidx(sendnodecnt)=neigh 02146 endif 02147 endif 02148 endif 02149 enddo 02150 enddo 02151 frontend(-1)=nfront 02152 endif 02153 M%ax_sendidx(k)%ninds=sendcnt 02154 allocate(M%ax_sendidx(k)%inds(sendcnt)) 02155 sendcnt=0 02156 if (ol>0) then 02157 M%ol_inner(k)%ninds=frontend(-lastlayer)-frontstart(-1)+1 02158 allocate(M%ol_inner(k)%inds(M%ol_inner(k)%ninds)) 02159 M%ol_inner(k)%inds=front(frontstart(-1):frontend(-lastlayer)) 02160 call quicksort(M%ol_inner(k)%ninds,M%ol_inner(k)%inds) 02161 if (sctls%verbose>3.and.A%nrows<200) then 02162 write(stream,*)myrank,'*** inner OL with: ',M%nghbrs(k),' is:',& 02163 M%ol_inner(k)%inds(1:M%ol_inner(k)%ninds) 02164 endif 02165 M%ol_solve(k)%ninds=M%ol_outer(k)%ninds+& 02166 ol_outer_on(k)%ninds+& 02167 M%ol_inner(k)%ninds 02168 allocate(M%ol_solve(k)%inds(M%ol_solve(k)%ninds)) 02169 j=M%ol_outer(k)%ninds 02170 M%ol_solve(k)%inds(1:j)=M%ol_outer(k)%inds(:) 02171 jj=j+ol_outer_on(k)%ninds 02172 M%ol_solve(k)%inds(j+1:jj)=ol_outer_on(k)%inds(:) 02173 deallocate(ol_outer_on(k)%inds) 02174 j=jj+M%ol_inner(k)%ninds 02175 M%ol_solve(k)%inds(jj+1:j)=M%ol_inner(k)%inds(:) 02176 call quicksort(M%ol_solve(k)%ninds,M%ol_solve(k)%inds) 02177 if (sctls%verbose>3.and.A%nrows<200) then 02178 write(stream,*)myrank,'*** solve OL with: ',M%nghbrs(k),' is:',& 02179 M%ol_solve(k)%inds(1:M%ol_solve(k)%ninds) 02180 endif 02181 endif 02182 do i=frontstart(-1),frontend(-lastlayer) 02183 node=front(i) 02184 if (onfront(node)==-numprocs-clrnode) then 02185 sendcnt=sendcnt+1 02186 M%ax_sendidx(k)%inds(sendcnt)=node 02187 endif 02188 enddo 02189 call quicksort(M%ax_sendidx(k)%ninds,M%ax_sendidx(k)%inds) 02190 if (sctls%verbose>3.and.A%nrows<200) then 02191 write(stream,*)myrank,'*** Ax:Sending to ',M%nghbrs(k),' nodes:',& 02192 M%ax_sendidx(k)%inds(1:M%ax_sendidx(k)%ninds) 02193 endif 02194 enddo ! loop over neighbours 02195 elseif (A%arrange_type==D_SpMtx_ARRNG_COLS) then 02196 !TODO 02197 endif 02198 deallocate(ol_outer_on) 02199 ! note: actually, interf/inner may contain also interf/receive_nodes 02200 ! connections 02201 ! --------- ---------- 02202 ! | interf. | interf./ | 02203 ! | 11| inner 12| 02204 ! ---------+---------- 02205 ! |^ inner/ | | 02206 ! | interf.| inner | 02207 ! | 21| 22| 02208 ! ---------- ---------- 02209 ! 02210 ! onfront>lastlayer => M%ax_recvidx 02211 ! 02212 ! Count A arrays size 02213 A%mtx_bbe(1,1)=0 02214 A%mtx_bbe(1,2)=0 02215 A%mtx_bbe(2,1)=0 02216 A%mtx_bbe(2,2)=0 02217 a_ghostsz=0 02218 maxleadind=0 02219 ol0connfrom_outside=0 02220 if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then 02221 ! start with the inner ones... 02222 do i=clrstarts(clr),clrstarts(clr+1)-1 02223 node=clrorder(i) 02224 if (sendnodes(node)==1) then ! node which will be sent, treat it as an 02225 ! interface node 02226 do j=A%M_bound(node),A%M_bound(node+1)-1 02227 neigh=A%indj(j) 02228 if (sendnodes(neigh)==1) then 02229 if (node>maxleadind) then 02230 maxleadind=node 02231 endif 02232 A%mtx_bbe(1,1)=A%mtx_bbe(1,1)+1 02233 elseif (M%eptnmap(neigh)==clr) then 02234 if (node>maxleadind) then 02235 maxleadind=node 02236 endif 02237 A%mtx_bbe(1,2)=A%mtx_bbe(1,2)+1 02238 elseif (ol==0) then ! this must be connection from outside 02239 if (node>maxleadind) then 02240 maxleadind=node 02241 endif 02242 ol0connfrom_outside=ol0connfrom_outside+1 02243 else 02244 if (node>maxleadind) then 02245 maxleadind=node 02246 endif 02247 A%mtx_bbe(1,2)=A%mtx_bbe(1,2)+1 02248 endif 02249 enddo 02250 else 02251 do j=A%M_bound(node),A%M_bound(node+1)-1 02252 neigh=A%indj(j) 02253 if (sendnodes(neigh)==1) then 02254 if (node>maxleadind) then 02255 maxleadind=node 02256 endif 02257 A%mtx_bbe(2,1)=A%mtx_bbe(2,1)+1 02258 !!elseif (M%eptnmap(neigh)==clr) then !!!! siin 02259 else 02260 if (node>maxleadind) then 02261 maxleadind=node 02262 endif 02263 A%mtx_bbe(2,2)=A%mtx_bbe(2,2)+1 02264 endif 02265 enddo 02266 endif 02267 enddo 02268 ! now go outside: 02269 do i=neigfstart(1,1),neigfend(lastlayer,M%nnghbrs) 02270 node=front(i) 02271 if (onfront(node)>numprocs) then ! gets value from comm. 02272 ! but we need to take it to the ghost matrix (if it is 02273 ! within the domain with overlap)! 02274 do j=A%M_bound(node),A%M_bound(node+1)-1 02275 neigh=A%indj(j) 02276 !if (ol>0.or.M%eptnmap(neigh)==clr) then 02277 if (ol>0) then 02278 !!if (onfront(neigh)>numprocs) then 02279 if (onfront(neigh)/=0) then 02280 a_ghostsz=a_ghostsz+1 02281 endif 02282 endif 02283 enddo 02284 else ! ordinary node 02285 ! may still have conn. from internal node to be sent out 02286 do j=A%M_bound(node),A%M_bound(node+1)-1 02287 neigh=A%indj(j) 02288 if (sendnodes(neigh)==1) then 02289 if (node>maxleadind) then 02290 maxleadind=node 02291 endif 02292 A%mtx_bbe(2,1)=A%mtx_bbe(2,1)+1 02293 elseif (onfront(neigh)/=0) then 02294 if (node>maxleadind) then 02295 maxleadind=node 02296 endif 02297 A%mtx_bbe(2,2)=A%mtx_bbe(2,2)+1 02298 endif! 02299 enddo ! 02300 endif 02301 enddo 02302 endif 02303 A%mtx_bbs(1,1)=1!; A%mtx_bbe(1,1) remains as it is 02304 A%mtx_bbs(1,2)=A%mtx_bbe(1,1)+1 02305 A%mtx_bbe(1,2)=A%mtx_bbs(1,2)+A%mtx_bbe(1,2)-1 02306 A%mtx_bbs(2,1)=A%mtx_bbe(1,2)+1 02307 A%mtx_bbe(2,1)=A%mtx_bbs(2,1)+A%mtx_bbe(2,1)-1 02308 A%mtx_bbs(2,2)=A%mtx_bbe(2,1)+1 02309 A%mtx_bbe(2,2)=A%mtx_bbs(2,2)+A%mtx_bbe(2,2)-1 02310 if (ol==0) then ! We put the additional part of the matrix that does not 02311 ! participate in solves but still in the Ax-operation between 02312 ! A%mtx_bbe(2,2) and A%nnz 02313 nnz=A%mtx_bbe(2,2)+ol0connfrom_outside 02314 ol0cfo=A%mtx_bbe(2,2) 02315 else 02316 nnz=A%mtx_bbe(2,2) 02317 a_gsz=0 02318 endif 02319 bbe(1,1)=A%mtx_bbs(1,1)-1 02320 bbe(1,2)=A%mtx_bbs(1,2)-1 02321 bbe(2,1)=A%mtx_bbs(2,1)-1 02322 bbe(2,2)=A%mtx_bbs(2,2)-1 02323 A%nnz=nnz 02324 allocate(itmp(nnz)) 02325 allocate(jtmp(nnz)) 02326 allocate(rtmp(nnz)) 02327 !write(stream,*)'maxleadind:',maxleadind 02328 allocate(btmp(maxleadind+1)) 02329 if (ol>0) then 02330 A_ghost=SpMtx_newInit(a_ghostsz) 02331 endif 02332 02333 if (A%arrange_type==D_SpMtx_ARRNG_ROWS) then 02334 ! start with the inner ones... 02335 do i=clrstarts(clr),clrstarts(clr+1)-1 02336 node=clrorder(i) 02337 if (sendnodes(node)==1) then ! node which will be sent, treat it as an 02338 ! interface node 02339 do j=A%M_bound(node),A%M_bound(node+1)-1 02340 neigh=A%indj(j) 02341 if (sendnodes(neigh)==1) then 02342 bbe(1,1)=bbe(1,1)+1 02343 btmp(node+1)=btmp(node+1)+1 02344 itmp(bbe(1,1))=node 02345 jtmp(bbe(1,1))=neigh 02346 rtmp(bbe(1,1))=A%val(j) 02347 elseif (M%eptnmap(neigh)==clr) then !!!! siin 02348 bbe(1,2)=bbe(1,2)+1 02349 btmp(node+1)=btmp(node+1)+1 02350 itmp(bbe(1,2))=node 02351 jtmp(bbe(1,2))=neigh 02352 rtmp(bbe(1,2))=A%val(j) 02353 elseif (ol==0) then ! this must be connection from outside 02354 ol0cfo=ol0cfo+1 02355 btmp(node+1)=btmp(node+1)+1 02356 itmp(ol0cfo)=node 02357 jtmp(ol0cfo)=neigh 02358 rtmp(ol0cfo)=A%val(j) 02359 else 02360 bbe(1,2)=bbe(1,2)+1 02361 btmp(node+1)=btmp(node+1)+1 02362 itmp(bbe(1,2))=node 02363 jtmp(bbe(1,2))=neigh 02364 rtmp(bbe(1,2))=A%val(j) 02365 endif 02366 enddo 02367 else 02368 do j=A%M_bound(node),A%M_bound(node+1)-1 02369 neigh=A%indj(j) 02370 if (sendnodes(neigh)==1) then 02371 bbe(2,1)=bbe(2,1)+1 02372 btmp(node+1)=btmp(node+1)+1 02373 itmp(bbe(2,1))=node 02374 jtmp(bbe(2,1))=neigh 02375 rtmp(bbe(2,1))=A%val(j) 02376 elseif (M%eptnmap(neigh)==clr) then !!!! siin 02377 bbe(2,2)=bbe(2,2)+1 02378 btmp(node+1)=btmp(node+1)+1 02379 itmp(bbe(2,2))=node 02380 jtmp(bbe(2,2))=neigh 02381 rtmp(bbe(2,2))=A%val(j) 02382 endif 02383 enddo 02384 endif 02385 enddo 02386 ! now go outside: 02387 do i=neigfstart(1,1),neigfend(lastlayer,M%nnghbrs) 02388 node=front(i) 02389 if (onfront(node)>numprocs) then ! gets value from comm. 02390 ! but we need to take it to the ghost matrix (if it is 02391 ! within the domain with overlap)! 02392 do j=A%M_bound(node),A%M_bound(node+1)-1 02393 neigh=A%indj(j) 02394 !!!!if (ol>0.or.M%eptnmap(neigh)==clr) then 02395 if (ol>0) then 02396 !!if (onfront(neigh)>numprocs) then 02397 if (onfront(neigh)/=0) then 02398 a_gsz=a_gsz+1 02399 A_ghost%indi(a_gsz)=node 02400 A_ghost%indj(a_gsz)=neigh 02401 A_ghost%val(a_gsz)=A%val(j) 02402 endif 02403 !elseif (M%eptnmap(neigh)==clr) then 02404 ! if (onfront(neigh)/=0) then 02405 ! a_gsz=a_gsz+1 02406 ! btmp(node+1)=btmp(node+1)+1 02407 ! itmp(a_gsz)=node 02408 ! jtmp(a_gsz)=neigh 02409 ! rtmp(a_gsz)=A%val(j) 02410 ! endif 02411 endif 02412 enddo 02413 else ! ordinary node 02414 ! may still have conn. from internal node to be sent out 02415 do j=A%M_bound(node),A%M_bound(node+1)-1 02416 neigh=A%indj(j) 02417 if (sendnodes(neigh)==1) then 02418 bbe(2,1)=bbe(2,1)+1 02419 btmp(node+1)=btmp(node+1)+1 02420 itmp(bbe(2,1))=node 02421 jtmp(bbe(2,1))=neigh 02422 rtmp(bbe(2,1))=A%val(j) 02423 elseif (onfront(neigh)/=0) then 02424 bbe(2,2)=bbe(2,2)+1 02425 btmp(node+1)=btmp(node+1)+1 02426 itmp(bbe(2,2))=node 02427 jtmp(bbe(2,2))=neigh 02428 rtmp(bbe(2,2))=A%val(j) 02429 endif 02430 enddo 02431 endif 02432 enddo 02433 endif 02434 ! write(stream,*)'bbe(1,1),A%mtx_bbe(1,1):',bbe(1,1),A%mtx_bbe(1,1) 02435 ! write(stream,*)'A%bbs(1,2):',A%mtx_bbs(1,2) 02436 ! write(stream,*)'bbe(1,2),A%mtx_bbe(1,2):',bbe(1,2),A%mtx_bbe(1,2) 02437 ! write(stream,*)'A%bbs(2,1):',A%mtx_bbs(2,1) 02438 ! write(stream,*)'bbe(2,1),A%mtx_bbe(2,1):',bbe(2,1),A%mtx_bbe(2,1) 02439 ! write(stream,*)'A%bbs(2,2):',A%mtx_bbs(2,2) 02440 ! write(stream,*)'bbe(2,2),A%mtx_bbe(2,2):',bbe(2,2),A%mtx_bbe(2,2) 02441 ! write(stream,*)'a_gsz,a_ghostsz:',a_gsz,a_ghostsz 02442 if (bbe(1,1)/=A%mtx_bbe(1,1)) then 02443 write(stream,*)'bbe(1,1),A%mtx_bbe(1,1):',bbe(1,1),A%mtx_bbe(1,1) 02444 call DOUG_abort('SpMtx_build_ghost -- bbe(1,1) wrong!',67) 02445 endif 02446 if (bbe(1,2)/=A%mtx_bbe(1,2)) then 02447 write(stream,*)'bbe(1,2),A%mtx_bbe(1,2):',bbe(1,2),A%mtx_bbe(1,2) 02448 call DOUG_abort('SpMtx_build_ghost -- bbe(1,2) wrong!',67) 02449 endif 02450 if (bbe(2,1)/=A%mtx_bbe(2,1)) then 02451 write(stream,*)'bbe(2,1),A%mtx_bbe(2,1):',bbe(2,1),A%mtx_bbe(2,1) 02452 call DOUG_abort('SpMtx_build_ghost -- bbe(2,1) wrong!',67) 02453 endif 02454 if (bbe(2,2)/=A%mtx_bbe(2,2)) then 02455 write(stream,*)'bbe(2,2),A%mtx_bbe(2,2):',bbe(2,2),A%mtx_bbe(2,2) 02456 call DOUG_abort('SpMtx_build_ghost -- bbe(2,2) wrong!',67) 02457 endif 02458 if (ol==0) then 02459 if (ol0cfo/=ol0connfrom_outside+A%mtx_bbe(2,2)) then 02460 write(stream,*)'ol0cfo,ol0connfrom_outside+A%mtx_bbe(2,2):',a_gsz,a_ghostsz+A%mtx_bbe(2,2) 02461 call DOUG_abort('SpMtx_build_ghost -- ol0cfo!',67) 02462 endif 02463 else 02464 if (a_gsz/=a_ghostsz) then 02465 write(stream,*)'a_gsz,a_ghostsz:',a_gsz,a_ghostsz 02466 call DOUG_abort('SpMtx_build_ghost -- a_gsz wrong!',67) 02467 endif 02468 endif 02469 02470 ! Resize actually the matrix: 02471 ! indi 02472 deallocate(A%indi) 02473 allocate(A%indi(1:A%nnz)) 02474 A%indi(1:A%nnz)=itmp(1:nnz) 02475 deallocate(itmp) 02476 ! indj 02477 deallocate(A%indj) 02478 allocate(A%indj(1:A%nnz)) 02479 A%indj(1:A%nnz)=jtmp(1:nnz) 02480 deallocate(jtmp) 02481 !M_bound 02482 btmp(1)=1 02483 do i=2,maxleadind+1 02484 btmp(i)=btmp(i-1)+btmp(i) 02485 enddo 02486 deallocate(A%M_bound) 02487 allocate(A%M_bound(maxleadind+1)) 02488 A%M_bound=btmp 02489 deallocate(btmp) 02490 ! val 02491 deallocate(A%val) 02492 allocate(A%val(1:A%nnz)) 02493 A%val(1:A%nnz)=rtmp(1:nnz) 02494 deallocate(rtmp) 02495 02496 deallocate(neigfend) 02497 deallocate(neigfstart) 02498 deallocate(directneigs) 02499 deallocate(ndirectneigs) 02500 deallocate(frontend) 02501 deallocate(frontstart) 02502 deallocate(onfront) 02503 deallocate(front) 02504 deallocate(neighmap) 02505 end subroutine SpMtx_build_ghost_v01 02506 02507 subroutine SpMtx_buildAdjncy(A,nedges,xadj,adjncy) 02508 use globals, only : stream, D_MSGLVL 02509 implicit none 02510 02511 type(SpMtx),intent(inout) :: A 02512 integer, intent(out) :: nedges 02513 integer, dimension(:), pointer :: xadj 02514 integer, dimension(:), pointer :: adjncy 02515 02516 integer :: i,k,s,s1,sadjncy 02517 integer, dimension(:), pointer :: counter 02518 02519 ! allocation for the adjacency data 02520 allocate(xadj(A%nrows+1)) 02521 xadj = 0 02522 ! count the lengths 02523 ! (NB! We are expecting matrix symmetric structure!!!) 02524 do k=1,A%nnz 02525 if (A%indi(k)<A%indj(k)) then 02526 xadj(A%indi(k))=xadj(A%indi(k))+1 02527 xadj(A%indj(k))=xadj(A%indj(k))+1 02528 endif 02529 enddo 02530 allocate(counter(A%nrows)) 02531 counter = 0 02532 s = xadj(1) 02533 xadj(1) = 1 02534 counter(1) = 1 02535 do i = 2,A%nrows 02536 s1 = xadj(i) 02537 xadj(i) = xadj(i-1)+s 02538 counter(i) = xadj(i) 02539 s = s1 02540 enddo 02541 xadj(A%nrows+1) = xadj(A%nrows) + s 02542 sadjncy = xadj(A%nrows+1) - 1 02543 allocate(adjncy(sadjncy)) 02544 ! pass 2 of the data 02545 do k=1,A%nnz 02546 if (A%indi(k)<A%indj(k)) then 02547 adjncy(counter(A%indi(k)))=A%indj(k) 02548 counter(A%indi(k))=counter(A%indi(k))+1 02549 adjncy(counter(A%indj(k)))=A%indi(k) 02550 counter(A%indj(k))=counter(A%indj(k))+1 02551 endif 02552 enddo 02553 nedges = sadjncy/2 02554 deallocate(counter) 02555 end subroutine SpMtx_buildAdjncy 02556 02558 subroutine SpMtx_buildAggrAdjncy(A,aggr,maxaggrsize,nedges,xadj,adjncy) 02559 use globals, only : stream, D_MSGLVL 02560 implicit none 02561 02562 type(SpMtx),intent(inout) :: A 02563 type(AggrInfo),intent(in) :: aggr !< aggregates 02564 integer, intent(in) :: maxaggrsize 02565 integer, intent(out) :: nedges 02566 integer, dimension(:), pointer :: xadj 02567 integer, dimension(:), pointer :: adjncy 02568 02569 integer :: i,k,s,s1,sadjncy,c1,c2 02570 integer, dimension(:), pointer :: counter,nneig 02571 integer, dimension(:,:), allocatable :: neigs 02572 integer :: nnodes 02573 02574 !At first, find, which are the neighbouring aggregates 02575 allocate(neigs(maxaggrsize,aggr%inner%nagr)) 02576 allocate(nneig(aggr%inner%nagr)) 02577 nnodes = size(aggr%inner%num) 02578 nneig=0 02579 do i=1,A%nnz 02580 if (A%indi(i)>nnodes.or.A%indj(i)>nnodes) cycle 02581 if (A%indi(i)<A%indj(i)) then 02582 c1=aggr%inner%num(A%indi(i)) 02583 c2=aggr%inner%num(A%indj(i)) 02584 if (c1/=c2) then 02585 k=1 02586 do while(k<=nneig(c1)) 02587 if (c2==neigs(k,c1)) then 02588 exit 02589 endif 02590 k=k+1 02591 enddo 02592 if (k>nneig(c1)) then 02593 neigs(k,c1)=c2 02594 nneig(c1)=k 02595 endif 02596 endif 02597 endif 02598 enddo 02599 ! allocation for the adjacency data 02600 allocate(xadj(aggr%inner%nagr+1)) 02601 xadj = 0 02602 ! count the lengths 02603 ! (NB! We are expecting matrix symmetric structure!!!) 02604 do c1=1,aggr%inner%nagr 02605 do k=1,nneig(c1) 02606 xadj(c1)=xadj(c1)+1 02607 c2=neigs(k,c1) 02608 xadj(c2)=xadj(c2)+1 02609 enddo 02610 enddo 02611 allocate(counter(aggr%inner%nagr)) 02612 counter = 0 02613 s = xadj(1) 02614 xadj(1) = 1 02615 counter(1) = 1 02616 do i = 2,aggr%inner%nagr 02617 s1 = xadj(i) 02618 xadj(i) = xadj(i-1)+s 02619 counter(i) = xadj(i) 02620 s = s1 02621 enddo 02622 xadj(aggr%inner%nagr+1) = xadj(aggr%inner%nagr) + s 02623 sadjncy = xadj(aggr%inner%nagr+1) - 1 02624 allocate(adjncy(sadjncy)) 02625 ! pass 2 of the data 02626 do c1=1,aggr%inner%nagr 02627 do k=1,nneig(c1) 02628 c2=neigs(k,c1) 02629 adjncy(counter(c1))=c2 02630 counter(c1)=counter(c1)+1 02631 adjncy(counter(c2))=c1 02632 counter(c2)=counter(c2)+1 02633 enddo 02634 enddo 02635 nedges = sadjncy/2 02636 deallocate(counter) 02637 deallocate(nneig) 02638 deallocate(neigs) 02639 end subroutine SpMtx_buildAggrAdjncy 02640 02641 subroutine SpMtx_SymmTest(A,eps) 02642 type(SpMtx),intent(in) :: A 02643 real(kind=rk),optional :: eps 02644 type(SpMtx) :: T 02645 integer :: i,ii,j,k,c 02646 logical :: found 02647 real(kind=rk) :: epsil=1.0E-15_rk 02648 real(kind=rk) :: aa,delta 02649 02650 if (present(eps)) then 02651 epsil=eps 02652 endif 02653 c=0 02654 T=SpMtx_Copy(A) 02655 call SpMtx_arrange(T,D_SpMtx_ARRNG_COLS) 02656 do k=1,A%nnz 02657 i=A%indi(k) 02658 j=A%indj(k) 02659 found=.false. 02660 doing:do ii=T%M_bound(i),T%M_bound(i+1)-1 02661 if (T%indj(ii)==i.and.T%indi(ii)==j) then ! found the transp. loc. 02662 found=.true. 02663 exit doing 02664 endif 02665 enddo doing 02666 if (.not.found) then 02667 write(stream,*)'SymmTest: element matching ',i,j,' not found!' 02668 c=c+1 02669 else 02670 aa=max(abs(A%val(k)),abs(T%val(k))) 02671 delta=abs(A%val(k)-T%val(ii)) 02672 if (delta>epsil.and.delta/aa>epsil) then 02673 write(stream,*)'SymmTest::',i,j, & 02674 delta,A%val(k),T%val(ii),delta,delta/aa 02675 c=c+1 02676 endif 02677 endif 02678 if (c>20) then 02679 write(stream,*)'Stopping after error count ',c 02680 stop 02681 endif 02682 enddo 02683 call SpMtx_Destroy(T) 02684 if (c==0) then 02685 write(stream,*) 'SymmTest successful with epsil ',epsil 02686 endif 02687 end subroutine SpMtx_SymmTest 02688 02689 !------------------------------------------------------ 02690 End Module SpMtx_arrangement 02691 !------------------------------------------------------ 02692
1.7.3-20110217