DOUG 0.2

CreateRestrict.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 
00043 module CoarseCreateRestrict
00044     use RealKind
00045 
00046 #include<doug_config.h>
00047 
00048 #ifdef D_COMPLEX
00049 #define float complex
00050 #else
00051 #define float real
00052 #endif
00053     ! For now
00054     real(kind=xyzk), parameter :: invdistpow=0.5_xyzk
00055 
00056 contains
00057 
00058     !! Create the Restriction matrix
00059     subroutine CreateRestrict(C,M,R)
00060         use CoarseGrid_class
00061         use GeomInterp
00062         use SpMtx_class
00063         use CoarseMtx_mod
00064         use globals
00065  
00066         implicit none
00067 
00068         !! Coarse Grid whose structure to use
00069         type(CoarseGrid), intent(inout) :: C
00070         !! Fine Mesh for which to build
00071         type(Mesh), intent(in) :: M
00072         !! Restriction matrix to create
00073         type(SpMtx), intent(out) :: R
00074 
00075   
00076         ! We can add a choice here to do different things
00077 
00078         select case(sctls%interpolation_type)
00079         
00080         case (COARSE_INTERPOLATION_INVDIST)
00081         call CreatePathRestrict(C,M,R,genNoData,&
00082                                 getInvDistVals,getSizeOne)
00083 
00084         case (COARSE_INTERPOLATION_KRIGING)
00085         call CreatePathRestrict(C,M,R,genKrigingData,&
00086                                 getKrigingVals,getKrigingSize)
00087                                 
00088         case (COARSE_INTERPOLATION_MULTLIN)
00089         !call CreatePathRestrict(C,M,Restrict,genMultiLinearData,&
00090         !                        getMultiLinearVals,getMultiLinearSize)
00091         call CreateGeneralRestrict(C,M,R,CalcMlinearInterp)
00092 
00093         end select
00094 
00095     end subroutine
00096 
00097     !! Calculates the node prolongation matrix taking the coarse element path
00098     !!   down to the fine node. Later rebuilds it into a freedom restrict mat.
00099     subroutine CreatePathRestrict(C,M,R,gendata,getvals,getsize)
00100         use RealKind
00101         use CoarseGrid_class
00102         use SpMtx_class
00103         use SpMtx_util
00104         use globals
00105         
00106         implicit none
00107 
00108         !! Coarse Grid whose structure to use
00109         type(CoarseGrid), intent(inout) :: C
00110         !! Fine Mesh for which to build
00111         type(Mesh), intent(in) :: M
00112         !! The restriction matrix to be created
00113         type(SpMtx), intent(out) :: R
00114   
00115         !! gendata is used to create aux data of size getsize 
00116         !!  for interpolation from points in pts
00117         interface
00118             subroutine gendata(cpts,csz,nsd,routp,ioutp)
00119                 use RealKind
00120                 real(kind=xyzk), intent(in) :: cpts(:,:) ! pts(nsd,csz)
00121                 integer :: csz, nsd
00122                 real(kind=xyzk), intent(out) :: routp(:) ! routp(undet.)
00123                 integer, intent(out) :: ioutp(:) ! ioutp(undet.)
00124             end subroutine gendata
00125         end interface
00126 
00127         !! getvals is used to create the prolongation data for a fine point
00128         interface
00129             subroutine getvals(cpts,csz,nsd,rinp,iinp,pt,outp)
00130                 use RealKind
00131                 real(kind=xyzk), intent(in) :: cpts(:,:) ! pts(nsd,csz)
00132                 integer :: csz, nsd
00133                 real(kind=xyzk), intent(in) :: rinp(:) ! rinp(undet.)
00134                 integer, intent(in) :: iinp(:) ! iinp(undet.)
00135                 real(kind=xyzk), intent(in) :: pt(:) ! pt(nsd)
00136                 float(kind=rk), intent(out) :: outp(:) ! outp(csz)
00137             end subroutine getvals
00138         end interface
00139 
00140         !! get the max size of aux data needed for ptcnt points
00141         interface
00142             subroutine getsize(ptcnt,nsd,rsz,isz)
00143                 integer, intent(in) :: ptcnt
00144                 integer, intent(in) :: nsd
00145                 integer, intent(out) :: rsz, isz
00146             end subroutine getsize
00147         end interface
00148        
00149        
00150         type(SpMtx) :: NP ! node prolongation matrix - restrict is transposed
00151        
00152         integer :: cnt
00153         integer :: i, j, k, pn, nd, f
00154         integer :: tnsd
00155         real(kind=xyzk) :: pts(M%nsd,2**M%nsd+C%mlvl)
00156         integer :: pinds(2**M%nsd+C%mlvl)
00157         real(kind=xyzk),pointer :: pt(:)
00158         real(kind=xyzk),pointer :: raux(:)
00159         integer,pointer :: iaux(:)
00160 
00161         ! For the reverse of freemap
00162         integer :: flist(C%nlfc)
00163         integer :: fbeg(C%nct+1)
00164         integer :: faux(C%nct)
00165 
00166         tnsd=2**M%nsd
00167 
00168         !****************************************************
00169         ! Create the Node Prolongation matrix structure
00170         !****************************************************
00171 
00172         ! Calculate nnz
00173         
00174         cnt=0
00175         
00176         ! Increment by initial grid
00177         do i=1, C%elnum
00178             cnt=cnt+tnsd*C%els(i)%nfs
00179         enddo
00180 
00181         ! Increment by refinements
00182         do i=1, C%refnum
00183             cnt = cnt + &
00184                      C%refels(i)%level*(C%refels(i)%lend-C%refels(i)%lbeg+1)
00185         enddo
00186 
00187         ! Allocate the node prolong matrix
00188         NP=SpMtx_newInit(cnt,nrows=M%lnnode,ncols=C%nct)
00189         allocate(NP%M_bound(M%lnnode+1))
00190 
00191         ! Create NP%M_bound - done in 2 phases
00192          NP%M_bound=0
00193        
00194         ! Vector in coarse non-refined elements
00195         do i=1, C%elnum
00196             if (C%els(i)%rbeg==-1) then
00197             do j=C%els(i)%lbeg,C%els(i)%lbeg+C%els(i)%nfs-1
00198                 NP%M_bound(C%elmap(j)+1)=tnsd
00199             enddo
00200             endif
00201         enddo
00202 
00203         !  Vector in refined elements
00204         do i=1, C%refnum
00205            cnt=tnsd+C%refels(i)%level
00206            do j=C%refels(i)%lbeg,C%refels(i)%lend
00207                 NP%M_bound(C%elmap(j)+1)=cnt
00208            enddo
00209         enddo
00210         
00211         ! And add the counts together to get the final M_bound
00212         NP%M_bound(1)=1
00213         do i=1,M%lnnode
00214             NP%M_bound(i+1)=NP%M_bound(i)+NP%M_bound(i+1)
00215         enddo
00216 
00217         !****************************************************
00218         ! Create the prolongation data for NP
00219         !****************************************************
00220 
00221         ! Init the aux buffer to max possible size
00222         call getsize(C%mlvl+tnsd,M%nsd,i, k)
00223         allocate(raux(i),iaux(k))
00224 
00225         do i=1,C%elnum
00226         
00227             ! Gather the elements corner points to an array
00228             do k=1,2**M%nsd
00229                 pts(:,k)=C%coords(:,C%els(i)%n(k))
00230             enddo
00231            
00232             ! If not refined loop over all the nodes contained in the element
00233             if (C%els(i)%rbeg==-1) then
00234                 ! Generate the auxilliary data
00235                 call gendata(pts,tnsd,M%nsd,raux,iaux)
00236                 pinds(1:tnsd)=C%els(i)%n
00237 
00238                 do j=C%els(i)%lbeg,C%els(i)%lbeg+C%els(i)%nfs-1
00239                     pn=C%elmap(j) ! The index of the fine node
00240 
00241                     ! Evaluate the functions in the point
00242                     call getvals(pts,tnsd,M%nsd,raux,iaux,M%lcoords(:,pn), &
00243                              NP%val(NP%M_bound(pn):NP%M_bound(pn+1)-1))
00244 
00245 !                    write(stream,*) "NP",pn,M%lnnode,NP%M_bound(pn),NP%nnz
00246                 
00247                     ! Set appropriate indi and indj
00248                     NP%indi(NP%M_bound(pn):NP%M_bound(pn+1)-1)=pn
00249                     NP%indj(NP%M_bound(pn):NP%M_bound(pn+1)-1)=pinds
00250                 enddo
00251             else ! Otherwise loop through the refinements
00252             pinds(1:tnsd)=C%els(i)%n
00253             j=C%els(i)%rbeg
00254             do while (j/=-1)
00255                 ! Set the node of the element into pts
00256                 pts(:,tnsd+C%refels(j)%level)=C%coords(:,C%refels(j)%node)
00257                 pinds(tnsd+C%refels(j)%level)=C%refels(j)%node
00258 
00259                 ! Check if we actually need to do anything here
00260                 if (C%refels(j)%lbeg<=C%refels(j)%lend) then
00261                     call gendata(pts,tnsd+C%refels(j)%level,M%nsd,raux,iaux)
00262 
00263                     ! Go through all the fine nodes that are in this element
00264                     do k=C%refels(j)%lbeg,C%refels(j)%lend
00265                         pn=C%elmap(k)
00266  
00267                         ! Evaluate the functions in the point
00268                         call getvals(pts,tnsd+C%refels(j)%level,&
00269                                 M%nsd,raux,iaux,M%lcoords(:,pn), &
00270                                 NP%val(NP%M_bound(pn):NP%M_bound(pn+1)-1))
00271                 
00272                         ! Set appropriate indi and indj
00273                         NP%indi(NP%M_bound(pn):NP%M_bound(pn+1)-1)=pn
00274                         NP%indj(NP%M_bound(pn):NP%M_bound(pn+1)-1)=pinds
00275                      enddo
00276                 endif
00277 
00278                 ! And go to the next refined element
00279                 j=C%refels(j)%next
00280             enddo
00281             endif
00282         enddo
00283 
00284         deallocate(raux,iaux)
00285         
00286 
00287         !****************************************************
00288         ! Create the reverse of freemap
00289         !****************************************************
00290 
00291         ! Find the counts of freedoms associated with nodes
00292         faux=0;
00293         do i=1,C%nlfc
00294             faux(C%cfreemap(i))=faux(C%cfreemap(i))+1
00295         enddo
00296 
00297         ! Create fbegs
00298         fbeg(1)=1;
00299         do i=1,C%nct
00300             fbeg(i+1)=faux(i)+fbeg(i)
00301         enddo
00302         
00303         ! Create flist
00304         faux=fbeg(1:C%nlfc)
00305         do i=1,C%nlfc
00306             nd=C%cfreemap(i)
00307             flist(faux(nd))=i
00308             faux(nd)=faux(nd)+1
00309         enddo
00310 
00311         !****************************************************
00312         ! Create R based on NP
00313         !****************************************************
00314 
00315         ! Calculate the upper bound for nnz
00316         cnt=0
00317         do i=1,NP%nnz
00318             cnt=cnt+fbeg(NP%indj(i)+1)-fbeg(NP%indj(i))
00319         enddo
00320         
00321         ! Allocate the matrix
00322         R=SpMtx_newInit(cnt,nrows=C%nlfc,ncols=M%nlf)
00323         allocate(R%M_bound(M%nlf+1))
00324 
00325         ! Fill it
00326         f=1;         
00327         do i=1,M%nlf    
00328             nd=M%lfreemap(i)
00329             R%M_bound(i)=f
00330             do j=NP%M_bound(nd),NP%M_bound(nd+1)-1
00331                 do k=fbeg(NP%indj(j)),fbeg(NP%indj(j)+1)-1
00332                 if (NP%val(j)>eps) then
00333                     R%val(f)=NP%val(j)
00334                     R%indi(f)=flist(k)
00335                     R%indj(f)=i
00336                     f=f+1
00337                 endif
00338                 enddo
00339             enddo
00340         enddo
00341         R%M_bound(M%nlf+1)=f
00342         
00343         ! Set the matrix to be row-sorted format
00344         R%arrange_type=D_SpMtx_ARRNG_COLS
00345 
00346         ! Destroy the node version
00347         call SpMtx_Destroy(NP)
00348         
00349         ! Give back the unneeded memory in P
00350         if (R%nnz>f-1) call SpMtx_resize(R,f-1)
00351         
00352         if (sctls%verbose>3) &
00353              write (stream,*) "Restriction matrix has ",f-1," elements"
00354 
00355    end subroutine CreatePathRestrict
00356 
00357     !************************************************
00358     ! Inverse Distances Interpolation
00359     !************************************************
00360     
00361     !! A function to pass to CreateProlong to return 1
00362     subroutine getSizeOne(ptcnt, nsd, rsz, isz)
00363         integer, intent(in) :: ptcnt, nsd
00364         integer, intent(out) :: rsz, isz
00365         rsz=1; isz=1;
00366     end subroutine getSizeOne
00367  
00368     !! A function to pass to CreateProlong to create no real data
00369     subroutine genNoData(cpts,csz,nsd,routp,ioutp)
00370         use RealKind
00371         real(kind=xyzk), intent(in) :: cpts(:,:) ! pts(nsd,csz)
00372         integer :: csz, nsd
00373         real(kind=xyzk), intent(out) :: routp(:) ! outp(undet.)
00374         integer, intent(out) :: ioutp(:)
00375         routp(1)=0.0_rk; ioutp(1)=0
00376     end subroutine genNoData
00377 
00378     !! Get the interpolation values by using inverse distances method
00379     subroutine getInvDistVals(cpts,csz,nsd,rinp,iinp,pt,outp)
00380         use RealKind
00381         use globals, only: eps
00382         real(kind=xyzk), intent(in) :: cpts(:,:) ! pts(nsd,csz)
00383         integer :: csz, nsd
00384         real(kind=xyzk), intent(in) :: rinp(:) ! rinp(1)
00385         integer, intent(in) :: iinp(:) ! iinp(1)
00386         real(kind=xyzk), intent(in) :: pt(:) ! pt(nsd)
00387         float(kind=rk), intent(out) :: outp(:) ! outp(csz)
00388 
00389         integer :: i
00390 
00391         ! Get the distances
00392         do i=1,csz
00393             outp(i)=(dot_product((cpts(:,i)-pt),(cpts(:,i)-pt)))**invdistpow
00394 
00395             if (outp(i)<=eps) then
00396                 outp(1:csz)=0.0_rk;
00397                 outp(i)=1.0_rk;
00398                 exit
00399             else
00400                 outp(i)=1.0_rk/outp(i)
00401             endif
00402         enddo
00403 
00404         ! Normalize to preserve constant
00405         outp(1:csz)=outp(1:csz)/sum(outp(1:csz))
00406         
00407     end subroutine getInvDistVals
00408     
00409     !************************************************
00410     ! Kriging Interpolation
00411     !************************************************
00412     
00413     !! A function to pass to CreateProlong for Kriging
00414     subroutine getKrigingSize(ptcnt, nsd, rsz, isz)
00415         integer, intent(in) :: ptcnt, nsd
00416         integer, intent(out) :: rsz, isz
00417         rsz=(ptcnt+1)**2 ! Matrix
00418         isz=ptcnt*2 ! Pivot data
00419     end subroutine getKrigingSize
00420  
00421     !! A function to pass to CreateProlong to create the Kriging aux data
00422     subroutine genKrigingData(cpts,csz,nsd,routp,ioutp)
00423         use RealKind
00424         real(kind=xyzk), intent(in) :: cpts(:,:) ! pts(nsd,csz)
00425         integer :: csz, nsd
00426         real(kind=xyzk), intent(out) :: routp(:) ! routp((csz+1)**2)
00427         integer, intent(out) :: ioutp(:) ! ioutp(2*csz)
00428        
00429         real(kind=xyzk) :: mat(csz+1,csz+1)
00430         integer :: i, k
00431         
00432         ! Create the matrix
00433         do i=1,csz
00434         mat(csz,i)=1.0_xyzk       ! Last row
00435         mat(i,csz)=1.0_xyzk       ! Last column
00436         mat(i,i)=0.0_xyzk         ! Diagonal
00437 
00438         ! Other elements
00439         do k=i+1,csz
00440         ! Set the distance
00441         mat(i,k)= &
00442          (dot_product((cpts(:,i)-cpts(:,k)),(cpts(:,i)-cpts(:,k))))**invdistpow
00443         ! Copy it
00444         mat(k,i)=mat(i,k)
00445         enddo
00446         enddo
00447 
00448         ! LUP-factorize mat
00449 
00450         !!!! call DGETRF(csz+1,csz+1,mat,csz+1,ioutp,i)
00451 
00452         ! if (i/=0) ERROR
00453 
00454         ! Pack mat into outp
00455         routp(1:(csz+1)**2)=reshape(mat,(/ (csz+1)**2 /) )
00456 
00457     end subroutine genKrigingData
00458 
00459     !! Get the interpolation values by using Kriging method
00460     subroutine getKrigingVals(cpts,csz,nsd,rinp,iinp,pt,outp)
00461         use RealKind
00462         real(kind=xyzk), intent(in) :: cpts(:,:) ! pts(nsd,csz)
00463         integer :: csz, nsd
00464         float(kind=rk), intent(in) :: rinp(:) ! rinp((csz+1)**2)
00465         integer, intent(in) :: iinp(:) ! iinp(2*csz)
00466         real(kind=xyzk), intent(in) :: pt(:) ! pt(nsd)
00467         float(kind=rk), intent(out) :: outp(:) ! outp(csz)
00468 
00469         float(kind=rk) :: vec(csz+1,1)
00470         integer :: i
00471         
00472         ! Calculate the distances vector for this point
00473         do i=1,csz
00474             vec(i,1)=(dot_product((cpts(:,i)-pt),(cpts(:,i)-pt)))**invdistpow
00475         enddo
00476         vec(csz+1,1)=1.0_xyzk
00477 
00478         ! Solve it
00479         !!!!call DGETRS('N', csz+1, 1, reshape(rinp,(/ csz+1, csz+1 /)),csz+1, &
00480         !!!!                iinp, vec, csz+1, i)
00481 
00482         ! if (i/=0) ERROR
00483 
00484         ! Copy the data over
00485         outp(1:csz)=vec(1:csz,1)
00486         
00487     end subroutine getKrigingVals
00488     
00489 
00490     !! Strip the restriction matrix
00491     !! useful for multiple processor case, where 
00492     !! the stripped variant is used for vector restrict/interpolate
00493     subroutine stripRestrict(M, R)
00494         use Mesh_class
00495         use SpMtx_class
00496         use globals, only : stream
00497 
00498         type(Mesh), intent(in) :: M
00499         type(SpMtx), intent(inout) :: R
00500 !        type(SpMtx), intent(out) :: Ra
00501 
00502         integer :: i,j,cnt
00503         integer :: unq(M%nlf)
00504 
00505         ! Calculate unq (has 1 if that node is "unique"
00506         ! a.k.a. this is the lowest ranked process which has it)
00507         unq=0; cnt=0
00508         do i=1,M%nell
00509            if (M%eptnmap(i)<=myrank) then ! lower ranked
00510            do j=1,M%nfrelt(i)  
00511                 if (M%gl_fmap(M%mhead(j,i))/=0) & ! non-unique
00512                         unq(M%gl_fmap(M%mhead(j,i)))=-1
00513            enddo
00514            elseif (M%eptnmap(i)==myrank+1) then ! my nodes
00515            do j=1,M%nfrelt(i)
00516                 if (unq(M%gl_fmap(M%mhead(j,i)))==0) & ! untagged
00517                         unq(M%gl_fmap(M%mhead(j,i)))=1
00518            enddo
00519            endif
00520         enddo
00521         
00522         ! Count the nnz in the stripped restrict matrix
00523 !        cnt=0
00524 !        if (R%arrange_type==D_SpMtx_ARRNG_COLS) then ! can use M_bound
00525 !            do i=1,R%ncols
00526 !                if (unq(i)==1) cnt=cnt+R%M_bound(i+1)-R%M_bound(i)
00527 !            enddo
00528 !        else ! count the matrix elements one by one
00529 !            do i=1,R%nnz
00530 !                if (unq(R%indj(i))==1) cnt=cnt+1
00531 !            enddo
00532 !        endif
00533 
00534 !        ! Initalize Ra       
00535 !        Ra=SpMtx_newInit(nnz=cnt,nrows=R%nrows,ncols=R%ncols)
00536 
00537         ! if we have something to actually do
00538 !        if (cnt/=0) then
00539             ! Fille Ra with its elements
00540             j=1
00541             do i=1,R%nnz
00542                 if (unq(R%indj(i))==1) then ! it stays
00543                     R%indi(j)=R%indi(i)
00544                     R%indj(j)=R%indj(i)
00545                     R%val(j)=R%val(i)
00546                     j=j+1
00547                 endif
00548             enddo
00549 !        endif
00550 
00551         call SpMtx_resize(R,j-1)
00552         R%arrange_type=D_SpMtx_ARRNG_NO
00553         deallocate(R%M_bound)
00554 
00555         if (sctls%verbose>3) &
00556              write (stream,*) "Stripped restriction matrix has ",j-1," elements"
00557 
00558     end subroutine stripRestrict
00559  
00560 end module CoarseCreateRestrict