1!------------------------------------------------------------------------------- 2! Copyright (c) 2020 FrontISTR Commons 3! This software is released under the MIT License, see LICENSE.txt 4!------------------------------------------------------------------------------- 5 6!---------------------------------------------------------------------- 7!> @brief HECMW_ORDERING_QMD is a program for the minimum degree 8! ordering using quotient graphs 9!---------------------------------------------------------------------- 10module hecmw_ordering_qmd 11 use hecmw_util 12 implicit none 13 14 private 15 public :: hecmw_ordering_genqmd 16 17contains 18 19 !======================================================================! 20 !> @brief hecmw_ordering_GENQMD 21 !======================================================================! 22 subroutine hecmw_ordering_GENQMD(Neqns,Nttbr,Xadj,Adj0,Perm,Invp) 23 implicit none 24 !------ 25 integer(kind=kint), intent(in):: Neqns 26 integer(kind=kint), intent(in):: Nttbr 27 integer(kind=kint), intent(in):: Adj0(:) 28 integer(kind=kint), intent(in):: Xadj(:) 29 integer(kind=kint), intent(out):: Perm(:) 30 integer(kind=kint), intent(out):: Invp(:) 31 !------ 32 integer(kind=kint), allocatable:: Deg(:) 33 integer(kind=kint), allocatable:: Marker(:) 34 integer(kind=kint), allocatable:: Rchset(:) 35 integer(kind=kint), allocatable:: Nbrhd(:) 36 integer(kind=kint), allocatable:: Qsize(:) 37 integer(kind=kint), allocatable:: Qlink(:) 38 integer(kind=kint), allocatable:: Adjncy(:) 39 integer(kind=kint):: inode 40 integer(kind=kint):: ip 41 integer(kind=kint):: irch 42 integer(kind=kint):: j 43 integer(kind=kint):: mindeg 44 integer(kind=kint):: ndeg 45 integer(kind=kint):: nhdsze 46 integer(kind=kint):: node 47 integer(kind=kint):: np 48 integer(kind=kint):: num 49 integer(kind=kint):: nump1 50 integer(kind=kint):: nxnode 51 integer(kind=kint):: rchsze 52 integer(kind=kint):: search 53 integer(kind=kint):: thresh 54 integer(kind=kint):: ierror 55 logical:: found 56 57 allocate (DEG(NEQns+1),stat=IERror) 58 if ( IERror/=0 ) stop "ALLOCATION ERROR, deg: SUB. genqmd" 59 allocate (MARker(NEQns+1),stat=IERror) 60 if ( IERror/=0 ) stop "ALLOCATION ERROR, marker: SUB. genqmd" 61 allocate (RCHset(NEQns+2),stat=IERror) 62 if ( IERror/=0 ) stop "ALLOCATION ERROR, rchset: SUB. genqmd" 63 allocate (NBRhd(NEQns+1),stat=IERror) 64 if ( IERror/=0 ) stop "ALLOCATION ERROR, nbrhd: SUB. genqmd" 65 allocate (QSIze(NEQns+1),stat=IERror) 66 if ( IERror/=0 ) stop "ALLOCATION ERROR, qsize: SUB. genqmd" 67 allocate (QLInk(NEQns+1),stat=IERror) 68 if ( IERror/=0 ) stop "ALLOCATION ERROR, qlink: SUB. genqmd" 69 allocate (ADJncy(2*NTTbr),stat=IERror) 70 if ( IERror/=0 ) stop "ALLOCATION ERROR, adjncy: SUB. genqmd" 71 72 mindeg = Neqns 73 Adjncy(1:Xadj(Neqns+1) - 1) = Adj0(1:Xadj(Neqns+1) - 1) 74 do node = 1, Neqns 75 Perm(node) = node 76 Invp(node) = node 77 Marker(node) = 0 78 Qsize(node) = 1 79 Qlink(node) = 0 80 ndeg = Xadj(node+1) - Xadj(node) 81 Deg(node) = ndeg 82 if ( ndeg<mindeg ) mindeg = ndeg 83 enddo 84 85 num = 0 86 loop1: do 87 search = 1 88 thresh = mindeg 89 mindeg = Neqns 90 loop2: do 91 nump1 = num + 1 92 if ( nump1>search ) search = nump1 93 found = .false. 94 do j = search, Neqns 95 node = Perm(j) 96 if ( Marker(node)>=0 ) then 97 ndeg = Deg(node) 98 if ( ndeg<=thresh ) then 99 found = .true. 100 exit 101 endif 102 if ( ndeg<mindeg ) mindeg = ndeg 103 endif 104 enddo 105 if (.not. found) cycle loop1 106 107 search = j 108 Marker(node) = 1 109 call QMDRCH(node,Xadj,Adjncy,Deg,Marker,rchsze,Rchset,nhdsze,Nbrhd) 110 nxnode = node 111 do 112 num = num + 1 113 np = Invp(nxnode) 114 ip = Perm(num) 115 Perm(np) = ip 116 Invp(ip) = np 117 Perm(num) = nxnode 118 Invp(nxnode) = num 119 Deg(nxnode) = -1 120 nxnode = Qlink(nxnode) 121 if ( nxnode<=0 ) then 122 if ( rchsze>0 ) then 123 call QMDUPD(Xadj,Adjncy,rchsze,Rchset,Deg,Qsize,Qlink,Marker,Rchset(rchsze+1:),Nbrhd(nhdsze+1:)) 124 Marker(node) = 0 125 do irch = 1, rchsze 126 inode = Rchset(irch) 127 if ( Marker(inode)>=0 ) then 128 Marker(inode) = 0 129 ndeg = Deg(inode) 130 if ( ndeg<mindeg ) mindeg = ndeg 131 if ( ndeg<=thresh ) then 132 mindeg = thresh 133 thresh = ndeg 134 search = Invp(inode) 135 endif 136 endif 137 enddo 138 if ( nhdsze>0 ) call QMDOT(node,Xadj,Adjncy,Marker,rchsze,Rchset,Nbrhd) 139 endif 140 if ( num>=Neqns ) exit 141 cycle loop2 142 endif 143 enddo 144 exit 145 enddo loop2 146 exit 147 enddo loop1 148 149 deallocate (DEG) 150 deallocate (MARker) 151 deallocate (RCHset) 152 deallocate (NBRhd) 153 deallocate (QSIze) 154 deallocate (QLInk) 155 deallocate (ADJncy) 156 end subroutine HECMW_ORDERING_GENQMD 157 158 !======================================================================! 159 !> @brief QMDRCH 160 !======================================================================! 161 subroutine QMDRCH(Root,Xadj,Adjncy,Deg,Marker,Rchsze,Rchset,Nhdsze,Nbrhd) 162 implicit none 163 !------ 164 integer(kind=kint), intent(in):: Root 165 integer(kind=kint), intent(in):: Adjncy(:) 166 integer(kind=kint), intent(in):: Deg(:) 167 integer(kind=kint), intent(in):: Xadj(:) 168 integer(kind=kint), intent(out):: Nhdsze 169 integer(kind=kint), intent(out):: Rchsze 170 integer(kind=kint), intent(inout):: Marker(:) 171 integer(kind=kint), intent(out):: Rchset(:) 172 integer(kind=kint), intent(out):: Nbrhd(:) 173 !------ 174 integer(kind=kint):: i 175 integer(kind=kint):: istrt 176 integer(kind=kint):: istop 177 integer(kind=kint):: j 178 integer(kind=kint):: jstrt 179 integer(kind=kint):: jstop 180 integer(kind=kint):: nabor 181 integer(kind=kint):: node 182 183 Nhdsze = 0 184 Rchsze = 0 185 istrt = Xadj(Root) 186 istop = Xadj(Root+1) - 1 187 if ( istop<istrt ) return 188 do i = istrt, istop 189 nabor = Adjncy(i) 190 if ( nabor==0 ) return 191 if ( Marker(nabor)==0 ) then 192 if ( Deg(nabor)<0 ) then 193 Marker(nabor) = -1 194 Nhdsze = Nhdsze + 1 195 Nbrhd(Nhdsze) = nabor 196 loop1: do 197 jstrt = Xadj(nabor) 198 jstop = Xadj(nabor+1) - 1 199 do j = jstrt, jstop 200 node = Adjncy(j) 201 nabor = -node 202 if ( node<0 ) cycle loop1 203 if ( node==0 ) exit 204 if ( Marker(node)==0 ) then 205 Rchsze = Rchsze + 1 206 Rchset(Rchsze) = node 207 Marker(node) = 1 208 endif 209 enddo 210 exit 211 enddo loop1 212 else 213 Rchsze = Rchsze + 1 214 Rchset(Rchsze) = nabor 215 Marker(nabor) = 1 216 endif 217 endif 218 enddo 219 end subroutine QMDRCH 220 221 !======================================================================! 222 !> @brief QMDUPD 223 !======================================================================! 224 subroutine QMDUPD(Xadj,Adjncy,Nlist,List,Deg,Qsize,Qlink,Marker,Rchset,Nbrhd) 225 implicit none 226 !------ 227 integer(kind=kint), intent(in):: Nlist 228 integer(kind=kint), intent(in):: Adjncy(:) 229 integer(kind=kint), intent(in):: List(:) 230 integer(kind=kint), intent(in):: Xadj(:) 231 integer(kind=kint), intent(inout):: Deg(:) 232 integer(kind=kint), intent(inout):: Marker(:) 233 integer(kind=kint), intent(out):: Rchset(:) 234 integer(kind=kint), intent(out):: Nbrhd(:) 235 integer(kind=kint), intent(inout):: Qsize(:) 236 integer(kind=kint), intent(inout):: Qlink(:) 237 !------ 238 integer(kind=kint):: deg0 239 integer(kind=kint):: deg1 240 integer(kind=kint):: il 241 integer(kind=kint):: inhd 242 integer(kind=kint):: inode 243 integer(kind=kint):: irch 244 integer(kind=kint):: j 245 integer(kind=kint):: jstrt 246 integer(kind=kint):: jstop 247 integer(kind=kint):: mark 248 integer(kind=kint):: nabor 249 integer(kind=kint):: nhdsze 250 integer(kind=kint):: node 251 integer(kind=kint):: rchsze 252 253 if ( Nlist<=0 ) return 254 deg0 = 0 255 nhdsze = 0 256 do il = 1, Nlist 257 node = List(il) 258 deg0 = deg0 + Qsize(node) 259 jstrt = Xadj(node) 260 jstop = Xadj(node+1) - 1 261 do j = jstrt, jstop 262 nabor = Adjncy(j) 263 if ( Marker(nabor)==0 .and. Deg(nabor)<0 ) then 264 Marker(nabor) = -1 265 nhdsze = nhdsze + 1 266 Nbrhd(nhdsze) = nabor 267 endif 268 enddo 269 enddo 270 271 if ( nhdsze>0 ) call QMDMRG(Xadj,Adjncy,Deg,Qsize,Qlink,Marker,deg0,nhdsze,Nbrhd,Rchset,Nbrhd(nhdsze+1:)) 272 do il = 1, Nlist 273 node = List(il) 274 mark = Marker(node) 275 if ( mark<=1 .and. mark>=0 ) then 276 call QMDRCH(node,Xadj,Adjncy,Deg,Marker,rchsze,Rchset,nhdsze,Nbrhd) 277 deg1 = deg0 278 if ( rchsze>0 ) then 279 do irch = 1, rchsze 280 inode = Rchset(irch) 281 deg1 = deg1 + Qsize(inode) 282 Marker(inode) = 0 283 enddo 284 endif 285 Deg(node) = deg1 - 1 286 if ( nhdsze>0 ) then 287 do inhd = 1, nhdsze 288 inode = Nbrhd(inhd) 289 Marker(inode) = 0 290 enddo 291 endif 292 endif 293 enddo 294 end subroutine QMDUPD 295 296 !======================================================================! 297 !> @brief QMDMRG 298 !======================================================================! 299 subroutine QMDMRG(Xadj,Adjncy,Deg,Qsize,Qlink,Marker,Deg0,Nhdsze,Nbrhd,Rchset,Ovrlp) 300 implicit none 301 !------ 302 integer(kind=kint), intent(in):: Deg0 303 integer(kind=kint), intent(in):: Nhdsze 304 integer(kind=kint), intent(in):: Adjncy(:) 305 integer(kind=kint), intent(in):: Nbrhd(:) 306 integer(kind=kint), intent(in):: Xadj(:) 307 integer(kind=kint), intent(inout):: Deg(:) 308 integer(kind=kint), intent(inout):: Qsize(:) 309 integer(kind=kint), intent(inout):: Qlink(:) 310 integer(kind=kint), intent(inout):: Marker(:) 311 integer(kind=kint), intent(out):: Rchset(:) 312 integer(kind=kint), intent(out):: Ovrlp(:) 313 !------ 314 integer(kind=kint):: deg1 315 integer(kind=kint):: head 316 integer(kind=kint):: inhd 317 integer(kind=kint):: iov 318 integer(kind=kint):: irch 319 integer(kind=kint):: j 320 integer(kind=kint):: jstrt 321 integer(kind=kint):: jstop 322 integer(kind=kint):: link 323 integer(kind=kint):: lnode 324 integer(kind=kint):: mark 325 integer(kind=kint):: mrgsze 326 integer(kind=kint):: nabor 327 integer(kind=kint):: node 328 integer(kind=kint):: novrlp 329 integer(kind=kint):: rchsze 330 integer(kind=kint):: root 331 332 if ( Nhdsze<=0 ) return 333 do inhd = 1, Nhdsze 334 root = Nbrhd(inhd) 335 Marker(root) = 0 336 enddo 337 do inhd = 1, Nhdsze 338 root = Nbrhd(inhd) 339 Marker(root) = -1 340 rchsze = 0 341 novrlp = 0 342 deg1 = 0 343 loop1: do 344 jstrt = Xadj(root) 345 jstop = Xadj(root+1) - 1 346 do j = jstrt, jstop 347 nabor = Adjncy(j) 348 root = -nabor 349 if ( nabor<0 ) cycle loop1 350 if ( nabor==0 ) exit 351 mark = Marker(nabor) 352 353 if ( mark>=0 ) then 354 if ( mark<=0 ) then 355 rchsze = rchsze + 1 356 Rchset(rchsze) = nabor 357 deg1 = deg1 + Qsize(nabor) 358 Marker(nabor) = 1 359 elseif ( mark<=1 ) then 360 novrlp = novrlp + 1 361 Ovrlp(novrlp) = nabor 362 Marker(nabor) = 2 363 endif 364 endif 365 enddo 366 exit 367 enddo loop1 368 head = 0 369 mrgsze = 0 370 loop2: do iov = 1, novrlp 371 node = Ovrlp(iov) 372 jstrt = Xadj(node) 373 jstop = Xadj(node+1) - 1 374 do j = jstrt, jstop 375 nabor = Adjncy(j) 376 if ( Marker(nabor)==0 ) then 377 Marker(node) = 1 378 cycle loop2 379 endif 380 enddo 381 mrgsze = mrgsze + Qsize(node) 382 Marker(node) = -1 383 lnode = node 384 do 385 link = Qlink(lnode) 386 if ( link<=0 ) then 387 Qlink(lnode) = head 388 head = node 389 exit 390 else 391 lnode = link 392 endif 393 enddo 394 enddo loop2 395 if ( head>0 ) then 396 Qsize(head) = mrgsze 397 Deg(head) = Deg0 + deg1 - 1 398 Marker(head) = 2 399 endif 400 root = Nbrhd(inhd) 401 Marker(root) = 0 402 if ( rchsze>0 ) then 403 do irch = 1, rchsze 404 node = Rchset(irch) 405 Marker(node) = 0 406 enddo 407 endif 408 enddo 409 end subroutine QMDMRG 410 411 !======================================================================! 412 !> @brief QMDOT 413 !======================================================================! 414 subroutine QMDOT(Root,Xadj,Adjncy,Marker,Rchsze,Rchset,Nbrhd) 415 implicit none 416 !------ 417 integer(kind=kint), intent(in):: Rchsze 418 integer(kind=kint), intent(in):: Root 419 integer(kind=kint), intent(in):: Marker(:) 420 integer(kind=kint), intent(in):: Rchset(:) 421 integer(kind=kint), intent(in):: Nbrhd(:) 422 integer(kind=kint), intent(in):: Xadj(:) 423 integer(kind=kint), intent(inout):: Adjncy(:) 424 !------ 425 integer(kind=kint):: inhd 426 integer(kind=kint):: irch 427 integer(kind=kint):: j 428 integer(kind=kint):: jstrt 429 integer(kind=kint):: jstop 430 integer(kind=kint):: link 431 integer(kind=kint):: nabor 432 integer(kind=kint):: node 433 434 irch = 0 435 inhd = 0 436 node = Root 437 loop1: do 438 jstrt = Xadj(node) 439 jstop = Xadj(node+1) - 2 440 if ( jstop>=jstrt ) then 441 do j = jstrt, jstop 442 irch = irch + 1 443 Adjncy(j) = Rchset(irch) 444 if ( irch>=Rchsze ) exit loop1 445 enddo 446 endif 447 link = Adjncy(jstop+1) 448 node = -link 449 if ( link>=0 ) then 450 inhd = inhd + 1 451 node = Nbrhd(inhd) 452 Adjncy(jstop+1) = -node 453 endif 454 enddo loop1 455 456 Adjncy(j+1) = 0 457 do irch = 1, Rchsze 458 node = Rchset(irch) 459 if ( Marker(node)>=0 ) then 460 jstrt = Xadj(node) 461 jstop = Xadj(node+1) - 1 462 do j = jstrt, jstop 463 nabor = Adjncy(j) 464 if ( Marker(nabor)<0 ) then 465 Adjncy(j) = Root 466 exit 467 endif 468 enddo 469 endif 470 enddo 471 end subroutine QMDOT 472 473end module hecmw_ordering_qmd 474