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,fname1,map_nshg,ione,"integer",iotype) 426 427 call readdatablock(igeom,fname1,ncorp_map,map_nshg, 428 & "integer",iotype) 429 430 call closefile(igeom,"read") 431 432 endif 433 434 end subroutine ! ramg_read_map 435 436!*********************************************** 437! ramg_check_diag 438!*********************************************** 439 subroutine ramg_check_diag(colm,rowp,lhs,nshg,nnz_tot) 440 implicit none 441 integer,intent(in) :: nshg,nnz_tot 442 integer,intent(in),dimension(nshg+1) :: colm 443 integer,intent(in),dimension(nnz_tot) :: rowp 444 real(kind=8),intent(in),dimension(nnz_tot) :: lhs 445 446 integer :: i,j,p 447 real(kind=8) :: diagrow 448 logical :: diagokay 449 450 diagokay = .true. 451 452 do i=1,nshg 453 p = colm(i) 454 if (rowp(p).ne.i) then 455 write(*,*)'matrix first row entry is not diagonal',i 456 diagokay = .false. 457 end if 458 diagrow = lhs(p) 459 if (diagrow.lt.0) then 460 write(*,*)'matrix diagonal < 0 at',i,diagrow 461 diagokay = .false. 462 end if 463 do j = colm(i)+1,colm(i+1)-1 464 p = rowp(j) 465 if (lhs(p).gt.diagrow) then 466 write(*,*)'matrix not diagonal dominant at row',i 467 diagokay = .false. 468 write(*,*)'diag:',diagrow,p,lhs(p) 469 exit 470 end if 471 end do 472 enddo 473 474 if (diagokay) then 475 write(*,*)'matrix check diagonal dominance okay' 476 end if 477 478 end subroutine 479 480 subroutine ramg_dot_p(level,v,u,redun,norm) 481 use ramg_data 482 include "common.h" 483 include "mpif.h" 484 include "auxmpi.h" 485 486 integer :: redun,level 487 real(kind=8),intent(in),dimension(amg_nshg(level),redun) ::v,u 488 real(kind=8),intent(inout) :: norm 489 integer :: i,k 490 real(kind=8) :: normt 491 normt = 0 492 do i=1,amg_nshg(level) 493 do k=1,redun 494 normt = normt + v(i,k)*u(i,k) 495 enddo 496 enddo 497 if (numpe.gt.1) then 498 call MPI_AllReduce(normt,norm,1,MPI_DOUBLE_PRECISION,MPI_SUM, 499 & MPI_COMM_WORLD,ierr) 500 else 501 norm=normt 502 endif 503 end subroutine !ramg_dot_p 504 505!************************************************************ 506! ramg_L2_norm 507! calculate norm = | r | 508!************************************************************ 509 subroutine ramg_L2_norm(nshg,v,norm) 510 implicit none 511 integer,intent(in) :: nshg 512 real(kind=8),intent(in),dimension(nshg) :: v 513 real(kind=8),intent(inout) :: norm 514 integer :: i 515 norm = 0 516 do i=1,nshg 517 norm = norm + v(i)*v(i) 518 enddo 519 norm = sqrt(norm) 520 end subroutine !ramg_L2_norm 521 522 subroutine ramg_L2_norm_p(level,v,vflow,norm) 523 use ramg_data 524 include "common.h" 525 include "mpif.h" 526 include "auxmpi.h" 527 528 integer,intent(in) :: vflow,level 529 real(kind=8),intent(in),dimension(amg_nshg(level),vflow) :: v 530 real(kind=8),intent(inout) :: norm 531 integer :: i,k 532 real(kind=8) :: normt 533 normt = 0 534 do i=1,amg_nshg(level) 535 do k=1,vflow 536 normt = normt + v(i,k)*v(i,k) 537 enddo 538 enddo 539 if (numpe.gt.1) then 540 call MPI_AllReduce(normt,norm,1,MPI_DOUBLE_PRECISION,MPI_SUM, 541 & MPI_COMM_WORLD,ierr) 542 else 543 norm = normt 544 endif 545 norm = sqrt(norm) 546 end subroutine !ramg_L2_norm_p 547 548!************************************************************ 549! ramg_allocate & deallocate 550! (de)allocate memory for level l, lhs matrix, rhs vector 551!************************************************************ 552 subroutine ramg_allocate( 553 & level,lnshg,lnnz_tot,n_sol) 554 use ramg_data 555 implicit none 556 557 integer,intent(in) :: level,lnshg,lnnz_tot 558 integer,intent(in) :: n_sol 559 integer :: mem_err,mem_err_s 560 mem_err_s = 0 561 if (lnnz_tot.ne.0) then ! zero means manual alloc later 562 amg_nnz(level) = lnnz_tot 563 allocate(amg_A_lhs(level)%p(amg_nnz(level),n_sol), 564 & stat=mem_err) 565 mem_err_s = mem_err_s + mem_err 566 allocate(amg_A_rowp(level)%p(amg_nnz(level)), 567 & stat=mem_err) 568 mem_err_s = mem_err_s + mem_err 569 endif 570 if (lnshg.ne.0) then 571 amg_nshg(level) = lnshg 572 allocate(amg_A_colm(level)%p(amg_nshg(level)+1), 573 & stat=mem_err) 574 mem_err_s = mem_err_s + mem_err 575 allocate(amg_A_rhs(level)%p(amg_nshg(level)), 576 & stat=mem_err) 577 mem_err_s = mem_err_s + mem_err 578 allocate(amg_ppeDiag(level)%p(amg_nshg(level)), 579 & stat=mem_err) 580 endif 581 mem_err_s = mem_err_s + mem_err 582 if (mem_err_s .ne. 0 ) then 583 write(6,7001)level 584 end if 5857001 format(/' **** ramg: Allocation error at level',i5) 586 ! zero out 587 if (lnnz_tot.ne.0) then 588 amg_A_lhs(level)%p(:,:) = 0 589 amg_A_rowp(level)%p(:) = 0 590 endif 591 if (lnshg.ne.0) then 592 amg_A_colm(level)%p(:) = 0 593 amg_A_rhs(level)%p(:) = 0 594 endif 595 return 596 597 end subroutine ! ramg_allocate 598 599 subroutine ramg_deallocate(level) 600 601 use ramg_data 602 include "common.h" 603 604 integer,intent(inout) :: level 605 integer :: mem_err,mem_err_s,i 606 607 608 if (level.eq.1) then 609 610 if (maxnev.gt.0) then 611 deallocate(ramg_ggb_ev) 612 deallocate(ramg_ggb_eA) 613 deallocate(cmtxA) 614 deallocate(cindx) 615 endif 616 617 618 if (iamg_reduce.gt.1) then 619 deallocate(reducemap) 620 deallocate(rmap1d) 621 endif 622 623 end if 624 625 626 level = abs(level) 627 628 mem_err_s = 0 629 if (associated(amg_A_lhs(level)%p)) then 630 deallocate(amg_A_lhs(level)%p, 631 & stat=mem_err) 632 mem_err_s = mem_err_s + mem_err 633 endif 634 if (associated(amg_A_rowp(level)%p)) then 635 deallocate(amg_A_rowp(level)%p, 636 & stat=mem_err) 637 mem_err_s = mem_err_s + mem_err 638 endif 639 if (associated(amg_A_colm(level)%p)) then 640 deallocate(amg_A_colm(level)%p, 641 & stat=mem_err) 642 mem_err_s = mem_err_s + mem_err 643 endif 644 if (associated(amg_ppeDiag(level)%p)) then 645 deallocate(amg_ppeDiag(level)%p, 646 & stat=mem_err) 647 mem_err_s = mem_err_s + mem_err 648 endif 649 if (associated(amg_A_rhs(level)%p)) then 650 deallocate(amg_A_rhs(level)%p, 651 & stat=mem_err) 652 mem_err_s = mem_err_s + mem_err 653 endif 654 if (associated(amg_ilwork(level)%p)) then 655 deallocate(amg_ilwork(level)%p, 656 & stat=mem_err) 657 mem_err_s = mem_err_s + mem_err 658 write(*,*)'mcheck deallocate ilwork,',level,myrank 659 endif 660 661 if (associated(amg_paramap(level)%p)) then 662 deallocate(amg_paramap(level)%p, 663 & stat=mem_err) 664 mem_err_s = mem_err_s + mem_err 665 endif 666 if (associated(amg_paraext(level)%p)) then 667 deallocate(amg_paraext(level)%p, 668 & stat=mem_err) 669 mem_err_s = mem_err_s + mem_err 670 endif 671 672 if (mem_err_s .ne. 0 ) then 673 write(6,7002)level 674 end if 675 676 if (associated(I_fc_colm(level)%p)) then 677 deallocate(CF_map(level)%p,stat=mem_err) 678 mem_err_s = mem_err_s + mem_err 679 deallocate(CF_revmap(level)%p,stat=mem_err) 680 mem_err_s = mem_err_s + mem_err 681 deallocate(I_cf_colm(level)%p,stat=mem_err) 682 mem_err_s = mem_err_s + mem_err 683 deallocate(I_cf_rowp(level)%p,stat=mem_err) 684 mem_err_s = mem_err_s + mem_err 685 deallocate(I_cf(level)%p,stat=mem_err) 686 mem_err_s = mem_err_s + mem_err 687 deallocate(I_fc_colm(level)%p,stat=mem_err) 688 mem_err_s = mem_err_s + mem_err 689 deallocate(I_fc_rowp(level)%p,stat=mem_err) 690 mem_err_s = mem_err_s + mem_err 691 deallocate(I_fc(level)%p,stat=mem_err) 692 mem_err_s = mem_err_s + mem_err 693 end if 694 amg_nnz(level) = 0 695 amg_nshg(level) = 0 6967002 format(/' **** ramg: Deallocation error at level',i5) 697 698 return 699 700 end subroutine ! ramg_deallocate 701 702 subroutine ramg_readin_i(iarray,rn1,rfname) 703 include "common.h" 704 include "mpif.h" 705 include "auxmpi.h" 706 integer rn1 707 integer,dimension(rn1) :: iarray 708 character*10 rfname 709 character*5 cname 710 711 character*20 mfname 712 713 integer i,t 714 mfname = trim(rfname)//'.dat'//cname(myrank+1) 715 716 write(*,*)'RAMG READ: ',mfname 717 open(264,file=trim(mfname),status='unknown') 718 do i=1,rn1 719 read(264,*)t,iarray(i) 720 enddo 721 close(264) 722 end subroutine ! ramg_readin_i 723 724 725 subroutine ramg_dump(rarray,rn1,rfname) 726 include "common.h" 727 include "mpif.h" 728 include "auxmpi.h" 729 integer rn1 730 real(kind=8),dimension(rn1) :: rarray 731 character*10 rfname 732 character*5 cname 733 734 character*20 mfname 735 736 integer i 737 mfname = trim(rfname)//'.dat'//cname(myrank+1) 738 739 write(*,*)'RAMG DUMP: ',mfname 740 open(264,file=trim(mfname),status='unknown') 741 do i=1,rn1 742 !write(264,'((I10)(A)(D10.3))')i,' ',rarray(i) 743 write(264,*)i,rarray(i) 744 enddo 745 close(264) 746 747 end subroutine ! ramg_dump 748 749 subroutine ramg_dump_ic(rarray,rn1,rn2,rfname) 750 include "common.h" 751 include "mpif.h" 752 include "auxmpi.h" 753 integer rn1,rn2 754 integer,dimension(rn1,rn2) :: rarray 755 character*10 rfname 756 character*5 cname 757 758 character*20 mfname 759 760 integer i 761 mfname = trim(rfname)//'.dat'//cname(myrank+1) 762 763 write(*,*)'RAMG DUMP: ',mfname 764 open(264,file=trim(mfname),status='unknown') 765 write(264,*)rn2 766 do j=1,rn2 767 write(264,*)j,(rarray(i,j),i=1,rn1) 768 enddo 769 close(264) 770 end subroutine 771 772 subroutine ramg_dump_i(rarray,rn1,rn2,rfname) 773 include "common.h" 774 include "mpif.h" 775 include "auxmpi.h" 776 integer rn1,rn2 777 integer,dimension(rn1,rn2) :: rarray 778 character*10 rfname 779 character*5 cname 780 781 character*20 mfname 782 783 integer i 784 mfname = trim(rfname)//'.dat'//cname(myrank+1) 785 786 write(*,*)'RAMG DUMP: ',mfname 787 open(264,file=trim(mfname),status='unknown') 788 write(264,*)rn1 789 do i=1,rn1 790 write(264,*)i,(rarray(i,j),j=1,rn2) 791 enddo 792 close(264) 793 794 end subroutine ! ramg_dump_i 795 796 subroutine ramg_dump_ir(iarray,rarray,rn1,rn2,rfname) 797 include "common.h" 798 include "mpif.h" 799 include "auxmpi.h" 800 integer rn1,rn2 801 integer,dimension(rn1) :: iarray 802 real(kind=8),dimension(rn1,rn2) :: rarray 803 character*10 rfname 804 character*5 cname 805 806 character*20 mfname 807 808 character*20 mformat 809 810 integer i 811 mfname = trim(rfname)//'.dat'//cname(myrank+1) 812 813 write(mformat,'((A6)(I1)(A7))')'((2I)(',rn2,'E14.3))' 814 mformat = trim(mformat) 815 write(*,*)'RAMG DUMP: ',mfname 816 open(264,file=trim(mfname),status='unknown') 817 do i=1,rn1 818 write(264,mformat) 819 & i,iarray(i),(rarray(i,j),j=1,rn2) 820 enddo 821 close(264) 822 823 end subroutine ! ramg_dump_ir 824 825 subroutine ramg_dump_rn_map(rarray,rn1,rn2,rfname) 826 use ramg_data 827 include "common.h" 828 include "mpif.h" 829 include "auxmpi.h" 830 integer rn1,rn2 831 real(kind=8),dimension(rn1,rn2) :: rarray 832 character*10 rfname 833 character*5 cname 834 835 character*20 mfname 836 character*20 mformat 837 838 integer i,j,ii 839 mfname = trim(rfname)//'.dat'//cname(myrank+1) 840 841 write(*,*)'RAMG DUMP: ',mfname 842 open(264,file=trim(mfname),status='unknown') 843 write(mformat,'((A6)(I1)(A7))')'((1I)(',rn2,'E14.4))' 844 mformat=trim(mformat) 845 do i=1,rn1 846 if (numpe.gt.1) then 847 ii = ncorp_map(i) 848 else 849 ii = i 850 endif 851 write(264,mformat)ii,(rarray(i,j),j=1,rn2) 852 enddo 853 close(264) 854 end subroutine ! ramg_dump_rn_map 855 856 subroutine ramg_dump_rn(rarray,rn1,rn2,rfname) 857 include "common.h" 858 include "mpif.h" 859 include "auxmpi.h" 860 integer rn1,rn2 861 real(kind=8),dimension(rn1,rn2) :: rarray 862 character*10 rfname 863 character*5 cname 864 865 character*20 mfname 866 character*20 mformat 867 868 integer i,j 869 mfname = trim(rfname)//'.dat'//cname(myrank+1) 870 871 write(*,*)'RAMG DUMP: ',mfname 872 open(264,file=trim(mfname),status='unknown') 873 write(mformat,'((A6)(I1)(A7))')'((1I)(',rn2,'E14.4))' 874 mformat=trim(mformat) 875 do i=1,rn1 876 write(264,mformat)i,(rarray(i,j),j=1,rn2) 877 enddo 878 close(264) 879 880 end subroutine ! ramg_dump_rn 881 882 subroutine ramg_dump_matlab_A(acolm,arowp,alhs,an,annz,redun, 883 & fname) 884 885 use ramg_data 886 include "common.h" 887 include "mpif.h" 888 include "auxmpi.h" 889 890 integer :: an,annz,redun 891 integer,dimension(an+1) :: acolm 892 integer,dimension(annz) :: arowp 893 real(kind=8),dimension(redun,annz) :: alhs 894 character fname*10,mtxname*5 895 character cname*5 896 897 character mfname*15,mAname*5 898 character mformat*20 899 900 integer i,j,p,k 901 902 mfname = trim(fname)//'.dat'//cname(myrank+1) 903 904 write(*,*)'RAMG Dump to matlab ',mfname,myrank 905 open(264,file=trim(mfname),status='unknown') 906 !write(264,*)an 907 !write(264,*)annz 908 !write(264,*)'1' 909 write(mformat,'((A6)(I1)(A7))')'((2I)(',redun,'E14.4))' 910 do i=1,an 911 do p=acolm(i),acolm(i+1)-1 912 j = arowp(p) 913 !write(264,"((I6)(I7)(A)(E8.2))")i,j,' ',alhs(p) 914! if ( (amg_paramap(1)%p(i).eq.amg_paramap(2)%p(j)).and. 915! if (iabs(amg_paramap(1)%p(i)).ne.(myrank+1)) then 916! & (amg_paramap(1)%p(i).ne.1)) then 917 write(264,mformat)i,j,(alhs(k,p),k=1,redun) 918! endif 919 enddo 920 enddo 921 close(264) 922 end subroutine 923 924 subroutine ramg_dump_matlab_map(acolm,arowp,alhs,an,annz,redun, 925 & fname) 926 use ramg_data 927 include "common.h" 928 include "mpif.h" 929 include "auxmpi.h" 930 931 integer :: an,annz,redun 932 integer,dimension(an+1) :: acolm 933 integer,dimension(annz) :: arowp 934 real(kind=8),dimension(redun,annz) :: alhs 935 !integer,dimension(an) :: ipmap 936 character fname*10 937 938 character mfname*20 939 character cname*5 940 character mformat*20 941 942 integer i,j,p,k 943 integer ii,jj 944 945 if (numpe.eq.1) then 946 call ramg_dump_matlab_A(acolm,arowp,alhs,an,annz,redun, 947 & fname) 948 return; 949 endif 950 951 mfname = trim(fname)//'.dat'//cname(myrank+1) 952 953 write(*,*)'RAMG Dump to matlab ',mfname,nshg,myrank 954 open(264,file=trim(mfname),status='unknown') 955 !write(264,*)an 956 !write(264,*)annz 957 !write(264,*)'1' 958 write(mformat,'((A6)(I1)(A7))')'((2I)(',redun,'E14.4))' 959 !call ramg_dump_i(ncorp_map,nshg,1,'pam_corp ') 960 do i=1,an 961 ii = ncorp_map(i)!ipmap(i)) 962 do p=acolm(i),acolm(i+1)-1 963 j = arowp(p) 964 jj = ncorp_map(j)!ipmap(j)) 965 !write(264,"((I6)(I7)(A)(E8.2))")i,j,' ',alhs(p) 966 write(264,mformat)ii,jj,(alhs(k,p),k=1,redun) 967 enddo 968 enddo 969 close(264) 970 end subroutine 971 972 973 subroutine ramg_dump_mupad_A(acolm,arowp,alhs,an,annz, 974 & fname,mtxname) 975 integer :: an,annz 976 integer,dimension(an+1) :: acolm 977 integer,dimension(annz) :: arowp 978 real(kind=8),dimension(annz) :: alhs 979 character fname*10,mtxname*5 980 981 character mfname*15,mAname*5 982 983 integer i,j,p,k 984 985 mfname = trim(fname)//'.mws' 986 mAname = trim(mtxname) 987 988 write(*,*)'RAMG Dump to Mupad ',mfname 989 open(264,file=trim(mfname),status='unknown') 990 write(264,*)mAname,':= matrix(',an,',',an,'):' 991 do i=1,an 992 do p=acolm(i),acolm(i+1)-1 993 j = arowp(p) 994 write(264,*)mAname,'[',i,',',j,']:=',alhs(p),':' 995 enddo 996 enddo 997 close(264) 998 999 end subroutine 1000 1001 subroutine ramg_print_time(str,v) 1002 1003 include "common.h" 1004 include "mpif.h" 1005 include "auxmpi.h" 1006 1007 character*(*) str 1008 real(kind=8) :: v,vave,vmax 1009 1010 if (numpe.gt.1) then 1011 call MPI_AllReduce(v,vmax,1, 1012 & MPI_DOUBLE_PRECISION,MPI_MAX,MPI_COMM_WORLD,ierr) 1013 call MPI_AllReduce(v,vave,1, 1014 & MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr) 1015 vave = vave/numpe 1016 else 1017 vave = v 1018 vmax = v 1019 endif 1020 1021 if ((iamg_verb.gt.1).and.(myrank.eq.master)) then 1022 write(*,7101)trim(str),vave,vmax 10237101 format(T1,A,T40,f9.2,' s (ave), ',f9.2,' s (max)') 1024 endif 1025 1026 end subroutine 1027 1028 subroutine ramg_output_coarsening 1029 use readarrays 1030 use ramg_data 1031 include "common.h" 1032 include "mpif.h" 1033 1034 character*20 mfname 1035 character*20 mformat 1036 1037 integer i 1038 1039 write(*,*)'ramg dump coarsening' 1040 mfname = 'amgcoarsen.dat' 1041 open(265,file=trim(mfname),status='unknown') 1042 write(265,*)nshg 1043 do i=1,nshg 1044 write(265,'((2I)(3E14.5))')i,amg_cfmap(i), 1045 & point2x(i,1),point2x(i,2),point2x(i,3) 1046 enddo 1047 close(265) 1048 1049 end subroutine 1050 1051!*********************************************************** 1052! <EOF> ramg TOOLS 1053!********************************************************** 1054