00001
00002
00003
00004 module ice_spacecurve
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 use ice_kinds_mod
00017
00018 implicit none
00019
00020
00021
00022 type, public :: factor_t
00023 integer(int_kind) :: numfact
00024 integer(int_kind), dimension(:),pointer :: factors
00025 integer(int_kind), dimension(:), pointer :: used
00026 end type
00027
00028
00029
00030 public :: GenSpaceCurve, &
00031 IsLoadBalanced
00032
00033 public :: Factor, &
00034 IsFactorable, &
00035 PrintFactor, &
00036 ProdFactor, &
00037 MatchFactor
00038
00039
00040
00041 private :: map, &
00042 PeanoM, &
00043 Hilbert, &
00044 Cinco, &
00045 GenCurve
00046
00047 private :: FirstFactor, &
00048 FindandMark
00049
00050 integer(int_kind), dimension(:,:), allocatable ::
00051 dir,
00052 ordered
00053 integer(int_kind), dimension(:), allocatable ::
00054 pos
00055
00056 integer(int_kind) ::
00057 maxdim,
00058 vcnt
00059
00060 logical :: verbose=.FALSE.
00061
00062 type (factor_t), public,save :: fact
00063
00064
00065
00066
00067
00068 contains
00069
00070
00071
00072
00073
00074
00075 recursive function Cinco(l,type,ma,md,ja,jd) result(ierr)
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090 integer(int_kind), intent(in) ::
00091 l,
00092 type,
00093 ma,
00094 md,
00095 ja,
00096 jd
00097
00098
00099
00100 integer(int_kind) :: ierr
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110 integer(int_kind) ::
00111 lma,
00112 lmd,
00113 lja,
00114 ljd,
00115 ltype,
00116 ll
00117
00118 logical :: debug = .FALSE.
00119
00120
00121 ll = l
00122 if(ll .gt. 1) ltype = fact%factors(ll-1)
00123
00124
00125
00126
00127 lma = ma
00128 lmd = md
00129 lja = lma
00130 ljd = lmd
00131
00132 if(ll .gt. 1) then
00133 if(debug) write(*,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00134 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00135 if(debug) call PrintCurve(ordered)
00136 else
00137 ierr = IncrementCurve(lja,ljd)
00138 if(debug) print *,'Cinco: After Position [0,0] ',pos
00139 endif
00140
00141
00142
00143
00144 lma = ma
00145 lmd = md
00146 lja = lma
00147 ljd = lmd
00148
00149 if(ll .gt. 1) then
00150 if(debug) write(*,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00151 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00152 if(debug) call PrintCurve(ordered)
00153 else
00154 ierr = IncrementCurve(lja,ljd)
00155 if(debug) print *,'After Position [1,0] ',pos
00156 endif
00157
00158
00159
00160
00161 lma = MOD(ma+1,maxdim)
00162 lmd = md
00163 lja = lma
00164 ljd = lmd
00165
00166 if(ll .gt. 1) then
00167 if(debug) write(*,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00168 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00169 if(debug) call PrintCurve(ordered)
00170 else
00171 ierr = IncrementCurve(lja,ljd)
00172 if(debug) print *,'After Position [2,0] ',pos
00173 endif
00174
00175
00176
00177
00178 lma = MOD(ma+1,maxdim)
00179 lmd = md
00180 lja = lma
00181 ljd = lmd
00182
00183 if(ll .gt. 1) then
00184 if(debug) write(*,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00185 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00186 if(debug) call PrintCurve(ordered)
00187 else
00188 ierr = IncrementCurve(lja,ljd)
00189 if(debug) print *,'After Position [2,1] ',pos
00190 endif
00191
00192
00193
00194
00195 lma = MOD(ma+1,maxdim)
00196 lmd = md
00197 lja = ma
00198 ljd = -md
00199
00200 if(ll .gt. 1) then
00201 if(debug) write(*,25) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00202 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00203 if(debug) call PrintCurve(ordered)
00204 else
00205 ierr = IncrementCurve(lja,ljd)
00206 if(debug) print *,'After Position [2,2] ',pos
00207 endif
00208
00209
00210
00211
00212 lma = MOD(ma+1,maxdim)
00213 lmd = -md
00214 lja = lma
00215 ljd = lmd
00216
00217 if(ll .gt. 1) then
00218 if(debug) write(*,26) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00219 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00220 if(debug) call PrintCurve(ordered)
00221 else
00222 ierr = IncrementCurve(lja,ljd)
00223 if(debug) print *,'After Position [1,2] ',pos
00224 endif
00225
00226
00227
00228
00229 lma = ma
00230 lmd = -md
00231 lja = lma
00232 ljd = lmd
00233
00234 if(ll .gt. 1) then
00235 if(debug) write(*,27) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00236 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00237 if(debug) call PrintCurve(ordered)
00238 else
00239 ierr = IncrementCurve(lja,ljd)
00240 if(debug) print *,'After Position [1,1] ',pos
00241 endif
00242
00243
00244
00245
00246 lma = ma
00247 lmd = -md
00248 lja = MOD(ma+1,maxdim)
00249 ljd = md
00250
00251 if(ll .gt. 1) then
00252 if(debug) write(*,28) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00253 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00254 if(debug) call PrintCurve(ordered)
00255 else
00256 ierr = IncrementCurve(lja,ljd)
00257 if(debug) print *,'After Position [0,1] ',pos
00258 endif
00259
00260
00261
00262
00263 lma = MOD(ma+1,maxdim)
00264 lmd = md
00265 lja = lma
00266 ljd = lmd
00267
00268 if(ll .gt. 1) then
00269 if(debug) write(*,29) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00270 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00271 if(debug) call PrintCurve(ordered)
00272 else
00273 ierr = IncrementCurve(lja,ljd)
00274 if(debug) print *,'After Position [0,2] ',pos
00275 endif
00276
00277
00278
00279
00280 lma = MOD(ma+1,maxdim)
00281 lmd = md
00282 lja = lma
00283 ljd = lmd
00284
00285 if(ll .gt. 1) then
00286 if(debug) write(*,30) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00287 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00288 if(debug) call PrintCurve(ordered)
00289 else
00290 ierr = IncrementCurve(lja,ljd)
00291 if(debug) print *,'After Position [0,3] ',pos
00292 endif
00293
00294
00295
00296
00297 lma = ma
00298 lmd = md
00299 lja = lma
00300 ljd = lmd
00301
00302 if(ll .gt. 1) then
00303 if(debug) write(*,31) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00304 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00305 if(debug) call PrintCurve(ordered)
00306 else
00307 ierr = IncrementCurve(lja,ljd)
00308 if(debug) print *,'After Position [0,4] ',pos
00309 endif
00310
00311
00312
00313
00314 lma = ma
00315 lmd = md
00316 lja = MOD(ma+1,maxdim)
00317 ljd = -md
00318
00319 if(ll .gt. 1) then
00320 if(debug) write(*,32) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00321 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00322 if(debug) call PrintCurve(ordered)
00323 else
00324 ierr = IncrementCurve(lja,ljd)
00325 if(debug) print *,'After Position [1,4] ',pos
00326 endif
00327
00328
00329
00330
00331 lma = MOD(ma+1,maxdim)
00332 lmd = -md
00333 lja = ma
00334 ljd = md
00335
00336 if(ll .gt. 1) then
00337 if(debug) write(*,33) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00338 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00339 if(debug) call PrintCurve(ordered)
00340 else
00341 ierr = IncrementCurve(lja,ljd)
00342 if(debug) print *,'After Position [1,3] ',pos
00343 endif
00344
00345
00346
00347
00348 lma = MOD(ma+1,maxdim)
00349 lmd = md
00350 lja = lma
00351 ljd = lmd
00352
00353 if(ll .gt. 1) then
00354 if(debug) write(*,34) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00355 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00356 if(debug) call PrintCurve(ordered)
00357 else
00358 ierr = IncrementCurve(lja,ljd)
00359 if(debug) print *,'After Position [2,3] ',pos
00360 endif
00361
00362
00363
00364
00365 lma = ma
00366 lmd = md
00367 lja = lma
00368 ljd = lmd
00369
00370 if(ll .gt. 1) then
00371 if(debug) write(*,35) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00372 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00373 if(debug) call PrintCurve(ordered)
00374 else
00375 ierr = IncrementCurve(lja,ljd)
00376 if(debug) print *,'After Position [2,4] ',pos
00377 endif
00378
00379
00380
00381
00382 lma = ma
00383 lmd = md
00384 lja = lma
00385 ljd = lmd
00386
00387 if(ll .gt. 1) then
00388 if(debug) write(*,36) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00389 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00390 if(debug) call PrintCurve(ordered)
00391 else
00392 ierr = IncrementCurve(lja,ljd)
00393 if(debug) print *,'After Position [3,4] ',pos
00394 endif
00395
00396
00397
00398
00399 lma = ma
00400 lmd = md
00401 lja = MOD(ma+1,maxdim)
00402 ljd = -md
00403
00404 if(ll .gt. 1) then
00405 if(debug) write(*,37) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00406 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00407 if(debug) call PrintCurve(ordered)
00408 else
00409 ierr = IncrementCurve(lja,ljd)
00410 if(debug) print *,'After Position [4,4] ',pos
00411 endif
00412
00413
00414
00415
00416 lma = ma
00417 lmd = -md
00418 lja = lma
00419 ljd = lmd
00420
00421 if(ll .gt. 1) then
00422 if(debug) write(*,38) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00423 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00424 if(debug) call PrintCurve(ordered)
00425 else
00426 ierr = IncrementCurve(lja,ljd)
00427 if(debug) print *,'After Position [4,3] ',pos
00428 endif
00429
00430
00431
00432
00433 lma = MOD(ma+1,maxdim)
00434 lmd = -md
00435 lja = lma
00436 ljd = lmd
00437
00438 if(ll .gt. 1) then
00439 if(debug) write(*,39) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00440 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00441 if(debug) call PrintCurve(ordered)
00442 else
00443 ierr = IncrementCurve(lja,ljd)
00444 if(debug) print *,'After Position [3,3] ',pos
00445 endif
00446
00447
00448
00449
00450 lma = MOD(ma+1,maxdim)
00451 lmd = -md
00452 lja = ma
00453 ljd = md
00454
00455 if(ll .gt. 1) then
00456 if(debug) write(*,40) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00457 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00458 if(debug) call PrintCurve(ordered)
00459 else
00460 ierr = IncrementCurve(lja,ljd)
00461 if(debug) print *,'After Position [3,2] ',pos
00462 endif
00463
00464
00465
00466
00467 lma = ma
00468 lmd = md
00469 lja = MOD(ma+1,maxdim)
00470 ljd = -md
00471
00472 if(ll .gt. 1) then
00473 if(debug) write(*,41) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00474 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00475 if(debug) call PrintCurve(ordered)
00476 else
00477 ierr = IncrementCurve(lja,ljd)
00478 if(debug) print *,'After Position [4,2] ',pos
00479 endif
00480
00481
00482
00483
00484 lma = ma
00485 lmd = -md
00486 lja = lma
00487 ljd = lmd
00488
00489 if(ll .gt. 1) then
00490 if(debug) write(*,42) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00491 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00492 if(debug) call PrintCurve(ordered)
00493 else
00494 ierr = IncrementCurve(lja,ljd)
00495 if(debug) print *,'After Position [4,1] ',pos
00496 endif
00497
00498
00499
00500
00501 lma = MOD(ma+1,maxdim)
00502 lmd = -md
00503 lja = lma
00504 ljd = lmd
00505
00506 if(ll .gt. 1) then
00507 if(debug) write(*,43) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00508 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00509 if(debug) call PrintCurve(ordered)
00510 else
00511 ierr = IncrementCurve(lja,ljd)
00512 if(debug) print *,'After Position [3,1] ',pos
00513 endif
00514
00515
00516
00517
00518 lma = MOD(ma+1,maxdim)
00519 lmd = -md
00520 lja = ma
00521 ljd = md
00522
00523 if(ll .gt. 1) then
00524 if(debug) write(*,44) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00525 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00526 if(debug) call PrintCurve(ordered)
00527 else
00528 ierr = IncrementCurve(lja,ljd)
00529 if(debug) print *,'After Position [3,0] ',pos
00530 endif
00531
00532
00533
00534
00535 lma = ma
00536 lmd = md
00537 lja = ja
00538 ljd = jd
00539
00540 if(ll .gt. 1) then
00541 if(debug) write(*,45) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00542 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00543 if(debug) call PrintCurve(ordered)
00544 else
00545 ierr = IncrementCurve(lja,ljd)
00546 if(debug) print *,'After Position [4,0] ',pos
00547 endif
00548
00549 21 format('Call Cinco Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
00550 22 format('Call Cinco Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
00551 23 format('Call Cinco Pos [2,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
00552 24 format('Call Cinco Pos [2,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
00553 25 format('Call Cinco Pos [2,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
00554 26 format('Call Cinco Pos [1,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
00555 27 format('Call Cinco Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
00556 28 format('Call Cinco Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
00557 29 format('Call Cinco Pos [0,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
00558 30 format('Call Cinco Pos [0,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
00559 31 format('Call Cinco Pos [0,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
00560 32 format('Call Cinco Pos [1,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
00561 33 format('Call Cinco Pos [1,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
00562 34 format('Call Cinco Pos [2,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
00563 35 format('Call Cinco Pos [2,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
00564 36 format('Call Cinco Pos [3,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
00565 37 format('Call Cinco Pos [4,4] Level ',i1,' at (',i2,',',i2,')',4(i3))
00566 38 format('Call Cinco Pos [4,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
00567 39 format('Call Cinco Pos [3,3] Level ',i1,' at (',i2,',',i2,')',4(i3))
00568 40 format('Call Cinco Pos [3,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
00569 41 format('Call Cinco Pos [4,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
00570 42 format('Call Cinco Pos [4,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
00571 43 format('Call Cinco Pos [3,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
00572 44 format('Call Cinco Pos [3,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
00573 45 format('Call Cinco Pos [4,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
00574
00575
00576
00577
00578 end function Cinco
00579
00580
00581
00582
00583
00584
00585 recursive function PeanoM(l,type,ma,md,ja,jd) result(ierr)
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600 integer(int_kind), intent(in) ::
00601 l,
00602 type,
00603 ma,
00604 md,
00605 ja,
00606 jd
00607
00608
00609
00610 integer(int_kind) :: ierr
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621 integer(int_kind) ::
00622 lma,
00623 lmd,
00624 lja,
00625 ljd,
00626 ltype,
00627 ll
00628
00629 logical :: debug = .FALSE.
00630
00631
00632
00633 ll = l
00634 if(ll .gt. 1) ltype = fact%factors(ll-1)
00635
00636
00637
00638 lma = MOD(ma+1,maxdim)
00639 lmd = md
00640 lja = lma
00641 ljd = lmd
00642
00643 if(ll .gt. 1) then
00644 if(debug) write(*,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00645 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00646 if(debug) call PrintCurve(ordered)
00647 else
00648 ierr = IncrementCurve(lja,ljd)
00649 if(debug) print *,'PeanoM: After Position [0,0] ',pos
00650 endif
00651
00652
00653
00654
00655
00656 lma = MOD(ma+1,maxdim)
00657 lmd = md
00658 lja = lma
00659 ljd = lmd
00660 if(ll .gt. 1) then
00661 if(debug) write(*,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00662 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00663 if(debug) call PrintCurve(ordered)
00664 else
00665 ierr = IncrementCurve(lja,ljd)
00666 if(debug) print *,'PeanoM: After Position [0,1] ',pos
00667 endif
00668
00669
00670
00671
00672 lma = ma
00673 lmd = md
00674 lja = lma
00675 ljd = lmd
00676 if(ll .gt. 1) then
00677 if(debug) write(*,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00678 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00679 if(debug) call PrintCurve(ordered)
00680 else
00681 ierr = IncrementCurve(lja,ljd)
00682 if(debug) print *,'PeanoM: After Position [0,2] ',pos
00683 endif
00684
00685
00686
00687
00688 lma = ma
00689 lmd = md
00690 lja = lma
00691 ljd = lmd
00692 if(ll .gt. 1) then
00693 if(debug) write(*,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00694 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00695 if(debug) call PrintCurve(ordered)
00696 else
00697 ierr = IncrementCurve(lja,ljd)
00698 if(debug) print *,'PeanoM: After Position [1,2] ',pos
00699 endif
00700
00701
00702
00703
00704
00705 lma = ma
00706 lmd = md
00707 lja = MOD(lma+1,maxdim)
00708 ljd = -lmd
00709
00710 if(ll .gt. 1) then
00711 if(debug) write(*,25) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00712 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00713 if(debug) call PrintCurve(ordered)
00714 else
00715 ierr = IncrementCurve(lja,ljd)
00716 if(debug) print *,'PeanoM: After Position [2,2] ',pos
00717 endif
00718
00719
00720
00721
00722 lma = ma
00723 lmd = -md
00724 lja = lma
00725 ljd = lmd
00726
00727 if(ll .gt. 1) then
00728 if(debug) write(*,26) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00729 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00730 if(debug) call PrintCurve(ordered)
00731 else
00732 ierr = IncrementCurve(lja,ljd)
00733 if(debug) print *,'PeanoM: After Position [2,1] ',pos
00734 endif
00735
00736
00737
00738
00739 lma = MOD(ma+1,maxdim)
00740 lmd = -md
00741 lja = lma
00742 ljd = lmd
00743
00744 if(ll .gt. 1) then
00745 if(debug) write(*,27) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00746 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00747 if(debug) call PrintCurve(ordered)
00748 else
00749 ierr = IncrementCurve(lja,ljd)
00750 if(debug) print *,'PeanoM: After Position [1,1] ',pos
00751 endif
00752
00753
00754
00755
00756
00757 lma = MOD(ma+1,maxdim)
00758 lmd = -md
00759 lja = MOD(lma+1,maxdim)
00760 ljd = -lmd
00761
00762 if(ll .gt. 1) then
00763 if(debug) write(*,28) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00764 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00765 if(debug) call PrintCurve(ordered)
00766 else
00767 ierr = IncrementCurve(lja,ljd)
00768 if(debug) print *,'PeanoM: After Position [1,0] ',pos
00769 endif
00770
00771
00772
00773
00774 lma = ma
00775 lmd = md
00776 lja = ja
00777 ljd = jd
00778
00779 if(ll .gt. 1) then
00780 if(debug) write(*,29) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00781 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00782 if(debug) call PrintCurve(ordered)
00783 else
00784 ierr = IncrementCurve(lja,ljd)
00785 if(debug) print *,'PeanoM: After Position [2,0] ',pos
00786 endif
00787
00788 21 format('Call PeanoM Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
00789 22 format('Call PeanoM Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
00790 23 format('Call PeanoM Pos [0,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
00791 24 format('Call PeanoM Pos [1,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
00792 25 format('Call PeanoM Pos [2,2] Level ',i1,' at (',i2,',',i2,')',4(i3))
00793 26 format('Call PeanoM Pos [2,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
00794 27 format('Call PeanoM Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
00795 28 format('Call PeanoM Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
00796 29 format('Call PeanoM Pos [2,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
00797
00798
00799
00800
00801 end function PeanoM
00802
00803
00804
00805
00806
00807
00808 recursive function Hilbert(l,type,ma,md,ja,jd) result(ierr)
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823 integer(int_kind), intent(in) ::
00824 l,
00825 type,
00826 ma,
00827 md,
00828 ja,
00829 jd
00830
00831
00832
00833 integer(int_kind) :: ierr
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844 integer(int_kind) ::
00845 lma,
00846 lmd,
00847 lja,
00848 ljd,
00849 ltype,
00850 ll
00851
00852 logical :: debug = .FALSE.
00853
00854
00855 ll = l
00856 if(ll .gt. 1) ltype = fact%factors(ll-1)
00857
00858
00859
00860 lma = MOD(ma+1,maxdim)
00861 lmd = md
00862 lja = lma
00863 ljd = lmd
00864
00865 if(ll .gt. 1) then
00866 if(debug) write(*,21) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00867 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00868 if(debug) call PrintCurve(ordered)
00869 else
00870 ierr = IncrementCurve(lja,ljd)
00871 if(debug) print *,'Hilbert: After Position [0,0] ',pos
00872 endif
00873
00874
00875
00876
00877
00878 lma = ma
00879 lmd = md
00880 lja = lma
00881 ljd = lmd
00882 if(ll .gt. 1) then
00883 if(debug) write(*,22) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00884 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00885 if(debug) call PrintCurve(ordered)
00886 else
00887 ierr = IncrementCurve(lja,ljd)
00888 if(debug) print *,'Hilbert: After Position [0,1] ',pos
00889 endif
00890
00891
00892
00893
00894
00895 lma = ma
00896 lmd = md
00897 lja = MOD(ma+1,maxdim)
00898 ljd = -md
00899
00900 if(ll .gt. 1) then
00901 if(debug) write(*,23) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00902 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00903 if(debug) call PrintCurve(ordered)
00904 else
00905 ierr = IncrementCurve(lja,ljd)
00906 if(debug) print *,'Hilbert: After Position [1,1] ',pos
00907 endif
00908
00909
00910
00911
00912 lma = MOD(ma+1,maxdim)
00913 lmd = -md
00914 lja = ja
00915 ljd = jd
00916
00917 if(ll .gt. 1) then
00918 if(debug) write(*,24) ll-1,pos(0),pos(1),lma,lmd,lja,ljd
00919 ierr = GenCurve(ll-1,ltype,lma,lmd,lja,ljd)
00920 if(debug) call PrintCurve(ordered)
00921 else
00922 ierr = IncrementCurve(lja,ljd)
00923 if(debug) print *,'Hilbert: After Position [1,0] ',pos
00924 endif
00925
00926 21 format('Call Hilbert Pos [0,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
00927 22 format('Call Hilbert Pos [0,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
00928 23 format('Call Hilbert Pos [1,1] Level ',i1,' at (',i2,',',i2,')',4(i3))
00929 24 format('Call Hilbert Pos [1,0] Level ',i1,' at (',i2,',',i2,')',4(i3))
00930
00931
00932
00933
00934 end function hilbert
00935
00936
00937
00938
00939
00940
00941 function IncrementCurve(ja,jd) result(ierr)
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953 integer(int_kind) ::
00954 ja,
00955 jd
00956
00957
00958 integer(int_kind) :: ierr
00959
00960
00961
00962
00963 ordered(pos(0)+1,pos(1)+1) = vcnt
00964
00965
00966
00967
00968 vcnt = vcnt + 1
00969 pos(ja) = pos(ja) + jd
00970
00971 ierr = 0
00972
00973
00974
00975 end function IncrementCurve
00976
00977
00978
00979
00980
00981
00982 function log2( n)
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993 integer(int_kind), intent(in) :: n
00994
00995
00996
00997 integer(int_kind) :: log2
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007 integer(int_kind) :: tmp
01008
01009
01010
01011
01012 log2 = 1
01013 tmp =n
01014 do while (tmp/2 .ne. 1)
01015 tmp=tmp/2
01016 log2=log2+1
01017 enddo
01018
01019
01020
01021
01022 end function log2
01023
01024
01025
01026
01027
01028
01029 function IsLoadBalanced(nelem,npart)
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040 integer(int_kind), intent(in) ::
01041 nelem,
01042 npart
01043
01044
01045 logical :: IsLoadBalanced
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055 integer(int_kind) :: tmp1
01056
01057
01058 tmp1 = nelem/npart
01059
01060 if(npart*tmp1 == nelem ) then
01061 IsLoadBalanced=.TRUE.
01062 else
01063 IsLoadBalanced=.FALSE.
01064 endif
01065
01066
01067
01068
01069 end function IsLoadBalanced
01070
01071
01072
01073
01074
01075
01076 function GenCurve(l,type,ma,md,ja,jd) result(ierr)
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088 integer(int_kind), intent(in) ::
01089 l,
01090 type,
01091 ma,
01092 md,
01093 ja,
01094 jd
01095
01096
01097
01098 integer(int_kind) :: ierr
01099
01100
01101
01102
01103
01104
01105
01106
01107 if(type == 2) then
01108 ierr = Hilbert(l,type,ma,md,ja,jd)
01109 elseif ( type == 3) then
01110 ierr = PeanoM(l,type,ma,md,ja,jd)
01111 elseif ( type == 5) then
01112 ierr = Cinco(l,type,ma,md,ja,jd)
01113 endif
01114
01115
01116
01117
01118 end function GenCurve
01119
01120
01121 function FirstFactor(fac) result(res)
01122 type (factor_t) :: fac
01123 integer :: res
01124 logical :: found
01125 integer (int_kind) :: i
01126
01127 found = .false.
01128 i=1
01129 do while (i<=fac%numfact .and. (.not. found))
01130 if(fac%used(i) == 0) then
01131 res = fac%factors(i)
01132 found = .true.
01133 endif
01134 i=i+1
01135 enddo
01136
01137 end function FirstFactor
01138
01139 function FindandMark(fac,val,f2) result(found)
01140 type (factor_t) :: fac
01141 integer :: val
01142 logical :: found
01143 logical :: f2
01144 integer (int_kind) :: i
01145
01146 found = .false.
01147 i=1
01148 do while (i<=fac%numfact .and. (.not. found))
01149 if(fac%used(i) == 0) then
01150 if(fac%factors(i) .eq. val) then
01151 if(f2) then
01152 fac%used(i) = 1
01153 found = .true.
01154 else if( .not. f2) then
01155 fac%used(i) = -1
01156 found = .false.
01157 endif
01158 endif
01159 endif
01160 i=i+1
01161 enddo
01162
01163 end function FindandMark
01164
01165
01166 subroutine MatchFactor(fac1,fac2,val,found)
01167 type (factor_t) :: fac1
01168 type (factor_t) :: fac2
01169 integer :: val
01170 integer :: val1
01171 logical :: found
01172 logical :: tmp
01173
01174 found = .false.
01175
01176 val1 = FirstFactor(fac1)
01177
01178 found = FindandMark(fac2,val1,.true.)
01179 tmp = FindandMark(fac1,val1,found)
01180 if (found) then
01181 val = val1
01182 else
01183 val = 1
01184 endif
01185
01186 end subroutine MatchFactor
01187
01188 function ProdFactor(fac) result(res)
01189
01190 type (factor_t) :: fac
01191 integer :: res
01192 integer (int_kind) :: i
01193
01194 res = 1
01195 do i=1,fac%numfact
01196 if(fac%used(i) <= 0) then
01197 res = res * fac%factors(i)
01198 endif
01199 enddo
01200
01201 end function ProdFactor
01202
01203 subroutine PrintFactor(msg,fac)
01204
01205
01206 character(len=*) :: msg
01207 type (factor_t) :: fac
01208 integer (int_kind) :: i
01209
01210 write(*,*) ' '
01211 write(*,*) 'PrintFactor: ',msg
01212 write(*,*) (fac%factors(i),i=1,fac%numfact)
01213 write(*,*) (fac%used(i),i=1,fac%numfact)
01214
01215
01216 end subroutine PrintFactor
01217
01218
01219
01220
01221
01222
01223 function Factor(num) result(res)
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235 integer(int_kind), intent(in) :: num
01236
01237
01238
01239 type (factor_t) :: res
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249 integer(int_kind) ::
01250 tmp,tmp2,tmp3,tmp5
01251 integer(int_kind) :: i,n
01252 logical :: found
01253
01254
01255
01256
01257 tmp = num
01258 tmp2 = log2(num)
01259 allocate(res%factors(tmp2))
01260 allocate(res%used(tmp2))
01261
01262 res%used = 0
01263 n=0
01264
01265
01266
01267
01268
01269 found=.TRUE.
01270 do while (found)
01271 found = .FALSE.
01272 tmp2 = tmp/2
01273 if( tmp2*2 == tmp ) then
01274 n = n + 1
01275 res%factors(n) = 2
01276 found = .TRUE.
01277 tmp = tmp2
01278 endif
01279 enddo
01280
01281
01282
01283
01284 found=.TRUE.
01285 do while (found)
01286 found = .FALSE.
01287 tmp3 = tmp/3
01288 if( tmp3*3 == tmp ) then
01289 n = n + 1
01290 res%factors(n) = 3
01291 found = .TRUE.
01292 tmp = tmp3
01293 endif
01294 enddo
01295
01296
01297
01298
01299 found=.TRUE.
01300 do while (found)
01301 found = .FALSE.
01302 tmp5 = tmp/5
01303 if( tmp5*5 == tmp ) then
01304 n = n + 1
01305 res%factors(n) = 5
01306 found = .TRUE.
01307 tmp = tmp5
01308 endif
01309 enddo
01310
01311
01312
01313
01314
01315
01316 tmp=1
01317 do i=1,n
01318 tmp = tmp * res%factors(i)
01319 enddo
01320 if(tmp == num) then
01321 res%numfact = n
01322 else
01323 res%numfact = -1
01324 endif
01325
01326
01327
01328 end function Factor
01329
01330
01331
01332
01333
01334
01335 function IsFactorable(n)
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347 integer(int_kind), intent(in) :: n
01348
01349
01350 logical :: IsFactorable
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360 type (factor_t) :: fact
01361
01362 fact = Factor(n)
01363 if(fact%numfact .ne. -1) then
01364 IsFactorable = .TRUE.
01365 else
01366 IsFactorable = .FALSE.
01367 endif
01368
01369
01370
01371
01372 end function IsFactorable
01373
01374
01375
01376
01377
01378
01379 subroutine map(l)
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389 integer(int_kind) :: l
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400 integer(int_kind) ::
01401 d,
01402 type,
01403 ierr
01404
01405 d = SIZE(pos)
01406
01407 pos=0
01408 maxdim=d
01409 vcnt=0
01410
01411 type = fact%factors(l)
01412 ierr = GenCurve(l,type,0,1,0,1)
01413
01414
01415
01416
01417
01418 end subroutine map
01419
01420
01421
01422
01423
01424
01425 subroutine PrintCurve(Mesh)
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437 integer(int_kind), intent(in), target :: Mesh(:,:)
01438
01439
01440
01441
01442
01443
01444
01445
01446 integer(int_kind) ::
01447 gridsize,
01448 i
01449
01450
01451
01452 gridsize = SIZE(Mesh,dim=1)
01453
01454 if(gridsize == 2) then
01455 write (*,*) "A Level 1 Hilbert Curve:"
01456 write (*,*) "------------------------"
01457 do i=1,gridsize
01458 write(*,2) Mesh(1,i),Mesh(2,i)
01459 enddo
01460 else if(gridsize == 3) then
01461 write (*,*) "A Level 1 Peano Meandering Curve:"
01462 write (*,*) "---------------------------------"
01463 do i=1,gridsize
01464 write(*,3) Mesh(1,i),Mesh(2,i),Mesh(3,i)
01465 enddo
01466 else if(gridsize == 4) then
01467 write (*,*) "A Level 2 Hilbert Curve:"
01468 write (*,*) "------------------------"
01469 do i=1,gridsize
01470 write(*,4) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i)
01471 enddo
01472 else if(gridsize == 5) then
01473 write (*,*) "A Level 1 Cinco Curve:"
01474 write (*,*) "------------------------"
01475 do i=1,gridsize
01476 write(*,5) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i),Mesh(5,i)
01477 enddo
01478 else if(gridsize == 6) then
01479 write (*,*) "A Level 1 Hilbert and Level 1 Peano Curve:"
01480 write (*,*) "------------------------------------------"
01481 do i=1,gridsize
01482 write(*,6) Mesh(1,i),Mesh(2,i),Mesh(3,i), &
01483 Mesh(4,i),Mesh(5,i),Mesh(6,i)
01484 enddo
01485 else if(gridsize == 8) then
01486 write (*,*) "A Level 3 Hilbert Curve:"
01487 write (*,*) "------------------------"
01488 do i=1,gridsize
01489 write(*,8) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
01490 Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i)
01491 enddo
01492 else if(gridsize == 9) then
01493 write (*,*) "A Level 2 Peano Meandering Curve:"
01494 write (*,*) "---------------------------------"
01495 do i=1,gridsize
01496 write(*,9) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
01497 Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
01498 Mesh(9,i)
01499 enddo
01500 else if(gridsize == 10) then
01501 write (*,*) "A Level 1 Hilbert and Level 1 Cinco Curve:"
01502 write (*,*) "---------------------------------"
01503 do i=1,gridsize
01504 write(*,10) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
01505 Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
01506 Mesh(9,i),Mesh(10,i)
01507 enddo
01508 else if(gridsize == 12) then
01509 write (*,*) "A Level 2 Hilbert and Level 1 Peano Curve:"
01510 write (*,*) "------------------------------------------"
01511 do i=1,gridsize
01512 write(*,12) Mesh(1,i),Mesh(2,i), Mesh(3,i), Mesh(4,i), &
01513 Mesh(5,i),Mesh(6,i), Mesh(7,i), Mesh(8,i), &
01514 Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i)
01515 enddo
01516 else if(gridsize == 15) then
01517 write (*,*) "A Level 1 Peano and Level 1 Cinco Curve:"
01518 write (*,*) "------------------------"
01519 do i=1,gridsize
01520 write(*,15) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
01521 Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
01522 Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i), &
01523 Mesh(13,i),Mesh(14,i),Mesh(15,i)
01524 enddo
01525 else if(gridsize == 16) then
01526 write (*,*) "A Level 4 Hilbert Curve:"
01527 write (*,*) "------------------------"
01528 do i=1,gridsize
01529 write(*,16) Mesh(1,i),Mesh(2,i),Mesh(3,i),Mesh(4,i), &
01530 Mesh(5,i),Mesh(6,i),Mesh(7,i),Mesh(8,i), &
01531 Mesh(9,i),Mesh(10,i),Mesh(11,i),Mesh(12,i), &
01532 Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i)
01533 enddo
01534 else if(gridsize == 18) then
01535 write (*,*) "A Level 1 Hilbert and Level 2 Peano Curve:"
01536 write (*,*) "------------------------------------------"
01537 do i=1,gridsize
01538 write(*,18) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
01539 Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
01540 Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
01541 Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
01542 Mesh(17,i),Mesh(18,i)
01543 enddo
01544 else if(gridsize == 20) then
01545 write (*,*) "A Level 2 Hilbert and Level 1 Cinco Curve:"
01546 write (*,*) "------------------------------------------"
01547 do i=1,gridsize
01548 write(*,20) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
01549 Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
01550 Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
01551 Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
01552 Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i)
01553 enddo
01554 else if(gridsize == 24) then
01555 write (*,*) "A Level 3 Hilbert and Level 1 Peano Curve:"
01556 write (*,*) "------------------------------------------"
01557 do i=1,gridsize
01558 write(*,24) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
01559 Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
01560 Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
01561 Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
01562 Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
01563 Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i)
01564 enddo
01565 else if(gridsize == 25) then
01566 write (*,*) "A Level 2 Cinco Curve:"
01567 write (*,*) "------------------------------------------"
01568 do i=1,gridsize
01569 write(*,25) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
01570 Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
01571 Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
01572 Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
01573 Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
01574 Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), &
01575 Mesh(25,i)
01576 enddo
01577 else if(gridsize == 27) then
01578 write (*,*) "A Level 3 Peano Meandering Curve:"
01579 write (*,*) "---------------------------------"
01580 do i=1,gridsize
01581 write(*,27) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
01582 Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
01583 Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
01584 Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
01585 Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
01586 Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), &
01587 Mesh(25,i),Mesh(26,i),Mesh(27,i)
01588 enddo
01589 else if(gridsize == 32) then
01590 write (*,*) "A Level 5 Hilbert Curve:"
01591 write (*,*) "------------------------"
01592 do i=1,gridsize
01593 write(*,32) Mesh(1,i), Mesh(2,i), Mesh(3,i), Mesh(4,i), &
01594 Mesh(5,i), Mesh(6,i), Mesh(7,i), Mesh(8,i), &
01595 Mesh(9,i), Mesh(10,i),Mesh(11,i),Mesh(12,i), &
01596 Mesh(13,i),Mesh(14,i),Mesh(15,i),Mesh(16,i), &
01597 Mesh(17,i),Mesh(18,i),Mesh(19,i),Mesh(20,i), &
01598 Mesh(21,i),Mesh(22,i),Mesh(23,i),Mesh(24,i), &
01599 Mesh(25,i),Mesh(26,i),Mesh(27,i),Mesh(28,i), &
01600 Mesh(29,i),Mesh(30,i),Mesh(31,i),Mesh(32,i)
01601 enddo
01602 endif
01603 2 format('|',2(i2,'|'))
01604 3 format('|',3(i2,'|'))
01605 4 format('|',4(i2,'|'))
01606 5 format('|',5(i2,'|'))
01607 6 format('|',6(i2,'|'))
01608 8 format('|',8(i2,'|'))
01609 9 format('|',9(i2,'|'))
01610 10 format('|',10(i2,'|'))
01611 12 format('|',12(i3,'|'))
01612 15 format('|',15(i3,'|'))
01613 16 format('|',16(i3,'|'))
01614 18 format('|',18(i3,'|'))
01615 20 format('|',20(i3,'|'))
01616 24 format('|',24(i3,'|'))
01617 25 format('|',25(i3,'|'))
01618 27 format('|',27(i3,'|'))
01619 32 format('|',32(i4,'|'))
01620
01621
01622
01623
01624 end subroutine PrintCurve
01625
01626
01627
01628
01629
01630
01631 subroutine GenSpaceCurve(Mesh)
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642 integer(int_kind), target,intent(inout) ::
01643 Mesh(:,:)
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653 integer(int_kind) ::
01654 level,
01655 dim
01656
01657 integer(int_kind) :: gridsize
01658
01659
01660
01661
01662
01663
01664 dim = 2
01665 gridsize = SIZE(Mesh,dim=1)
01666 fact = factor(gridsize)
01667 level = fact%numfact
01668
01669 if(verbose) print *,'GenSpacecurve: level is ',level
01670 allocate(ordered(gridsize,gridsize))
01671
01672
01673
01674
01675 allocate(pos(0:dim-1))
01676
01677
01678
01679
01680 ordered(:,:) = 0
01681
01682 call map(level)
01683
01684 Mesh(:,:) = ordered(:,:)
01685
01686 deallocate(pos,ordered)
01687
01688
01689
01690
01691 end subroutine GenSpaceCurve
01692
01693 recursive subroutine qsort(a)
01694
01695 integer, intent(inout) :: a(:)
01696 integer :: split
01697
01698 if(SIZE(a) > 1) then
01699 call partition(a,split)
01700 call qsort(a(:split-1))
01701 call qsort(a(split:))
01702 endif
01703
01704 end subroutine qsort
01705
01706 subroutine partition(a,marker)
01707
01708 INTEGER, INTENT(IN OUT) :: a(:)
01709 INTEGER, INTENT(OUT) :: marker
01710 INTEGER :: left, right, pivot, temp
01711
01712 pivot = (a(1) + a(size(a))) / 2
01713 left = 0
01714 right = size(a) + 1
01715
01716 DO WHILE (left < right)
01717 right = right - 1
01718 DO WHILE (a(right) > pivot)
01719 right = right-1
01720 END DO
01721 left = left + 1
01722 DO WHILE (a(left) < pivot)
01723 left = left + 1
01724 END DO
01725 IF (left < right) THEN
01726 temp = a(left)
01727 a(left) = a(right)
01728 a(right) = temp
01729 END IF
01730 END DO
01731
01732 IF (left == right) THEN
01733 marker = left + 1
01734 ELSE
01735 marker = left
01736 END IF
01737
01738 end subroutine partition
01739
01740
01741
01742 end module ice_spacecurve
01743
01744