SUBROUTINE METIS() USE PRE_GLOBAL C---------------------------------------------------------------------- C INTERFACE ROUTINE FOR PADCIRC TO USE THE METIS 4.0 LIBRARY C A GRAPH PARTITION TOOL, FOR DOMAIN DECOMPOSITION C Version 2.1 vjp 6/7/2006 C Redimensioned CO_NODES and IDUALS and added bounds test C Added check that adjacency matrix is symmetric C---------------------------------------------------------------------- C LOGICAL FOUND, SYMMETRIC INTEGER MNED, MNEDLOC, IDUMP INTEGER I, J, K, IEL, INODE, JNODE, NCOUNT, ITOT, NEDGETOT INTEGER MAXDUALS INTEGER WEIGHTFLAG, NUMFLAG, NPARTS, OPTIONS(5) INTEGER EDGECUT, OPTYPE, NBYTES C INTEGER,ALLOCATABLE :: IDUALS(:,:),ITVECT(:),ITVECT2(:) INTEGER,ALLOCATABLE :: XADJ(:), ADJNCY(:) INTEGER,ALLOCATABLE :: VWGTS(:), EWGTS(:) INTEGER,ALLOCATABLE :: CO_NODES(:,:),NEDGES(:), NEDLOC(:) INTEGER,ALLOCATABLE :: NUMDUALS(:) EXTERNAL metis_partgraphkway,metis_estimatememory C ALLOCATE ( ITVECT(MNP),ITVECT2(MNP) ) ALLOCATE ( NUMDUALS(MNP) ) ALLOCATE ( XADJ(MNP+1), VWGTS(MNP), NEDGES(MNP) ) ALLOCATE ( NEDLOC(MNP) ) C C--COMPUTE INDEX OF WEIR DUALS WHICH IS ZERO IF NOT A WEIR NODE C DO INODE=1, MNP NUMDUALS(INODE) = 0 ENDDO DO J=1, NWEIR NUMDUALS(WEIR(J)) = NUMDUALS(WEIR(J))+1 NUMDUALS(WEIRD(J)) = NUMDUALS(WEIRD(J))+1 ENDDO MAXDUALS = 0 DO J=1, MNP IF (NUMDUALS(J) .ge. MAXDUALS) MAXDUALS = NUMDUALS(J) ENDDO print *, "This model has ",NWEIR," weir node pairs" print *, "maximum number duals for any weir node = ", maxduals ALLOCATE ( IDUALS(MAXDUALS,MNP) ) DO INODE=1, MNP DO K=1, MAXDUALS IDUALS(K,INODE) = 0 ENDDO ENDDO DO J=1, NWEIR DO K=1, MAXDUALS IF (IDUALS(K,WEIR(J)) == 0) THEN IDUALS(K,WEIR(J)) = WEIRD(J) EXIT ENDIF ENDDO ENDDO DO J=1, NWEIR DO K=1, MAXDUALS IF (IDUALS(K,WEIRD(J)) == 0) THEN IDUALS(K,WEIRD(J)) = WEIR(J) EXIT ENDIF ENDDO ENDDO C------------------------------------------------------------- C COMPUTES THE TOTAL NUMBER OF EDGES --> MNED C AND THE MAX NUMBER OF EDGES FOR ANY NODE --> MNEDLOC C BOTH COUNTS INCLUDE WEIR-NODE PAIRS C------------------------------------------------------------- MNED = 0 DO INODE = 1,MNP NEDLOC(INODE) = 0 ENDDO DO J=1, 3 DO IEL=1, MNE INODE = NNEG(J,IEL) NCOUNT = NEDLOC(INODE) + 2 MNED = MNED + 2 DO K=1, MAXDUALS IF (IDUALS(K,INODE).NE.0) THEN NCOUNT = NCOUNT + 1 MNED = MNED + 1 ENDIF ENDDO NEDLOC(INODE) = NCOUNT ENDDO ENDDO MNEDLOC = 0 DO INODE=1, MNP IF (NEDLOC(INODE) .ge. MNEDLOC) MNEDLOC = NEDLOC(INODE) ENDDO c print *, "total number of edges = ", MNED print *, "maximum co-nodes for any node = ", MNEDLOC ALLOCATE ( ADJNCY(MNED), EWGTS(MNED) ) ALLOCATE ( CO_NODES(MNEDLOC,MNP) ) C C--COMPUTE CO_NODES LISTS AND NUMBER OF EDGES CONTAINING A NODE C DO INODE = 1,MNP NEDGES(INODE) = 0 ENDDO C DO IEL=1, MNE INODE = NNEG(1,IEL) CO_NODES(NEDGES(INODE)+1,INODE) = NNEG(2,IEL) CO_NODES(NEDGES(INODE)+2,INODE) = NNEG(3,IEL) NCOUNT = NEDGES(INODE) + 2 DO K=1, MAXDUALS IF (IDUALS(K,INODE).NE.0) THEN NCOUNT = NCOUNT + 1 CO_NODES(NCOUNT,INODE) = IDUALS(K,INODE) ENDIF ENDDO NEDGES(INODE) = NCOUNT ENDDO C DO IEL=1, MNE INODE = NNEG(2,IEL) CO_NODES(NEDGES(INODE)+1,INODE) = NNEG(3,IEL) CO_NODES(NEDGES(INODE)+2,INODE) = NNEG(1,IEL) NCOUNT = NEDGES(INODE) + 2 DO K=1, MAXDUALS IF (IDUALS(K,INODE).NE.0) THEN NCOUNT = NCOUNT + 1 CO_NODES(NCOUNT,INODE) = IDUALS(K,INODE) ENDIF ENDDO NEDGES(INODE) = NCOUNT ENDDO C DO IEL=1, MNE INODE = NNEG(3,IEL) CO_NODES(NEDGES(INODE)+1,INODE) = NNEG(1,IEL) CO_NODES(NEDGES(INODE)+2,INODE) = NNEG(2,IEL) NCOUNT = NEDGES(INODE) + 2 DO K=1, MAXDUALS IF (IDUALS(K,INODE).NE.0) THEN NCOUNT = NCOUNT + 1 CO_NODES(NCOUNT,INODE) = IDUALS(K,INODE) ENDIF ENDDO NEDGES(INODE) = NCOUNT ENDDO C C REMOVE REDUNDANCY IN NODE LISTS C NEDGETOT = 0 ! This will be twice number of edges DO INODE = 1,MNP DO J=1, NEDGES(INODE) ITVECT(J) = CO_NODES(J,INODE) ENDDO IF (NEDGES(INODE).GT.1) THEN NCOUNT = NEDGES(INODE) CALL SORT(NCOUNT,ITVECT) JNODE = ITVECT(1) CO_NODES(1,INODE) = JNODE NCOUNT = 1 DO J=2, NEDGES(INODE) IF (ITVECT(J).NE.JNODE) THEN NCOUNT = NCOUNT + 1 JNODE = ITVECT(J) CO_NODES(NCOUNT,INODE) = JNODE ENDIF ENDDO ELSE print *, "node = ",INODE," is isolated" stop 'vic' ENDIF NEDGES(INODE) = NCOUNT NEDGETOT = NEDGETOT + NCOUNT if (nedges(inode) == 0) then print *, "inode = ", inode, " belongs to no edges" stop 'vic' endif ENDDO NEDGETOT = NEDGETOT/2 print *, "edge count = ",nedgetot C check that adjacency matrix is symmetric C SYMMETRIC = .true. DO INODE = 1, MNP DO J = 1, NEDGES(INODE) JNODE = CO_NODES(J,INODE) FOUND = .false. DO K= 1, NEDGES(JNODE) IF (CO_NODES(K,JNODE) == INODE) THEN FOUND = .true. EXIT ENDIF ENDDO IF (.not. FOUND) THEN SYMMETRIC = .false. print *, "node ",inode," adjacent to ",jnode," but not visa-versa" ENDIF ENDDO ENDDO IF (.not. SYMMETRIC) THEN stop 'bad adjacency matrix: not symmetric!' ENDIF C C COMPUTE WEIGHTS OF THE GRAPH VERTICES C DO INODE = 1,MNP VWGTS(INODE) = NEDGES(INODE) ENDDO C C--COMPUTE ADJACENCY LIST OF GRAPH AND ITS EDGE WEIGHTS C XADJ(1) = 1 ITOT = 0 DO INODE = 1,MNP DO J = 1, NEDGES(INODE) ITOT = ITOT + 1 JNODE = CO_NODES(J,INODE) ADJNCY(ITOT) = JNODE EWGTS(ITOT) = (VWGTS(JNODE)+VWGTS(INODE)) ENDDO XADJ(INODE+1) = ITOT+1 ENDDO C C Dump graph to a file for debugging C IDUMP = 1 IF (IDUMP.EQ.1) THEN OPEN(FILE='metis_graph.txt',UNIT=99) WRITE(99,100) MNP, NEDGETOT, 11, 1 DO INODE=1, MNP WRITE(99,200) VWGTS(INODE), & (CO_NODES(J,INODE), EWGTS(XADJ(INODE)+J-1),J=1,NEDGES(INODE)) ENDDO CLOSE(99) ENDIF C C--CALL K-WAY METIS FOR PARTITIONING C NUMFLAG = 1 NPARTS = MNPROC OPTIONS(1) = 1 OPTIONS(2) = 3 OPTIONS(3) = 1 OPTIONS(4) = 3 ! minimize number of co-domains OPTIONS(5) = 0 c WEIGHTFLAG = 3 ! use weights for nodes and edges C OPTYPE = 2 CALL metis_estimatememory(MNP,XADJ,ADJNCY,NUMFLAG, & OPTYPE,NBYTES) print *, "" print *, "Grid Partition Data" print *, "METIS 4.0 will require approximately ",nbytes," bytes" C CALL metis_partgraphkway( MNP,XADJ,ADJNCY,VWGTS,EWGTS, & WEIGHTFLAG,NUMFLAG,NPARTS,OPTIONS,EDGECUT,PROC) C print *, "Total Edges Cut = ",EDGECUT print *, "writing mesh partition to file: partmesh.txt" OPEN(990,FILE='partmesh.txt') DO I=1, MNP WRITE(990,*) PROC(I) ENDDO CLOSE(990) 100 FORMAT(4I10) 200 FORMAT(100I10) C RETURN END