[cig-commits] r12913 - seismo/3D/SPECFEM3D_GLOBE/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Sep 17 12:37:13 PDT 2008
Author: dkomati1
Date: 2008-09-17 12:37:13 -0700 (Wed, 17 Sep 2008)
New Revision: 12913
Removed:
seismo/3D/SPECFEM3D_GLOBE/trunk/get_perm_cuthill_mckee.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90
Log:
suppressed some old routines that used to be called for Cuthill-McKee sorting in previous versions
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90 2008-09-17 19:17:36 UTC (rev 12912)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90 2008-09-17 19:37:13 UTC (rev 12913)
@@ -426,91 +426,3 @@
end subroutine compute_jacobian_2D
-
-
-subroutine sort_arrays_for_cuthill (ispecb,xstore,ystore,zstore,ibelm,normal,jacobian2D,nspec2D,NGLL1,NGLL2,nspec)
-
- implicit none
-
- include "constants.h"
-
- integer :: ispecb,nspec2D,NGLL1,NGLL2,nspec,ispec_tmp,dummy_var,i
-
- integer ibelm(nspec2D)
- real(kind=CUSTOM_REAL) jacobian2D(NGLL1,NGLL2,NSPEC2D)
- real(kind=CUSTOM_REAL) normal(NDIM,NGLL1,NGLL2,NSPEC2D)
-
- double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
-
-! arrays for sorting routine
- integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
- logical, dimension(:), allocatable :: ifseg
- double precision, dimension(:), allocatable :: work
- double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
- integer, dimension(:), allocatable :: perm
- integer, dimension(:), allocatable :: ibelm_tmp
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: normal_tmp
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: jacobian2D_tmp
-
-! get permutation
- allocate (xstore_selected(ispecb))
- allocate (ystore_selected(ispecb))
- allocate (zstore_selected(ispecb))
- allocate(ind(ispecb))
- allocate(ninseg(ispecb))
- allocate(iglob(ispecb))
- allocate(locval(ispecb))
- allocate(ifseg(ispecb))
- allocate(iwork(ispecb))
- allocate(work(ispecb))
- allocate(perm(ispecb))
-
- do ispec_tmp=1,ispecb
- xstore_selected(ispec_tmp) = xstore(1,1,1,ibelm(ispec_tmp))
- ystore_selected(ispec_tmp) = ystore(1,1,1,ibelm(ispec_tmp))
- zstore_selected(ispec_tmp) = zstore(1,1,1,ibelm(ispec_tmp))
- perm(ispec_tmp) = ispec_tmp
- enddo
-
- call sort_array_coordinates(ispecb,xstore_selected,ystore_selected,zstore_selected, &
- perm,iglob,locval,ifseg,dummy_var,ind,ninseg,iwork,work)
-
- deallocate (xstore_selected)
- deallocate (ystore_selected)
- deallocate (zstore_selected)
- deallocate(ind)
- deallocate(ninseg)
- deallocate(iglob)
- deallocate(locval)
- deallocate(ifseg)
- deallocate(iwork)
- deallocate(work)
-
-! permutation of ibelm
- allocate(ibelm_tmp(ispecb))
- ibelm_tmp(1:ispecb) = ibelm(1:ispecb)
- do i = 1,ispecb
- ibelm(i) = ibelm_tmp(perm(i))
- enddo
- deallocate(ibelm_tmp)
-
-! permutation of normal
- allocate(normal_tmp(NDIM,NGLL1,NGLL2,ispecb))
- normal_tmp(:,:,:,1:ispecb) = normal(:,:,:,1:ispecb)
- do i = 1,ispecb
- normal(:,:,:,i) = normal_tmp(:,:,:,perm(i))
- enddo
- deallocate(normal_tmp)
-
-! permutation of jacobian2D
- allocate(jacobian2D_tmp(NGLL1,NGLL2,ispecb))
- jacobian2D_tmp(:,:,1:ispecb) = jacobian2D(:,:,1:ispecb)
- do i = 1,ispecb
- jacobian2D(:,:,i) = jacobian2D_tmp(:,:,perm(i))
- enddo
- deallocate(jacobian2D_tmp)
- deallocate(perm)
-
-end subroutine sort_arrays_for_cuthill
Deleted: seismo/3D/SPECFEM3D_GLOBE/trunk/get_perm_cuthill_mckee.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_perm_cuthill_mckee.f90 2008-09-17 19:17:36 UTC (rev 12912)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_perm_cuthill_mckee.f90 2008-09-17 19:37:13 UTC (rev 12913)
@@ -1,789 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
-! --------------------------------------------------
-!
-! Main authors: Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory, California Institute of Technology, USA
-! and University of Pau / CNRS / INRIA, France
-! (c) California Institute of Technology and University of Pau / CNRS / INRIA
-! February 2008
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
-! implement reverse Cuthill-McKee (1969) ordering, introduced in
-! E. Cuthill and J. McKee. Reducing the bandwidth of sparse symmetric matrices.
-! In Proceedings of the 1969 24th national conference, pages 157-172,
-! New-York, New-York, USA, 1969. ACM Press.
-! see for instance http://en.wikipedia.org/wiki/Cuthill%E2%80%93McKee_algorithm
-
- subroutine get_perm(ibool,perm,limit,nspec,nglob,INVERSE,FACE)
-
- implicit none
-
- include "constants.h"
-
- logical :: INVERSE,FACE
-
-! input
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
-
-! output
- integer, dimension(nspec) :: perm
-
-! local variables
- integer nspec,nglob_GLL_full
-
-! maximum number of neighbors of a spectral element (in principle, it could be any value)
- integer, parameter :: MAX_NUMBER_OF_NEIGHBORS = 100
-
-! global corner numbers that need to be created
- integer, dimension(nglob) :: global_corner_number
-
- integer mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1)
- integer, dimension(:), allocatable :: ne,np,adj
- integer xadj(nspec+1)
-
-! arrays to store the permutation and inverse permutation of the Cuthill-McKee algorithm
- integer, dimension(nspec) :: invperm
-
- logical maskel(nspec)
-
- integer i,istart,istop,number_of_neighbors
-
- integer nglob_eight_corners_only,nglob
-
-! only count the total size of the array that will be created, or actually create it
- logical count_only
- integer total_size_ne,total_size_adj,limit
-
-!
-!-----------------------------------------------------------------------
-!
- if(PERFORM_CUTHILL_MCKEE) then
-
- ! total number of points in the mesh
- nglob_GLL_full = nglob
-
- !---- call Charbel Farhat's routines
- call form_elt_connectivity_foelco(mn,mp,nspec,global_corner_number,nglob_GLL_full,ibool,nglob_eight_corners_only)
- do i=1,nspec
- istart = mp(i)
- istop = mp(i+1) - 1
- enddo
-
- allocate(np(nglob_eight_corners_only+1))
- count_only = .true.
- total_size_ne = 1
- allocate(ne(total_size_ne))
- call form_node_connectivity_fonoco(mn,mp,ne,np,nglob_eight_corners_only,nspec,count_only,total_size_ne)
- deallocate(ne)
- allocate(ne(total_size_ne))
- count_only = .false.
- call form_node_connectivity_fonoco(mn,mp,ne,np,nglob_eight_corners_only,nspec,count_only,total_size_ne)
- do i=1,nglob_eight_corners_only
- istart = np(i)
- istop = np(i+1) - 1
- enddo
-
- count_only = .true.
- total_size_adj = 1
- allocate(adj(total_size_adj))
- call create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec,nglob_eight_corners_only,&
- count_only,total_size_ne,total_size_adj,FACE)
- deallocate(adj)
- allocate(adj(total_size_adj))
- count_only = .false.
- call create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec,nglob_eight_corners_only,&
- count_only,total_size_ne,total_size_adj,FACE)
- do i=1,nspec
- istart = xadj(i)
- istop = xadj(i+1) - 1
- number_of_neighbors = istop-istart+1
- if(number_of_neighbors < 1) stop 'error: your mesh seems to have at least one element not connected to any other'
- if(number_of_neighbors > MAX_NUMBER_OF_NEIGHBORS) stop 'error: your mesh seems to have an unlikely high valence'
- enddo
- deallocate(ne,np)
-
-! call the Cuthill-McKee sorting algorithm
- call cuthill_mckee(adj,xadj,perm,invperm,nspec,total_size_adj,limit,INVERSE)
- deallocate(adj)
- else
-! create identity permutation in order to do nothing
- do i=1,nspec
- perm(i) = i
- enddo
- endif
-
- end subroutine get_perm
-
-!=======================================================================
-!
-! Charbel Farhat's FEM topology routines
-!
-! Dimitri Komatitsch, February 1996 - Code based on Farhat's original version
-! described in his technical report from 1987
-!
-! modified and adapted by Dimitri Komatitsch, May 2006
-!
-!=======================================================================
-
- subroutine form_elt_connectivity_foelco(mn,mp,nspec,global_corner_number, &
- nglob_GLL_full,ibool,nglob_eight_corners_only)
-
-!-----------------------------------------------------------------------
-!
-! Forms the MN and MP arrays
-!
-! Input :
-! -------
-! ibool Array needed to build the element connectivity table
-! nspec Number of elements in the domain
-! NGNOD_HEXAHEDRA number of nodes per hexahedron (brick with 8 corners)
-!
-! Output :
-! --------
-! MN, MP This is the element connectivity array pair.
-! Array MN contains the list of the element
-! connectivity, that is, the nodes contained in each
-! element. They are stored in a stacked fashion.
-!
-! Pointer array MP stores the location of each
-! element list. Its length is equal to the number
-! of elements plus one.
-!
-!-----------------------------------------------------------------------
-
- implicit none
-
- include "constants.h"
-
- integer nspec,nglob_GLL_full
-
-! arrays with mesh parameters per slice
- integer, intent(in), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
-
-! global corner numbers that need to be created
- integer, intent(out), dimension(nglob_GLL_full) :: global_corner_number
- integer, intent(out) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1)
- integer, intent(out) :: nglob_eight_corners_only
-
- integer ninter,nsum,ispec,node,k,inumcorner,ix,iy,iz
-
- ninter = 1
- nsum = 1
- mp(1) = 1
-
-!---- define topology of the elements in the mesh
-!---- we need to define adjacent numbers from the sub-mesh consisting of the corners only
- nglob_eight_corners_only = 0
- global_corner_number(:) = -1
-
- do ispec=1,nspec
-
- inumcorner = 0
- do iz = 1,NGLLZ,NGLLZ-1
- do iy = 1,NGLLY,NGLLY-1
- do ix = 1,NGLLX,NGLLX-1
-
- inumcorner = inumcorner + 1
- if(inumcorner > NGNOD_HEXAHEDRA) stop 'corner number too large'
-
-! check if this point was already assigned a number previously, otherwise create one and store it
- if(global_corner_number(ibool(ix,iy,iz,ispec)) == -1) then
- nglob_eight_corners_only = nglob_eight_corners_only + 1
- global_corner_number(ibool(ix,iy,iz,ispec)) = nglob_eight_corners_only
- endif
-
- node = global_corner_number(ibool(ix,iy,iz,ispec))
- do k=nsum,ninter-1
- if(node == mn(k)) goto 200
- enddo
-
- mn(ninter) = node
- ninter = ninter + 1
- 200 continue
-
- enddo
- enddo
- enddo
-
- nsum = ninter
- mp(ispec + 1) = nsum
-
- enddo
-
- end subroutine form_elt_connectivity_foelco
-
-!
-!----------------------------------------------------
-!
-
- subroutine form_node_connectivity_fonoco(mn,mp,ne,np,nglob_eight_corners_only, &
- nspec,count_only,total_size_ne)
-
-!-----------------------------------------------------------------------
-!
-! Forms the NE and NP arrays
-!
-! Input :
-! -------
-! MN, MP, nspec
-! nglob_eight_corners_only Number of nodes in the domain
-!
-! Output :
-! --------
-! NE, NP This is the node-connected element array pair.
-! Integer array NE contains a list of the
-! elements connected to each node, stored in stacked fashion.
-!
-! Array NP is the pointer array for the
-! location of a node's element list in the NE array.
-! Its length is equal to the number of points plus one.
-!
-!-----------------------------------------------------------------------
-
- implicit none
-
- include "constants.h"
-
-! only count the total size of the array that will be created, or actually create it
- logical count_only
- integer total_size_ne
-
- integer nglob_eight_corners_only,nspec
-
- integer, intent(in) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1)
-
- integer, intent(out) :: ne(total_size_ne),np(nglob_eight_corners_only+1)
-
- integer nsum,inode,ispec,j
-
- nsum = 1
- np(1) = 1
-
- do inode=1,nglob_eight_corners_only
- do 200 ispec=1,nspec
-
- do j=mp(ispec),mp(ispec + 1) - 1
- if (mn(j) == inode) then
- if(count_only) then
- total_size_ne = nsum
- else
- ne(nsum) = ispec
- endif
- nsum = nsum + 1
- goto 200
- endif
- enddo
- 200 continue
-
- np(inode + 1) = nsum
-
- enddo
-
- end subroutine form_node_connectivity_fonoco
-
-!
-!----------------------------------------------------
-!
-
- subroutine create_adjacency_table_adjncy(mn,mp,ne,np,adj,xadj,maskel,nspec, &
- nglob_eight_corners_only,count_only,total_size_ne,total_size_adj,FACE)
-
-!-----------------------------------------------------------------------
-!
-! Establishes the element adjacency information of the mesh
-! Two elements are considered adjacent if they share a face.
-!
-! Input :
-! -------
-! MN, MP, NE, NP, nspec
-! MASKEL logical mask (length = nspec)
-!
-! Output :
-! --------
-! ADJ, XADJ This is the element adjacency array pair. Array
-! ADJ contains the list of the elements adjacent to
-! element i. They are stored in a stacked fashion.
-! Pointer array XADJ stores the location of each element list.
-!
-!-----------------------------------------------------------------------
-
- implicit none
-
- include "constants.h"
-
-! only count the total size of the array that will be created, or actually create it
- logical count_only,FACE
- integer total_size_ne,total_size_adj
-
- integer nglob_eight_corners_only
-
- integer, intent(in) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1),ne(total_size_ne),np(nglob_eight_corners_only+1)
-
- integer, intent(out) :: adj(total_size_adj),xadj(nspec+1)
-
- logical maskel(nspec)
- integer countel(nspec)
-
- integer nspec,iad,ispec,istart,istop,ino,node,jstart,jstop,nelem,jel
-
- xadj(1) = 1
- iad = 1
-
- do ispec=1,nspec
-
-! reset mask
- maskel(:) = .false.
-
-! mask current element
- maskel(ispec) = .true.
- if (FACE) countel(:) = 0
-
- istart = mp(ispec)
- istop = mp(ispec+1) - 1
- do ino=istart,istop
- node = mn(ino)
- jstart = np(node)
- jstop = np(node + 1) - 1
- do 120 jel=jstart,jstop
- nelem = ne(jel)
- if(maskel(nelem)) goto 120
- if (FACE) then
- ! if 2 elements share at least 3 corners, therefore they share a face
- countel(nelem) = countel(nelem) + 1
- if (countel(nelem)>=3) then
- if(count_only) then
- total_size_adj = iad
- else
- adj(iad) = nelem
- endif
- maskel(nelem) = .true.
- iad = iad + 1
- endif
- else
- if(count_only) then
- total_size_adj = iad
- else
- adj(iad) = nelem
- endif
- maskel(nelem) = .true.
- iad = iad + 1
- endif
- 120 continue
- enddo
-
- xadj(ispec+1) = iad
-
- enddo
-
- end subroutine create_adjacency_table_adjncy
-
-!
-!----------------------------------------------------
-!
-
- subroutine cuthill_mckee(adj,xadj,mask,invperm_all,nspec,total_size_adj,limit,INVERSE)
-
- implicit none
- include "constants.h"
-
- integer, intent(in) :: nspec,total_size_adj, limit
- integer, intent(in) :: adj(total_size_adj),xadj(nspec+1)
-
- integer, intent(out), dimension(nspec) :: mask,invperm_all
- integer, dimension(nspec) :: invperm_sub
- logical :: INVERSE
- integer ispec,gsize,counter,nspec_sub,root,total_ordered_elts, next_root
-
-! fill the mask with ones
- mask(:) = 1
- invperm_all(:) = 0
- counter = 0
- nspec_sub = limit
- root = 1
- total_ordered_elts = 0
-
- do while(total_ordered_elts < nspec)
- ! creation of a sublist of sorted elements which fit in the cache (the criterion of size is limit)
- ! limit = nb of element that can fit in the L2 cache
- call Cut_McK( root, nspec, total_size_adj, xadj, adj, mask, gsize, invperm_sub, limit, nspec_sub, next_root)
- ! add the sublist in the main permutation list
- invperm_all(total_ordered_elts+1:total_ordered_elts+nspec_sub) = invperm_sub(1:nspec_sub)
- total_ordered_elts = total_ordered_elts + nspec_sub
- ! seek for a new root to build the new sublist
- if (next_root > 0) then
- root = next_root
- else
- if (total_ordered_elts /= nspec) &
- call find_next_root(next_root,xadj,adj,total_size_adj,mask,invperm_all,total_ordered_elts,nspec)
- root = next_root
- endif
- enddo
-
- if (INVERSE) then
- do ispec=1,nspec
- mask(invperm_all(ispec)) = ispec
- enddo
- else
- mask(:) = invperm_all(:)
- endif
-
- end subroutine cuthill_mckee
-
-
-!*******************************************************************************
-! Objective: Cuthill-McKee ordering
-! The algorithm is:
-!
-! X(1) = ROOT.
-! for ( I = 1 to N-1)
-! Find all unlabeled neighbors of X(I),
-! assign them the next available labels, in order of increasing degree.
-!
-! Parameters:
-! root the starting point for the cm ordering.
-! nbnodes the number of nodes.
-! nnz the number of adjacency entries.
-!
-! xadj/adj the graph
-! mask only those nodes with nonzero mask are considered
-!
-! gsize the number of the connected component
-! invp Inverse permutation (from new order to old order)
-!*******************************************************************************
-
-subroutine find_next_root(next_root,xadj,adj,total_size_adj,mask,invperm_all,total_ordered_elts,nspec)
-
- implicit none
-
- include "constants.h"
-
-! input
- integer, intent(in) :: total_size_adj,total_ordered_elts,nspec
- integer, intent(in) :: adj(total_size_adj),xadj(nspec+1)
- integer, intent(in), dimension(nspec) :: mask,invperm_all
-! output
- integer, intent(out) :: next_root
-! variables
- integer :: cur_node,neighbor_node,i,j
-
- do i=total_ordered_elts, 1, -1
- cur_node = invperm_all(i)
- do j= xadj(cur_node), xadj(cur_node+1)-1
- neighbor_node = adj(j)
- if (mask(neighbor_node)/=0) then
- next_root=neighbor_node
- return
- endif
- enddo
- enddo
-
-end subroutine find_next_root
-
-!*******************************************************************************
-! Objective: Cuthill-McKee ordering
-! The algorithm is:
-!
-! X(1) = ROOT.
-! for ( I = 1 to N-1)
-! Find all unlabeled neighbors of X(I),
-! assign them the next available labels, in order of increasing degree.
-!
-! Parameters:
-! root the starting point for the cm ordering.
-! nbnodes the number of nodes.
-! nnz the number of adjacency entries.
-!
-! xadj/adj the graph
-! mask only those nodes with nonzero mask are considered
-!
-! gsize the number of the connected component
-! invp Inverse permutation (from new order to old order)
-!*******************************************************************************
-
-subroutine Cut_McK( root, nbnodes, nnz, xadj, adj, mask, gsize, invp, limit, nspec_sub, next_root)
-
- implicit none
-
- include "constants.h"
-
-!--------------------------------------------------------------- Input Variables
- integer root, nnz, nbnodes, limit, nspec_sub, next_root
-
- integer xadj(nbnodes+1), adj(nnz), mask(nbnodes)
-
-!-------------------------------------------------------------- Output Variables
- integer gsize
- integer invp(nbnodes)
-
-!--------------------------------------------------------------- Local Variables
- integer i, j, k, l, lbegin, lnbr, linvp, lvlend, nbr, node, fnbr
- integer deg(nbnodes)
-
-! Find the degrees of the nodes in the subgraph specified by mask and root
-! Here invp is used to store a levelization of the subgraph
- invp(:)=0
- deg(:)=0
- call degree ( root, nbnodes, nnz, xadj, adj, mask, gsize, deg, invp)
-
- mask(root) = 0
-
- IF (gsize > 1) THEN
- !If there is at least 2 nodes in the subgraph
- lvlend = 0
- lnbr = 1
-
- DO while (lvlend < lnbr)
- !lbegin/lvlend point to the begin/end of the present level
- lbegin = lvlend + 1
- lvlend = lnbr
-
- do i= lbegin, lvlend
- node = invp(i)
-
- !Find the unnumbered neighbours of node.
- !fnbr/lnbr point to the first/last neighbors of node
- fnbr = lnbr + 1
- do j= xadj(node), xadj(node+1)-1
- nbr = adj(j)
-
- if (mask(nbr) /= 0) then
- lnbr = lnbr + 1
- mask(nbr) = 0
- invp(lnbr) = nbr
- endif
- enddo
-
- !If no neighbors, go to next node in this level.
- IF (lnbr > fnbr) THEN
- !Sort the neighbors of NODE in increasing order by degree.
- !Linear insertion is used.
- k = fnbr
- do while (k < lnbr)
- l = k
- k = k + 1
- nbr = invp(k)
-
- DO WHILE (fnbr < l)
- linvp = invp(l)
-
- if (deg(linvp) <= deg(nbr)) then
- exit
- endif
-
- invp(l+1) = linvp
- l = l-1
- ENDDO
-
- invp(l+1) = nbr
- enddo
- ENDIF
- enddo
- ENDDO
-
- ENDIF
-
- if (gsize > limit) then
- do i = limit + 1 , nbnodes
- node=invp(i)
- if (node /=0) mask(node) = 1
- enddo
- next_root = invp(limit +1)
- nspec_sub = limit
- else
- next_root = -1
- nspec_sub = gsize
- endif
-
-END subroutine Cut_McK
-
-
-!*******************************************************************************
-! Objective: computes the degrees of the nodes in the connected graph
-!
-! Parameters:
-! root the root node
-! nbnodes the number of nodes in the graph
-! nnz the graph size
-! xadj/adj the whole graph
-! mask Only nodes with mask == 0 are considered
-!
-! gsize the number of nodes in the connected graph
-! deg degree for all the nodes in the connected graph
-! level levelization of the connected graph
-!
-!*******************************************************************************
-
-subroutine degree( root, nbnodes, nnz, xadj, adj, mask, gsize, deg, level )
-
- implicit none
-
-!--------------------------------------------------------------- Input Variables
- integer root, nbnodes, nnz
- integer xadj(nbnodes+1), adj(nnz), mask(nbnodes)
-
-!-------------------------------------------------------------- Output Variables
- integer gsize
- integer deg(nbnodes), level(nbnodes)
-
-!--------------------------------------------------------------- Local Variables
- integer i, j, ideg, lbegin, lvlend, lvsize, nxt, nbr, node
-
-! added a test to detect disconnected subsets in the mesh
-! (in which case Cuthill-McKee fails and should be turned off)
- if(root > nbnodes+1) stop 'error: root > nbnodes+1 in Cuthill-McKee'
- if(root < 1) then
- print *,'error: root < 1 in Cuthill-McKee; you probably have a mesh composed of'
- print *,'two disconnected subsets of elements, in which case Cuthill-McKee fails and should be turned off.'
- print *,'please set PERFORM_CUTHILL_MCKEE = .false. in constants.h and recompile.'
- print *,'please also doublecheck that you indeed want to run two separate meshes simultaneously,'
- print *,'which is extremely unusual (but formally not incorrect).'
- stop 'fatal error in Cuthill-McKee'
- endif
-
-! The sign of xadj(I) is used to indicate if node i has been considered
- xadj(root) = -xadj(root)
- level(1) = root
- nxt = 1
- lvlend = 0
- lvsize = 1
-
- DO WHILE (lvsize > 0)
- ! lbegin/lvlend points the begin/end of the present level
- lbegin = lvlend + 1
- lvlend = nxt
-
- ! Find the degrees of nodes in the present level and generate the next level
- DO i= lbegin, lvlend
- node = level(i)
- ideg = 0
- do j= ABS( xadj(node) ), ABS( xadj(node+1) )-1
- nbr = adj(j)
-
- if (mask(nbr) /= 0) then
- ideg = ideg + 1
-
- if (xadj(nbr) >= 0) then
- xadj(nbr) = -xadj(nbr)
- nxt = nxt + 1
- level(nxt) = nbr
- endif
- endif
- enddo
-
- deg(node) = ideg
- ENDDO
-
- !Compute the level size of the next level
- lvsize = nxt - lvlend
- ENDDO
-
- !Reset xadj to its correct sign
- do i = 1, nxt
- node = level(i)
- xadj(node) = -xadj(node)
- enddo
-
- gsize = nxt
-
-END subroutine degree
-
-!
-!-----------------------------------------------------------------------
-!
-
- subroutine permute_elements_real(array_to_permute,temp_array,perm,nspec)
-
- implicit none
-
- include "constants.h"
-
- integer, intent(in) :: nspec
- integer, intent(in), dimension(nspec) :: perm
-
- real(kind=CUSTOM_REAL), intent(inout), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: array_to_permute,temp_array
-
- integer old_ispec,new_ispec
-
-! copy the original array
- temp_array(:,:,:,:) = array_to_permute(:,:,:,:)
-
- do old_ispec = 1,nspec
- new_ispec = perm(old_ispec)
- array_to_permute(:,:,:,new_ispec) = temp_array(:,:,:,old_ispec)
- enddo
-
- end subroutine permute_elements_real
-
-!
-!-----------------------------------------------------------------------
-!
-
-! implement permutation of elements for arrays of integer type
- subroutine permute_elements_integer(array_to_permute,temp_array,perm,nspec)
-
- implicit none
-
- include "constants.h"
-
- integer, intent(in) :: nspec
- integer, intent(in), dimension(nspec) :: perm
-
- integer, intent(inout), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: array_to_permute,temp_array
-
- integer old_ispec,new_ispec
-
-! copy the original array
- temp_array(:,:,:,:) = array_to_permute(:,:,:,:)
-
- do old_ispec = 1,nspec
- new_ispec = perm(old_ispec)
- array_to_permute(:,:,:,new_ispec) = temp_array(:,:,:,old_ispec)
- enddo
-
- end subroutine permute_elements_integer
-
-!
-!-----------------------------------------------------------------------
-!
-
-! implement permutation of elements for arrays of double precision type
- subroutine permute_elements_dble(array_to_permute,temp_array,perm,nspec)
-
- implicit none
-
- include "constants.h"
-
- integer, intent(in) :: nspec
- integer, intent(in), dimension(nspec) :: perm
-
- double precision, intent(inout), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: array_to_permute,temp_array
-
- integer old_ispec,new_ispec
-
-! copy the original array
- temp_array(:,:,:,:) = array_to_permute(:,:,:,:)
-
- do old_ispec = 1,nspec
- new_ispec = perm(old_ispec)
- array_to_permute(:,:,:,new_ispec) = temp_array(:,:,:,old_ispec)
- enddo
-
- end subroutine permute_elements_dble
-
More information about the cig-commits
mailing list