1!********************************************************** 2! RAMG tools, functions needed in ramg_driver and 3! ramg_interface. 4! Module defined in ramg_data.f 5! 6! ramg_calcIAI 7! ramg_calcIvFC/CF 8! ramg_direct 9! ramg_allocate 10! ramg_deallocate 11! ramg_dump 12! 13!*********************************************************** 14 15 16!!*********************************************************** 17! ramg_calcIvCF 18! calculate V coarse to fine, 19! v = I * VC 20!*********************************************************** 21 subroutine ramg_calcIvCF(level1,level2,VC,vf, 22 & ilwork,BC,iBC,iper) 23 use ramg_data 24 include "common.h" 25 include "mpif.h" 26 27 !implicit none 28 integer,intent(in) :: level1,level2 29 real(kind=8),intent(inout),dimension(amg_nshg(level1)) :: vf 30 real(kind=8),intent(in),dimension(amg_nshg(level2)) :: VC 31 32 integer,intent(in),dimension(nlwork) :: ilwork 33 integer,intent(in),dimension(nshg) :: iBC,iper 34 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 35 36 integer :: i,j,k,p 37 38 real(kind=8) :: cpusec(2) 39 40 call cpu_time(cpusec(1)) 41 42 vf = 0 43 do i=1,amg_nshg(level1) 44 do k=I_cf_colm(level1)%p(i),I_cf_colm(level1)%p(i+1)-1 45 j = I_cf_rowp(level1)%p(k) 46 vf(i) = vf(i) + VC(j)*I_cf(level1)%p(k) 47 enddo 48 enddo 49 50 if ( level1 .eq. -1 ) then 51 call commIn(vf,ilwork,1,iper,iBC,BC) 52 call commOut(vf,ilwork,1,iper,iBC,BC) 53 end if 54 55 call cpu_time(cpusec(2)) 56 ramg_time(13)=ramg_time(13)+cpusec(2)-cpusec(1) 57 58 end subroutine ! ramg_calcIvCF 59 60!*********************************************************** 61! ramg_calcIvFC 62! calculate v fine to coarse, 63! VC = IT * v 64!*********************************************************** 65 subroutine ramg_calcIvFC(level1,level2,vf,VC) 66 use ramg_data 67 implicit none 68 integer,intent(in) :: level1,level2 69 real(kind=8),intent(in),dimension(amg_nshg(level1)) :: vf 70 real(kind=8),intent(inout),dimension(amg_nshg(level2)) :: VC 71 72 integer :: i,j,k,p 73 real(kind=8) :: cpusec(2) 74 75 call cpu_time(cpusec(1)) 76 77 78 VC = 0 79 do i=1,amg_nshg(level2) 80 do k=I_fc_colm(level1)%p(i),I_fc_colm(level1)%p(i+1)-1 81 j = I_fc_rowp(level1)%p(k) 82 VC(i) = VC(i) + vf(j)*I_fc(level1)%p(k) 83 enddo 84 enddo 85 call cpu_time(cpusec(2)) 86 ramg_time(13)=ramg_time(13)+cpusec(2)-cpusec(1) 87 88 end subroutine ! ramg_calcIvFC 89 90!!*********************************************************** 91! ramg_calcSvCF 92! calculate V coarse to fine, 93! v = S * VC 94!*********************************************************** 95 96!*********************************************************** 97! ramg_calcIvFC 98! calculate v fine to coarse, 99! VC = ST * v 100!*********************************************************** 101 102!************************************************************ 103! ramg_calcAv 104! calculate u = A*v 105!************************************************************ 106 subroutine ramg_calcAv(gcolm,growp,glhs,gnshg,gnnz_tot, 107 & u,v,gcomm) 108 use ramg_data 109 integer,intent(in) :: gnshg,gnnz_tot 110 integer,intent(in),dimension(gnshg+1) :: gcolm 111 integer,intent(in),dimension(gnnz_tot) :: growp 112 real(kind=8),intent(in),dimension(gnnz_tot) :: glhs 113 real(kind=8),intent(in),dimension(gnshg) :: v 114 real(kind=8),intent(inout),dimension(gnshg) :: u 115 integer,intent(in) :: gcomm 116 117 integer :: i,j,k,p 118 119 u = 0 120 do i=1,gnshg 121 do k=gcolm(i),gcolm(i+1)-1 122 j = growp(k) 123 u(i) = u(i)+glhs(k)*v(j) 124 enddo 125 enddo 126 127! if (gcomm.eq.1) then 128! end if 129 130 end subroutine ! ramg_calcAv 131 132!************************************************************ 133! ramg_calcAv_g 134! calculate u = A*v 135! Globally: commout, do product, commin, (zeroout) 136!************************************************************ 137 subroutine ramg_calcAv_g(level,u,v,colm,rowp,lhsK,lhsP, 138 & ilwork,BC,iBC,iper,gcomm) 139 use ramg_data 140 include "common.h" 141 include "mpif.h" 142 include "auxmpi.h" 143 144 integer,intent(in),dimension(nlwork) :: ilwork 145 integer,intent(in),dimension(nshg) :: iBC,iper 146 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 147 integer,intent(in),dimension(nshg+1) :: colm 148 integer,intent(in),dimension(nnz_tot) :: rowp 149 150 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 151 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 152 integer,intent(in) :: level 153 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: v 154 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: u 155 integer,intent(in) :: gcomm 156 157 integer :: i,j,k,p,mi,mj 158 real(kind=8) :: cpusec(2) 159 160 call cpu_time(cpusec(1)) 161 162 IF (level.eq.1) THEN 163 !IF (.FALSE.) THEN 164 call ramg_PPEAp(u,v,colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 165 ELSE 166 if (gcomm.eq.1) then 167 call ramg_commOut(v,level,ilwork,1,iper,iBC,BC) 168 endif 169 u = 0 170 do i=1,amg_nshg(level) 171 mi = amg_paramap(level)%p(i) 172 do k=amg_A_colm(level)%p(i),amg_A_colm(level)%p(i+1)-1 173 j = amg_A_rowp(level)%p(k) 174 mj = amg_paramap(level)%p(j) 175 if (.not.( (mi.eq.mj).and.(mi.lt.0) )) then 176 u(i) = u(i)+amg_A_lhs(level)%p(k,1)*v(j) 177 endif 178 enddo 179 enddo 180 if (gcomm.eq.1) then 181 call ramg_commIn(u,level,ilwork,1,iper,iBC,BC) 182 endif 183 ENDIF 184 call cpu_time(cpusec(2)) 185 if (level.eq.1) then 186 ramg_time(11) = ramg_time(11) + cpusec(2)-cpusec(1) 187 else 188 ramg_time(12) = ramg_time(12) + cpusec(2)-cpusec(1) 189 endif 190 191 end subroutine ! ramg_calcAv_g 192 193!************************************************************ 194! ramg_calcAv_b: same as calcAv_g, but only apply 195! on boundary nodes. 196! calculate u = A*v 197! Globally: commout, do product, commin, (zeroout) 198!************************************************************ 199 subroutine ramg_calcAv_b(level,u,v, 200 & ilwork,BC,iBC,iper,gcomm,diag) 201 use ramg_data 202 include "common.h" 203 include "mpif.h" 204 include "auxmpi.h" 205 206 integer,intent(in),dimension(nlwork) :: ilwork 207 integer,intent(in),dimension(nshg) :: iBC,iper 208 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 209 210 integer,intent(in) :: level 211 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: v 212 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: u 213 integer,intent(in) :: gcomm 214 integer,intent(in) :: diag !=1: NOT include diagonal A_(i,i) 215 !=0: include diagonal 216 217 integer :: i,j,k,p,mi,mj,mk 218 real(kind=8) :: cpusec(2) 219 220 call cpu_time(cpusec(1)) 221 222 if (gcomm.eq.1) then 223 call ramg_commOut(v,level,ilwork,1,iper,iBC,BC) 224 endif 225 u = 0 226 do i=1,amg_nshg(level) 227 mi = amg_paramap(level)%p(i) 228 mk = amg_paraext(level)%p(i) 229 IF (mk.ne.(myrank+1)) THEN ! only treat boundary nodes 230 do k=amg_A_colm(level)%p(i)+diag,amg_A_colm(level)%p(i+1)-1 231 j = amg_A_rowp(level)%p(k) 232 mj = amg_paramap(level)%p(j) 233 if (.not.( (mi.eq.mj).and.(mi.lt.0) )) then 234 u(i) = u(i)+amg_A_lhs(level)%p(k,1)*v(j) 235 endif 236 enddo 237 ELSE 238 u(i) = v(i) 239 ENDIF 240 enddo 241 if (gcomm.eq.1) then 242 call ramg_commIn(u,level,ilwork,1,iper,iBC,BC) 243 endif 244 245 end subroutine ! ramg_calcAv_b 246 247 248!************************************************* 249! check_paracomm 250!************************************************* 251 subroutine ramg_check_paracomm(ilwork,BC,iBC,iper) 252 use ramg_data 253 include "common.h" 254 include "mpif.h" 255 include "auxmpi.h" 256 257 integer,intent(in),dimension(nlwork) :: ilwork 258 integer,intent(in),dimension(nshg) :: iBC,iper 259 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 260 261 integer i,j 262 type(r1d),dimension(ramg_levelx) :: tarray 263 character :: fname*80 264 265 do i=1,ramg_levelx 266 allocate(tarray(i)%p(amg_nshg(i))) 267 do j=1,amg_nshg(i) 268 call random_number(tarray(i)%p) 269 enddo 270 enddo 271 272 do i=1,ramg_levelx 273 write(fname,'((A8)(I1))')'tarray_a',i 274 call ramg_dump(tarray(i)%p,amg_nshg(i),fname) 275 enddo 276 277 do i=1,ramg_levelx 278 call ramg_commIn(tarray(i)%p,i,ilwork,1,iper,iBC,BC) 279 write(fname,'((A8)(I1))')'tarray_i',i 280 call ramg_dump(tarray(i)%p,amg_nshg(i),fname) 281 enddo 282 283 do i=1,ramg_levelx 284 call ramg_commOut(tarray(i)%p,i,ilwork,1,iper,iBC,BC) 285 write(fname,'((A8)(I1))')'tarray_b',i 286 call ramg_dump(tarray(i)%p,amg_nshg(i),fname) 287 enddo 288 289! call commOut(tarray(1)%p,ilwork,1,iper,iBC,BC) 290 291! do i=1,ramg_levelx 292! write(fname,'((A8)(I1))')'tarray_i',i 293! call ramg_dump_rn_map(tarray(i)%p,amg_nshg(i),1,fname) 294! enddo 295 296 297 do i=1,ramg_levelx 298 deallocate(tarray(i)%p) 299 enddo 300 301 end subroutine ! ramg_check_paracomm 302 303!************************************************** 304! ramg_zeroOut: zero out slave values 305!************************************************** 306 subroutine ramg_zeroOut(u,ilwork,BC,iBC,iper) 307 include "common.h" 308 include "auxmpi.h" 309 include "mpif.h" 310 311 real(kind=8),dimension(nshg),intent(inout) :: u 312 integer,intent(in),dimension(nlwork) :: ilwork 313 integer,intent(in),dimension(nshg) :: iBC,iper 314 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 315 316 integer i,j,k,p 317 integer :: itag,iacc,iother,numseg,isgbeg,itkbeg 318 integer :: numtask,lenseg 319 320 if (numpe.lt.2) return 321 !call commOut(u,ilwork,1,iper,iBC,BC) 322 numtask = ilwork(1) 323 itkbeg = 1 324 do i=1,numtask 325 iacc = ilwork(itkbeg+2) 326 numseg = ilwork(itkbeg+4) 327 if (iacc.eq.0) then 328 do j=1,numseg 329 isgbeg = ilwork(itkbeg+3+j*2) 330 lenseg = ilwork(itkbeg+4+j*2) 331 isgend = isgbeg + lenseg -1 332 u(isgbeg:isgend) = 0 333 enddo 334 endif 335 itkbeg = itkbeg+4+2*numseg 336 enddo 337 338 end subroutine ! ramg_zeroOut 339 340!************************************************ 341! ramg_freeBC: set "fine" on boundary nodes 342!************************************************ 343 subroutine ramg_freeBC(amg_F,ilwork,BC,iBC,iper) 344 use ramg_data 345 include "common.h" 346 include "auxmpi.h" 347 include "mpif.h" 348 integer,intent(inout),dimension(nshg) :: amg_F 349 integer,intent(in),dimension(nlwork) :: ilwork 350 integer,intent(in),dimension(nshg) :: iBC,iper 351 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 352 integer :: itag,iacc,iother,numseg,isgbeg,itkbeg,numtask 353 integer :: i,p 354 integer,dimension(nshg,2) :: tmpout 355 356 character :: filename*80 357 integer :: nentry,nshg,iglobal,ilocal,icf,iproc 358 359 return 360 361 if (numpe.ge.2) then 362 363 numtask = ilwork(1) 364 itkbeg = 1 365 do i=1,numtask 366 itag = ilwork(itkbeg+1) 367 iacc = ilwork(itkbeg+2) 368 iother = ilwork(itkbeg+3) 369 numseg = ilwork(itkbeg+4) 370 isgbeg = ilwork(itkbeg+5) 371 do j=1,numseg 372 p = ilwork(itkbeg+3+j*2) 373 amg_F(p) = 2 374 enddo 375 enddo 376 377 tmpout(:,1) = ncorp_map 378 tmpout(:,2) = amg_F 379 380 call ramg_dump_i(tmpout,nshg,2,'CFsplit ') 381 else ! DEBUGGING PART, READ 2-PROC CASE INTO 1-proc 382 filename = "cfs.dat" 383 filename = trim(filename) 384 open(unit=999,file=filename,status='unknown') 385 read(999,*)nentry 386 do i=1,nentry 387 read(999,*)iglobal,ilocal,icf,iproc 388 !amg_F(iglobal) = icf 389 enddo 390 close(999) 391 end if 392 393 end subroutine ! ramg_freeBC 394 395 396!********************************************** 397! ramg_read_map: 398! read in ncorp array from local to global 399! mapping 400!********************************************** 401 subroutine ramg_prep_reduce_map 402 use ramg_data 403 include "common.h" 404 include "mpif.h" 405 include "auxmpi.h" 406 407 character*5 cname 408 character*255 fname1 409 integer igeom,map_nshg,i 410 411 allocate(ncorp_map(nshg)) 412 if (numpe.lt.2) then! construct dummy-array 413 do i=1,nshg 414 ncorp_map(i) = i 415 enddo 416 else 417 fname1='geombc.dat' 418 fname1=trim(fname1)//cname(myrank+1) 419 420 call openfile(fname1,"read",igeom) 421 !write(*,*)'mcheck:',myrank,fname1,igeom 422 423 fname1="mode number map from partition to global" 424 ione=1 425 call readheader(igeom, 426 & 'mode number map from partition to global' // char(0), 427 & map_nshg,ione,"integer",iotype) 428 429 call readdatablock(igeom, 430 & 'mode number map from partition to global' // char(0), 431 & ncorp_map,map_nshg, "integer",iotype) 432 433 call closefile(igeom,"read") 434 435 endif 436 437 end subroutine ! ramg_read_map 438 439!*********************************************** 440! ramg_check_diag 441!*********************************************** 442 subroutine ramg_check_diag(colm,rowp,lhs,nshg,nnz_tot) 443 implicit none 444 integer,intent(in) :: nshg,nnz_tot 445 integer,intent(in),dimension(nshg+1) :: colm 446 integer,intent(in),dimension(nnz_tot) :: rowp 447 real(kind=8),intent(in),dimension(nnz_tot) :: lhs 448 449 integer :: i,j,p 450 real(kind=8) :: diagrow 451 logical :: diagokay 452 453 diagokay = .true. 454 455 do i=1,nshg 456 p = colm(i) 457 if (rowp(p).ne.i) then 458 write(*,*)'matrix first row entry is not diagonal',i 459 diagokay = .false. 460 end if 461 diagrow = lhs(p) 462 if (diagrow.lt.0) then 463 write(*,*)'matrix diagonal < 0 at',i,diagrow 464 diagokay = .false. 465 end if 466 do j = colm(i)+1,colm(i+1)-1 467 p = rowp(j) 468 if (lhs(p).gt.diagrow) then 469 write(*,*)'matrix not diagonal dominant at row',i 470 diagokay = .false. 471 write(*,*)'diag:',diagrow,p,lhs(p) 472 exit 473 end if 474 end do 475 enddo 476 477 if (diagokay) then 478 write(*,*)'matrix check diagonal dominance okay' 479 end if 480 481 end subroutine 482 483 subroutine ramg_dot_p(level,v,u,redun,norm) 484 use ramg_data 485 include "common.h" 486 include "mpif.h" 487 include "auxmpi.h" 488 489 integer :: redun,level 490 real(kind=8),intent(in),dimension(amg_nshg(level),redun) ::v,u 491 real(kind=8),intent(inout) :: norm 492 integer :: i,k 493 real(kind=8) :: normt 494 normt = 0 495 do i=1,amg_nshg(level) 496 do k=1,redun 497 normt = normt + v(i,k)*u(i,k) 498 enddo 499 enddo 500 if (numpe.gt.1) then 501 call MPI_AllReduce(normt,norm,1,MPI_DOUBLE_PRECISION,MPI_SUM, 502 & MPI_COMM_WORLD,ierr) 503 else 504 norm=normt 505 endif 506 end subroutine !ramg_dot_p 507 508!************************************************************ 509! ramg_L2_norm 510! calculate norm = | r | 511!************************************************************ 512 subroutine ramg_L2_norm(nshg,v,norm) 513 implicit none 514 integer,intent(in) :: nshg 515 real(kind=8),intent(in),dimension(nshg) :: v 516 real(kind=8),intent(inout) :: norm 517 integer :: i 518 norm = 0 519 do i=1,nshg 520 norm = norm + v(i)*v(i) 521 enddo 522 norm = sqrt(norm) 523 end subroutine !ramg_L2_norm 524 525 subroutine ramg_L2_norm_p(level,v,vflow,norm) 526 use ramg_data 527 include "common.h" 528 include "mpif.h" 529 include "auxmpi.h" 530 531 integer,intent(in) :: vflow,level 532 real(kind=8),intent(in),dimension(amg_nshg(level),vflow) :: v 533 real(kind=8),intent(inout) :: norm 534 integer :: i,k 535 real(kind=8) :: normt 536 normt = 0 537 do i=1,amg_nshg(level) 538 do k=1,vflow 539 normt = normt + v(i,k)*v(i,k) 540 enddo 541 enddo 542 if (numpe.gt.1) then 543 call MPI_AllReduce(normt,norm,1,MPI_DOUBLE_PRECISION,MPI_SUM, 544 & MPI_COMM_WORLD,ierr) 545 else 546 norm = normt 547 endif 548 norm = sqrt(norm) 549 end subroutine !ramg_L2_norm_p 550 551!************************************************************ 552! ramg_allocate & deallocate 553! (de)allocate memory for level l, lhs matrix, rhs vector 554!************************************************************ 555 subroutine ramg_allocate( 556 & level,lnshg,lnnz_tot,n_sol) 557 use ramg_data 558 implicit none 559 560 integer,intent(in) :: level,lnshg,lnnz_tot 561 integer,intent(in) :: n_sol 562 integer :: mem_err,mem_err_s 563 mem_err_s = 0 564 if (lnnz_tot.ne.0) then ! zero means manual alloc later 565 amg_nnz(level) = lnnz_tot 566 allocate(amg_A_lhs(level)%p(amg_nnz(level),n_sol), 567 & stat=mem_err) 568 mem_err_s = mem_err_s + mem_err 569 allocate(amg_A_rowp(level)%p(amg_nnz(level)), 570 & stat=mem_err) 571 mem_err_s = mem_err_s + mem_err 572 endif 573 if (lnshg.ne.0) then 574 amg_nshg(level) = lnshg 575 allocate(amg_A_colm(level)%p(amg_nshg(level)+1), 576 & stat=mem_err) 577 mem_err_s = mem_err_s + mem_err 578 allocate(amg_A_rhs(level)%p(amg_nshg(level)), 579 & stat=mem_err) 580 mem_err_s = mem_err_s + mem_err 581 allocate(amg_ppeDiag(level)%p(amg_nshg(level)), 582 & stat=mem_err) 583 endif 584 mem_err_s = mem_err_s + mem_err 585 if (mem_err_s .ne. 0 ) then 586 write(6,7001)level 587 end if 5887001 format(/' **** ramg: Allocation error at level',i5) 589 ! zero out 590 if (lnnz_tot.ne.0) then 591 amg_A_lhs(level)%p(:,:) = 0 592 amg_A_rowp(level)%p(:) = 0 593 endif 594 if (lnshg.ne.0) then 595 amg_A_colm(level)%p(:) = 0 596 amg_A_rhs(level)%p(:) = 0 597 endif 598 return 599 600 end subroutine ! ramg_allocate 601 602 subroutine ramg_deallocate(level) 603 604 use ramg_data 605 include "common.h" 606 607 integer,intent(inout) :: level 608 integer :: mem_err,mem_err_s,i 609 610 611 if (level.eq.1) then 612 613 if (maxnev.gt.0) then 614 deallocate(ramg_ggb_ev) 615 deallocate(ramg_ggb_eA) 616 deallocate(cmtxA) 617 deallocate(cindx) 618 endif 619 620 621 if (iamg_reduce.gt.1) then 622 deallocate(reducemap) 623 deallocate(rmap1d) 624 endif 625 626 end if 627 628 629 level = abs(level) 630 631 mem_err_s = 0 632 if (associated(amg_A_lhs(level)%p)) then 633 deallocate(amg_A_lhs(level)%p, 634 & stat=mem_err) 635 mem_err_s = mem_err_s + mem_err 636 endif 637 if (associated(amg_A_rowp(level)%p)) then 638 deallocate(amg_A_rowp(level)%p, 639 & stat=mem_err) 640 mem_err_s = mem_err_s + mem_err 641 endif 642 if (associated(amg_A_colm(level)%p)) then 643 deallocate(amg_A_colm(level)%p, 644 & stat=mem_err) 645 mem_err_s = mem_err_s + mem_err 646 endif 647 if (associated(amg_ppeDiag(level)%p)) then 648 deallocate(amg_ppeDiag(level)%p, 649 & stat=mem_err) 650 mem_err_s = mem_err_s + mem_err 651 endif 652 if (associated(amg_A_rhs(level)%p)) then 653 deallocate(amg_A_rhs(level)%p, 654 & stat=mem_err) 655 mem_err_s = mem_err_s + mem_err 656 endif 657 if (associated(amg_ilwork(level)%p)) then 658 deallocate(amg_ilwork(level)%p, 659 & stat=mem_err) 660 mem_err_s = mem_err_s + mem_err 661 write(*,*)'mcheck deallocate ilwork,',level,myrank 662 endif 663 664 if (associated(amg_paramap(level)%p)) then 665 deallocate(amg_paramap(level)%p, 666 & stat=mem_err) 667 mem_err_s = mem_err_s + mem_err 668 endif 669 if (associated(amg_paraext(level)%p)) then 670 deallocate(amg_paraext(level)%p, 671 & stat=mem_err) 672 mem_err_s = mem_err_s + mem_err 673 endif 674 675 if (mem_err_s .ne. 0 ) then 676 write(6,7002)level 677 end if 678 679 if (associated(I_fc_colm(level)%p)) then 680 deallocate(CF_map(level)%p,stat=mem_err) 681 mem_err_s = mem_err_s + mem_err 682 deallocate(CF_revmap(level)%p,stat=mem_err) 683 mem_err_s = mem_err_s + mem_err 684 deallocate(I_cf_colm(level)%p,stat=mem_err) 685 mem_err_s = mem_err_s + mem_err 686 deallocate(I_cf_rowp(level)%p,stat=mem_err) 687 mem_err_s = mem_err_s + mem_err 688 deallocate(I_cf(level)%p,stat=mem_err) 689 mem_err_s = mem_err_s + mem_err 690 deallocate(I_fc_colm(level)%p,stat=mem_err) 691 mem_err_s = mem_err_s + mem_err 692 deallocate(I_fc_rowp(level)%p,stat=mem_err) 693 mem_err_s = mem_err_s + mem_err 694 deallocate(I_fc(level)%p,stat=mem_err) 695 mem_err_s = mem_err_s + mem_err 696 end if 697 amg_nnz(level) = 0 698 amg_nshg(level) = 0 6997002 format(/' **** ramg: Deallocation error at level',i5) 700 701 return 702 703 end subroutine ! ramg_deallocate 704 705 subroutine ramg_readin_i(iarray,rn1,rfname) 706 include "common.h" 707 include "mpif.h" 708 include "auxmpi.h" 709 integer rn1 710 integer,dimension(rn1) :: iarray 711 character*10 rfname 712 character*5 cname 713 714 character*20 mfname 715 716 integer i,t 717 mfname = trim(rfname)//'.dat'//cname(myrank+1) 718 719 write(*,*)'RAMG READ: ',mfname 720 open(264,file=trim(mfname),status='unknown') 721 do i=1,rn1 722 read(264,*)t,iarray(i) 723 enddo 724 close(264) 725 end subroutine ! ramg_readin_i 726 727 728 subroutine ramg_dump(rarray,rn1,rfname) 729 include "common.h" 730 include "mpif.h" 731 include "auxmpi.h" 732 integer rn1 733 real(kind=8),dimension(rn1) :: rarray 734 character*10 rfname 735 character*5 cname 736 737 character*20 mfname 738 739 integer i 740 mfname = trim(rfname)//'.dat'//cname(myrank+1) 741 742 write(*,*)'RAMG DUMP: ',mfname 743 open(264,file=trim(mfname),status='unknown') 744 do i=1,rn1 745 !write(264,'((I10)(A)(D10.3))')i,' ',rarray(i) 746 write(264,*)i,rarray(i) 747 enddo 748 close(264) 749 750 end subroutine ! ramg_dump 751 752 subroutine ramg_dump_ic(rarray,rn1,rn2,rfname) 753 include "common.h" 754 include "mpif.h" 755 include "auxmpi.h" 756 integer rn1,rn2 757 integer,dimension(rn1,rn2) :: rarray 758 character*10 rfname 759 character*5 cname 760 761 character*20 mfname 762 763 integer i 764 mfname = trim(rfname)//'.dat'//cname(myrank+1) 765 766 write(*,*)'RAMG DUMP: ',mfname 767 open(264,file=trim(mfname),status='unknown') 768 write(264,*)rn2 769 do j=1,rn2 770 write(264,*)j,(rarray(i,j),i=1,rn1) 771 enddo 772 close(264) 773 end subroutine 774 775 subroutine ramg_dump_i(rarray,rn1,rn2,rfname) 776 include "common.h" 777 include "mpif.h" 778 include "auxmpi.h" 779 integer rn1,rn2 780 integer,dimension(rn1,rn2) :: rarray 781 character*10 rfname 782 character*5 cname 783 784 character*20 mfname 785 786 integer i 787 mfname = trim(rfname)//'.dat'//cname(myrank+1) 788 789 write(*,*)'RAMG DUMP: ',mfname 790 open(264,file=trim(mfname),status='unknown') 791 write(264,*)rn1 792 do i=1,rn1 793 write(264,*)i,(rarray(i,j),j=1,rn2) 794 enddo 795 close(264) 796 797 end subroutine ! ramg_dump_i 798 799 subroutine ramg_dump_ir(iarray,rarray,rn1,rn2,rfname) 800 include "common.h" 801 include "mpif.h" 802 include "auxmpi.h" 803 integer rn1,rn2 804 integer,dimension(rn1) :: iarray 805 real(kind=8),dimension(rn1,rn2) :: rarray 806 character*10 rfname 807 character*5 cname 808 809 character*20 mfname 810 811 character*20 mformat 812 813 integer i 814 mfname = trim(rfname)//'.dat'//cname(myrank+1) 815 816 write(mformat,'((A6)(I1)(A7))')'((2I)(',rn2,'E14.3))' 817 mformat = trim(mformat) 818 write(*,*)'RAMG DUMP: ',mfname 819 open(264,file=trim(mfname),status='unknown') 820 do i=1,rn1 821 write(264,mformat) 822 & i,iarray(i),(rarray(i,j),j=1,rn2) 823 enddo 824 close(264) 825 826 end subroutine ! ramg_dump_ir 827 828 subroutine ramg_dump_rn_map(rarray,rn1,rn2,rfname) 829 use ramg_data 830 include "common.h" 831 include "mpif.h" 832 include "auxmpi.h" 833 integer rn1,rn2 834 real(kind=8),dimension(rn1,rn2) :: rarray 835 character*10 rfname 836 character*5 cname 837 838 character*20 mfname 839 character*20 mformat 840 841 integer i,j,ii 842 mfname = trim(rfname)//'.dat'//cname(myrank+1) 843 844 write(*,*)'RAMG DUMP: ',mfname 845 open(264,file=trim(mfname),status='unknown') 846 write(mformat,'((A6)(I1)(A7))')'((1I)(',rn2,'E14.4))' 847 mformat=trim(mformat) 848 do i=1,rn1 849 if (numpe.gt.1) then 850 ii = ncorp_map(i) 851 else 852 ii = i 853 endif 854 write(264,mformat)ii,(rarray(i,j),j=1,rn2) 855 enddo 856 close(264) 857 end subroutine ! ramg_dump_rn_map 858 859 subroutine ramg_dump_rn(rarray,rn1,rn2,rfname) 860 include "common.h" 861 include "mpif.h" 862 include "auxmpi.h" 863 integer rn1,rn2 864 real(kind=8),dimension(rn1,rn2) :: rarray 865 character*10 rfname 866 character*5 cname 867 868 character*20 mfname 869 character*20 mformat 870 871 integer i,j 872 mfname = trim(rfname)//'.dat'//cname(myrank+1) 873 874 write(*,*)'RAMG DUMP: ',mfname 875 open(264,file=trim(mfname),status='unknown') 876 write(mformat,'((A6)(I1)(A7))')'((1I)(',rn2,'E14.4))' 877 mformat=trim(mformat) 878 do i=1,rn1 879 write(264,mformat)i,(rarray(i,j),j=1,rn2) 880 enddo 881 close(264) 882 883 end subroutine ! ramg_dump_rn 884 885 subroutine ramg_dump_matlab_A(acolm,arowp,alhs,an,annz,redun, 886 & fname) 887 888 use ramg_data 889 include "common.h" 890 include "mpif.h" 891 include "auxmpi.h" 892 893 integer :: an,annz,redun 894 integer,dimension(an+1) :: acolm 895 integer,dimension(annz) :: arowp 896 real(kind=8),dimension(redun,annz) :: alhs 897 character fname*10,mtxname*5 898 character cname*5 899 900 character mfname*15,mAname*5 901 character mformat*20 902 903 integer i,j,p,k 904 905 mfname = trim(fname)//'.dat'//cname(myrank+1) 906 907 write(*,*)'RAMG Dump to matlab ',mfname,myrank 908 open(264,file=trim(mfname),status='unknown') 909 !write(264,*)an 910 !write(264,*)annz 911 !write(264,*)'1' 912 write(mformat,'((A6)(I1)(A7))')'((2I)(',redun,'E14.4))' 913 do i=1,an 914 do p=acolm(i),acolm(i+1)-1 915 j = arowp(p) 916 !write(264,"((I6)(I7)(A)(E8.2))")i,j,' ',alhs(p) 917! if ( (amg_paramap(1)%p(i).eq.amg_paramap(2)%p(j)).and. 918! if (iabs(amg_paramap(1)%p(i)).ne.(myrank+1)) then 919! & (amg_paramap(1)%p(i).ne.1)) then 920 write(264,mformat)i,j,(alhs(k,p),k=1,redun) 921! endif 922 enddo 923 enddo 924 close(264) 925 end subroutine 926 927 subroutine ramg_dump_matlab_map(acolm,arowp,alhs,an,annz,redun, 928 & fname) 929 use ramg_data 930 include "common.h" 931 include "mpif.h" 932 include "auxmpi.h" 933 934 integer :: an,annz,redun 935 integer,dimension(an+1) :: acolm 936 integer,dimension(annz) :: arowp 937 real(kind=8),dimension(redun,annz) :: alhs 938 !integer,dimension(an) :: ipmap 939 character fname*10 940 941 character mfname*20 942 character cname*5 943 character mformat*20 944 945 integer i,j,p,k 946 integer ii,jj 947 948 if (numpe.eq.1) then 949 call ramg_dump_matlab_A(acolm,arowp,alhs,an,annz,redun, 950 & fname) 951 return; 952 endif 953 954 mfname = trim(fname)//'.dat'//cname(myrank+1) 955 956 write(*,*)'RAMG Dump to matlab ',mfname,nshg,myrank 957 open(264,file=trim(mfname),status='unknown') 958 !write(264,*)an 959 !write(264,*)annz 960 !write(264,*)'1' 961 write(mformat,'((A6)(I1)(A7))')'((2I)(',redun,'E14.4))' 962 !call ramg_dump_i(ncorp_map,nshg,1,'pam_corp ') 963 do i=1,an 964 ii = ncorp_map(i)!ipmap(i)) 965 do p=acolm(i),acolm(i+1)-1 966 j = arowp(p) 967 jj = ncorp_map(j)!ipmap(j)) 968 !write(264,"((I6)(I7)(A)(E8.2))")i,j,' ',alhs(p) 969 write(264,mformat)ii,jj,(alhs(k,p),k=1,redun) 970 enddo 971 enddo 972 close(264) 973 end subroutine 974 975 976 subroutine ramg_dump_mupad_A(acolm,arowp,alhs,an,annz, 977 & fname,mtxname) 978 integer :: an,annz 979 integer,dimension(an+1) :: acolm 980 integer,dimension(annz) :: arowp 981 real(kind=8),dimension(annz) :: alhs 982 character fname*10,mtxname*5 983 984 character mfname*15,mAname*5 985 986 integer i,j,p,k 987 988 mfname = trim(fname)//'.mws' 989 mAname = trim(mtxname) 990 991 write(*,*)'RAMG Dump to Mupad ',mfname 992 open(264,file=trim(mfname),status='unknown') 993 write(264,*)mAname,':= matrix(',an,',',an,'):' 994 do i=1,an 995 do p=acolm(i),acolm(i+1)-1 996 j = arowp(p) 997 write(264,*)mAname,'[',i,',',j,']:=',alhs(p),':' 998 enddo 999 enddo 1000 close(264) 1001 1002 end subroutine 1003 1004 subroutine ramg_print_time(str,v) 1005 1006 include "common.h" 1007 include "mpif.h" 1008 include "auxmpi.h" 1009 1010 character*(*) str 1011 real(kind=8) :: v,vave,vmax 1012 1013 if (numpe.gt.1) then 1014 call MPI_AllReduce(v,vmax,1, 1015 & MPI_DOUBLE_PRECISION,MPI_MAX,MPI_COMM_WORLD,ierr) 1016 call MPI_AllReduce(v,vave,1, 1017 & MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr) 1018 vave = vave/numpe 1019 else 1020 vave = v 1021 vmax = v 1022 endif 1023 1024 if ((iamg_verb.gt.1).and.(myrank.eq.master)) then 1025 write(*,7101)trim(str),vave,vmax 10267101 format(T1,A,T40,f9.2,' s (ave), ',f9.2,' s (max)') 1027 endif 1028 1029 end subroutine 1030 1031 subroutine ramg_output_coarsening 1032 use readarrays 1033 use ramg_data 1034 include "common.h" 1035 include "mpif.h" 1036 1037 character*20 mfname 1038 character*20 mformat 1039 1040 integer i 1041 1042 write(*,*)'ramg dump coarsening' 1043 mfname = 'amgcoarsen.dat' 1044 open(265,file=trim(mfname),status='unknown') 1045 write(265,*)nshg 1046 do i=1,nshg 1047 write(265,'((2I)(3E14.5))')i,amg_cfmap(i), 1048 & point2x(i,1),point2x(i,2),point2x(i,3) 1049 enddo 1050 close(265) 1051 1052 end subroutine 1053 1054!*********************************************************** 1055! <EOF> ramg TOOLS 1056!********************************************************** 1057