DOUG 0.2

SpMtx_arrangement.F90

Go to the documentation of this file.
00001 ! DOUG - Domain decomposition On Unstructured Grids
00002 ! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
00003 ! Department of Mathematics, University of Bath
00004 !
00005 ! This library is free software; you can redistribute it and/or
00006 ! modify it under the terms of the GNU Lesser General Public
00007 ! License as published by the Free Software Foundation; either
00008 ! version 2.1 of the License, or (at your option) any later version.
00009 !
00010 ! This library is distributed in the hope that it will be useful,
00011 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013 ! Lesser General Public License for more details.
00014 !
00015 ! You should have received a copy of the GNU Lesser General Public
00016 ! License along with this library; if not, write to the Free Software
00017 ! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00018 ! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
00019 ! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
00020 ! mailto:info(at)dougdevel.org)
00021 
00022 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