1!C Copyright (C) 2005 The Scalable Software Infrastructure Project. All rights reserved. 2!C 3!C Redistribution and use in source and binary forms, with or without 4!C modification, are permitted provided that the following conditions are met: 5!C 1. Redistributions of source code must retain the above copyright 6!C notice, this list of conditions and the following disclaimer. 7!C 2. Redistributions in binary form must reproduce the above copyright 8!C notice, this list of conditions and the following disclaimer in the 9!C documentation and/or other materials provided with the distribution. 10!C 3. Neither the name of the project nor the names of its contributors 11!C may be used to endorse or promote products derived from this software 12!C without specific prior written permission. 13!C 14!C THIS SOFTWARE IS PROVIDED BY THE SCALABLE SOFTWARE INFRASTRUCTURE PROJECT 15!C ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 16!C TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 17!C PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE SCALABLE SOFTWARE INFRASTRUCTURE 18!C PROJECT BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, 19!C OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 20!C SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 21!C INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 22!C CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 23!C ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 24!C POSSIBILITY OF SUCH DAMAGE. 25 26#ifndef ZERO_ORIGIN 27#define ZERO_ORIGIN 0 28#endif 29 30!C ************************************************ 31!C * MODULE solver_AMGCG 32!C CONTAINS 33!C * SUBROUTINE v_cycle 34!C * SUBROUTINE sgs 35!C * SUBROUTINE matrix_counting 36!C * SUBROUTINE matrix_arrange 37!C * SUBROUTINE AMGCG 38!C * SUBROUTINE clear_matrix 39!C ************************************************ 40 41MODULE solver_AMGCG 42CONTAINS 43 44 SUBROUTINE v_cycle(N, problem_B, problem_X, LEVEL_NUM, Temp) 45 USE data_structure_for_AMG 46 IMPLICIT NONE 47 48 INTEGER(kind=kint), INTENT(in) :: LEVEL_NUM 49 INTEGER(kind=kint) :: N, NP 50 INTEGER(kind=kint) :: TMP_INT 51 REAL (kind=kreal), DIMENSION(:), pointer :: D 52 !C DIMENSION(N) 53 REAL (kind=kreal), DIMENSION(N), target :: problem_B 54 !C DIMENSION(N) 55 REAL (kind=kreal), DIMENSION(N), target:: problem_X 56 !C DIMENSION(N) 57 58 REAL (kind=kreal), DIMENSION(: ), pointer:: B 59 !C DIMENSION(N) 60 REAL (kind=kreal), DIMENSION(: ), pointer:: X 61 !C DIMENSION(N) 62 63 REAL (kind=kreal), DIMENSION(:), pointer:: AU 64 REAL (kind=kreal), DIMENSION(:), pointer:: AL 65 66 INTEGER(kind=kint ), DIMENSION(:), pointer :: INU 67 !C DIMENSION(0:N) 68 INTEGER(kind=kint ), DIMENSION(:), pointer :: IAU 69 !C DIMENSION(NPU) 70 INTEGER(kind=kint ), DIMENSION(:), pointer :: INL 71 !C DIMENSION(0:N) 72 INTEGER(kind=kint ), DIMENSION(:), pointer :: IAL 73 74 REAL (kind=kreal), DIMENSION(:), pointer :: coarser_B 75 REAL (kind=kreal), DIMENSION(:), pointer :: coarser_X 76 REAL (kind=kreal), DIMENSION(N) :: Temp 77 78 INTEGER(kind=kint ) :: nth_lev,i,inod,j,ieL,isL,is,ie,coarser_N,isU,ieU,k,coarser_NP 79 INTEGER(kind=kint ) :: coarser_NEIBPETOT 80 REAL (kind=kreal) :: R_val,B_val,w,X_val,R_norm,GR_norm 81 82 83 type(INTER_LEVEL_OPERATOR),pointer :: R, P 84 85 86 87 HIERARCHICAL_DATA(1) % B =>problem_B 88 HIERARCHICAL_DATA(1) % X =>problem_X 89 90 do i = 1, N 91 problem_X(i) = 0.0D0 92 END do 93 94 DO nth_lev=1,LEVEL_NUM-1 95 N= HIERARCHICAL_DATA(nth_lev) % N 96 NP=HIERARCHICAL_DATA(nth_lev) % NP 97 INU=>HIERARCHICAL_DATA(nth_lev) % INU 98 INL=>HIERARCHICAL_DATA(nth_lev) % INL 99 X => HIERARCHICAL_DATA(nth_lev) % X 100 B => HIERARCHICAL_DATA(nth_lev) % B 101 D => HIERARCHICAL_DATA(nth_lev) % D 102 IAU=>HIERARCHICAL_DATA(nth_lev) % IAU 103 IAL=>HIERARCHICAL_DATA(nth_lev) % IAL 104 AU =>HIERARCHICAL_DATA(nth_lev) % AU 105 AL =>HIERARCHICAL_DATA(nth_lev) % AL 106 107 coarser_B => HIERARCHICAL_DATA(nth_lev+1) % B 108 coarser_X => HIERARCHICAL_DATA(nth_lev+1) % X 109 110 TMP_INT=1 111 CALL sgs(N,NP,D,AL,INL,IAL,AU,INU,IAU,B,X,TMP_INT,nth_lev) 112 113 DO j= 1, N 114 R_val= B(j) - D(j)*X(j) 115 isU= INU(j-1)+1 116 ieU= INU(j) 117 DO i= isU, ieU 118 inod=IAU(i)+ZERO_ORIGIN 119 R_val= R_val - AU(i) * X(inod) 120 END DO 121 isL= INL(j-1)+1 122 ieL= INL(j) 123 DO i= isL, ieL 124 inod= IAL(i)+ZERO_ORIGIN 125 R_val= R_val - AL(i) * X(inod) 126 END DO 127 Temp(j)=R_val 128 END DO 129 130 !C restrict the residual vector 131 R => HIERARCHICAL_DATA(nth_lev+1) % R 132 coarser_N= HIERARCHICAL_DATA(nth_lev+1) % N 133 coarser_NP=HIERARCHICAL_DATA(nth_lev+1) % NP 134 135 DO j= 1, coarser_NP 136 is= R % IN(j-1)+1 137 ie= R % IN(j) 138 B_val=0.0 139 DO i=is, ie 140 inod = R % CN(i) 141 B_val= B_val + R % V(i) * Temp(inod) 142 END DO 143 coarser_B(j)= B_val 144 END DO 145 146 DO j=1,coarser_NP 147 coarser_X(j)= 0.0 148 END DO 149 END DO 150 151 N=HIERARCHICAL_DATA(LEVEL_NUM) % N 152 NP=HIERARCHICAL_DATA(LEVEL_NUM) % NP 153 154 X => HIERARCHICAL_DATA(LEVEL_NUM) % X 155 B => HIERARCHICAL_DATA(LEVEL_NUM) % B 156 AL=> HIERARCHICAL_DATA(LEVEL_NUM) % AL 157 INL=>HIERARCHICAL_DATA(LEVEL_NUM) % INL 158 159 INU => HIERARCHICAL_DATA(LEVEL_NUM) % INU 160 D => HIERARCHICAL_DATA(LEVEL_NUM) % D 161 IAU => HIERARCHICAL_DATA(LEVEL_NUM) % IAU 162 IAL => HIERARCHICAL_DATA(LEVEL_NUM) % IAL 163 AU => HIERARCHICAL_DATA(LEVEL_NUM) % AU 164 TMP_INT=30 165 CALL sgs(N,NP,D,AL,INL,IAL,AU,INU,IAU,B,X,TMP_INT,LEVEL_NUM) 166 167 168 DO nth_lev=LEVEL_NUM-1,1,-1 169 N= HIERARCHICAL_DATA(nth_lev) % N 170 NP=HIERARCHICAL_DATA(nth_lev) % NP 171 172 INU => HIERARCHICAL_DATA(nth_lev) % INU 173 INL => HIERARCHICAL_DATA(nth_lev) % INL 174 X => HIERARCHICAL_DATA(nth_lev) % X 175 B => HIERARCHICAL_DATA(nth_lev) % B 176 D => HIERARCHICAL_DATA(nth_lev) % D 177 IAU => HIERARCHICAL_DATA(nth_lev) % IAU 178 IAL => HIERARCHICAL_DATA(nth_lev) % IAL 179 AU => HIERARCHICAL_DATA(nth_lev) % AU 180 AL => HIERARCHICAL_DATA(nth_lev) % AL 181 coarser_X => HIERARCHICAL_DATA(nth_lev+1) % X 182 P => HIERARCHICAL_DATA(nth_lev+1) % P 183 184 185 DO j= 1, N 186 is= P % IN(j-1)+1 187 ie= P % IN(j) 188 X_val=0.0 189 DO i=is, ie 190 inod = P % CN(i) 191 X_val= X_val+P % V(i) * coarser_X(inod) 192 END DO 193 X(j)=X(j)+X_val 194 END DO 195 196 TMP_INT=1 197 CALL sgs(N,NP,D,AL,INL,IAL,AU,INU,IAU,B,X,TMP_INT,nth_lev) 198 END DO 199 END SUBROUTINE v_cycle 200 201 202 SUBROUTINE sgs(N, NP, D, AL, INL, IAL, AU, INU, IAU, B, X, count, LEVEL_NO) 203 204 USE data_structure_for_AMG 205 206 207 IMPLICIT NONE 208 209 INTEGER(kind=kint),INTENT(in ) :: N,NP,count 210 211 REAL (kind=kreal), DIMENSION(:) :: D 212 REAL (kind=kreal), DIMENSION(:) :: B 213 REAL (kind=kreal), DIMENSION(:) :: X 214 215 REAL (kind=kreal), DIMENSION(:) :: AU 216 REAL (kind=kreal), DIMENSION(:) :: AL 217 218 INTEGER(kind=kint ), DIMENSION(0:) :: INU 219 INTEGER(kind=kint ), DIMENSION(:) :: IAU 220 INTEGER(kind=kint ), DIMENSION(0:) :: INL 221 INTEGER(kind=kint ), DIMENSION(:) :: IAL 222 223 INTEGER(kind=kint), intent(in) :: LEVEL_NO 224 225!!$ INTEGER(kind=kint), DIMENSION(:),pointer :: NEIBPE 226!!$ INTEGER(kind=kint), DIMENSION(:),pointer :: STACK_IMPORT,STACK_EXPORT 227!!$ INTEGER(kind=kint), DIMENSION(:),pointer :: NOD_EXPORT,NOD_IMPORT 228!!$ INTEGER(kind=kint) :: NEIBPETOT 229!!$ 230 231 232 233 INTEGER(kind=kint) :: cnt,isL,ieL,isU,ieU,i,j,inod 234 REAL (kind=kreal) :: R_val,R_norm,GR_norm,v 235 236!!$ NEIBPETOT = HIERARCHICAL_DATA(LEVEL_NO) % COMM_TABLE % NEIBPETOT 237!!$ NEIBPE => HIERARCHICAL_DATA(LEVEL_NO) % COMM_TABLE % NEIBPE 238!!$ STACK_IMPORT => HIERARCHICAL_DATA(LEVEL_NO) % COMM_TABLE % STACK_IMPORT 239!!$ STACK_EXPORT => HIERARCHICAL_DATA(LEVEL_NO) % COMM_TABLE % STACK_EXPORT 240!!$ NOD_IMPORT => HIERARCHICAL_DATA(LEVEL_NO) % COMM_TABLE % NOD_IMPORT 241!!$ NOD_EXPORT => HIERARCHICAL_DATA(LEVEL_NO) % COMM_TABLE % NOD_EXPORT 242 243 244 DO cnt=1,count 245 246 DO i=1,N 247 v = B(i) 248 249 isL=INL(i-1)+1 250 ieL=INL(i) 251 DO j=isL,ieL 252 inod=IAL(j)+ZERO_ORIGIN 253 v = v - AL(j)*X(inod) 254 END DO 255 256 isU=INU(i-1)+1 257 ieU=INU(i) 258 DO j=isU,ieU 259 inod=IAU(j)+ZERO_ORIGIN 260 v = v - AU(j)*X(inod) 261 END DO 262 263 X(i) = v / D(i) 264 END DO 265 DO i=N,1,-1 266 v = B(i) 267 268 isL=INL(i-1)+1 269 ieL=INL(i) 270 DO j=isL,ieL 271 inod=IAL(j)+ZERO_ORIGIN 272 v = v - AL(j)*X(inod) 273 END DO 274 275 isU=INU(i-1)+1 276 ieU=INU(i) 277 DO j=isU,ieU 278 inod=IAU(j)+ZERO_ORIGIN 279 v = v - AU(j)*X(inod) 280 END DO 281 282 X(i) = v / D(i) 283 END DO 284 END DO 285 END SUBROUTINE sgs 286 287 SUBROUTINE matrix_counting(LEVEL_NUM,SOLVER_COMM,my_rank) 288 USE data_structure_for_AMG 289 IMPLICIT NONE 290 291 integer(kind=kint) ,intent(in) :: LEVEL_NUM,SOLVER_COMM,my_rank 292 integer(kind=kint) :: i,j,nonzero_l,nonzero_g,level,n,np 293 REAL,DIMENSION(:),allocatable::comm_rate 294 REAL:: rate,crate_l,crate_g 295 296 allocate(comm_rate(LEVEL_NUM)) 297 DO level=1,LEVEL_NUM 298 n= HIERARCHICAL_DATA(level) % N 299 np=HIERARCHICAL_DATA(level) % NP 300 !C in the case of n == 0 301 if(np == 0) then 302 crate_l=1 303 else 304 crate_l=REAL(n)/np 305 end if 306 crate_g=0 307 308 crate_g=crate_l 309 comm_rate(level) = crate_g 310 END DO 311 312 nonzero_l=0 313 314 DO level=1,LEVEL_NUM 315 n= HIERARCHICAL_DATA(level) % N 316 nonzero_l=nonzero_l+1 + HIERARCHICAL_DATA(level)%INU(n) 317 nonzero_l=nonzero_l + HIERARCHICAL_DATA(level)%INL(n) 318 END DO 319 320 nonzero_g=0 321 322 nonzero_g=nonzero_l 323 i=nonzero_g 324 325 326 n=HIERARCHICAL_DATA(1)%N 327 328 !C in the case of n == 0 329 nonzero_l = 0 330 if(n>0) then 331 nonzero_l=1+HIERARCHICAL_DATA(1)%INU(n)+HIERARCHICAL_DATA(1)%INL(n) 332 end if 333 334 335 nonzero_g = 0 336 337 nonzero_g=nonzero_l 338 j=nonzero_g 339 340 if(my_rank==0)then 341 rate=REAL(i)/j 342 write(*,*) "nonzeros in all level:", i, "rate of level:",rate 343 write(*,*) "worst rate of own nodes of vector" 344 write(*,*) comm_rate 345 end if 346 END SUBROUTINE matrix_counting 347 348 !C Symmetrize the matrix especially the rows: N+1..NP 349 SUBROUTINE matrix_arrange (N,NP,D,AL,INL,IAL,AU,INU,IAU) 350 351 IMPLICIT NONE 352#include "precision.inc" 353 354 INTEGER(kind=kint ), intent(in) :: N,NP 355 REAL (kind=kreal), intent(inout),DIMENSION(:) :: D 356 REAl (kind=kreal), intent(inout),DIMENSION(:) :: AU,AL 357 INTEGER(kind=kint ), intent(in),DIMENSION(:) :: IAU,IAL 358 INTEGER(kind=kint ), intent(in),DIMENSION(0:) :: INU,INL 359 360 INTEGER(kind=kint ) :: i,j,isU,ieU,isL,ieL,k,column 361 logical :: flag 362 363 364 DO j=N+1,NP 365 isL = INL(j-1)+1 366 ieL = INL(j) 367 DO i=isL,ieL 368 column = IAL(i)+ZERO_ORIGIN 369 if(column>N) then 370 AL(i)=0 371 else if(column>0 .and. column<=N) then 372 isU = INU(column-1)+1 373 ieU = INU(column) 374 flag = .false. 375 DO k=isU,ieU 376 if(IAU(k)+ZERO_ORIGIN==j)then 377 AL(i)=AU(k) 378 flag= .true. 379 exit 380 end if 381 end DO 382 if(.not.flag) then 383 AL(i) =0.0 384 end if 385 end if 386 END DO 387 388 isU = INU(j-1)+1 389 ieU = INU(j) 390 DO i=isU,ieU 391 AU(i)= 0.0 392 END DO 393 END DO 394 END SUBROUTINE matrix_arrange 395 396 397 !C 398 !C*** CG 399 !C 400 SUBROUTINE AMGCG (N, NP, NPL, NPU, & 401 & D, AL, INL, IAL, AU, INU, IAU, & 402 & B, X, resid,iter) 403 404 USE data_structure_for_AMG 405 USE data_creation_AMGCG 406 IMPLICIT NONE 407 INTEGER(kind=kint ), intent(in) :: N,NP,NPL,NPU 408 REAL (kind=kreal) :: RESID 409 REAL (kind=kreal) :: SIGMA_DIAG 410 REAL (kind=kreal) :: SIGMA 411 INTEGER(kind=kint ) :: ITER 412 INTEGER(kind=kint ) :: ERROR 413 INTEGER(kind=kint ) :: my_rank 414 415 INTEGER(kind=kint ) :: NPROCS,SOLVER_COMM 416 417 REAL (kind=kreal), DIMENSION(: ), pointer :: D 418 REAL (kind=kreal), DIMENSION(: ), pointer :: B 419 REAL (kind=kreal), DIMENSION(NP ) :: X 420 421 REAL (kind=kreal), DIMENSION(:), pointer :: AU 422 REAL (kind=kreal), DIMENSION(:), pointer :: AL 423 424 INTEGER(kind=kint ), DIMENSION(:), pointer :: INU 425 INTEGER(kind=kint ), DIMENSION(:), pointer :: IAU 426 INTEGER(kind=kint ), DIMENSION(:), pointer :: INL 427 INTEGER(kind=kint ), DIMENSION(:), pointer :: IAL 428 429 430 REAL (kind=kreal), DIMENSION(:,:), ALLOCATABLE :: WW 431 REAL (kind=kreal), DIMENSION(:), ALLOCATABLE :: Temp 432 433 REAL (kind=kreal), DIMENSION(:), ALLOCATABLE, SAVE :: DD 434 REAL (kind=kreal), DIMENSION(:), ALLOCATABLE, SAVE :: SCALE 435 436 INTEGER(kind=kint ) :: P, Q, R, Z, MAXIT, IFLAG=0 437 REAL (kind=kreal) :: TOL, W, SS 438 439 INTEGER(kind=kint) :: LEVEL_NUM,cnt,NEIBPETOT,I,J,ISU,IEU,ISL,IEL,WSIZE 440 INTEGER(kind=kint) :: INOD 441 REAL (kind=kreal) :: WVAL,BNRM20,BNRM2,RHO0,RHO,BETA,RHO1,C10,C1,ALPHA,DNRM20 442 REAL (kind=kreal) :: DNRM2,X1,X2 443 REAL (kind=kreal) :: STARTTIME,ENDTIME,RTIME,STARTTIME2,ENDTIME2,RTIME2 444 445 REAL (kind=kreal) :: theta 446 !C-- theta means the critierion for strong connection. 447 !C-- if theta=0 then all nonzero elements are recognized as 448 !C-- strong connection. 449 450 !C-- INIT. 451 ERROR= 0 452 NPROCS=0 453 my_rank=0 454 455 ALLOCATE (WW(NP,3)) 456 457 458 ALLOCATE ( Temp(NP) ) 459 460 R = 1 461 Z = 2 462 Q = 2 463 P = 3 464 465 MAXIT = ITER 466 TOL = RESID 467 468 CALL data_creation(NP, NPL, NPU, D, AL, INL, IAL, AU, INU, IAU, LEVEL_NUM, theta) 469 470 471 !C +-----------------------+ 472 !C | {r0}= {b} - [A]{xini} | 473 !C +-----------------------+ 474 !C=== 475 !C-- BEGIN calculation 476 do j= 1, N 477 WVAL= B(j) - D(j) * X(j) 478 isU= INU(j-1) + 1 479 ieU= INU(j ) 480 481 do i= isU, ieU 482 inod= IAU(i)+ZERO_ORIGIN 483 WVAL= WVAL - AU(i) * X(inod) 484 enddo 485 486 isL= INL(j-1) + 1 487 ieL= INL(j ) 488 do i= isL, ieL 489 inod= IAL(i)+ZERO_ORIGIN 490 WVAL= WVAL - AL(i) * X(inod) 491 enddo 492 WW(j,R)= WVAL 493 enddo 494 495 BNRM20 = 0.d0 496 do i= 1, N 497 BNRM20 = BNRM20 + B(i)**2 498 enddo 499 500 BNRM2=BNRM20 501 502 if (BNRM2.eq.0.d0) BNRM2= 1.d0 503 ITER = 0 504 505 !C=== 506 do iter= 1, MAXIT 507 !C 508 !C************************************************* Conjugate Gradient Iteration 509 510 !C 511 !C +----------------+ 512 !C | {z}= [Minv]{r} | 513 !C +----------------+ 514 !C=== 515 do i=1,NP 516 WW(i,Z)=0.0 517 end do 518 CALL v_cycle(N, WW(:,R), WW(:,Z), LEVEL_NUM, Temp) 519!!$ WW(:,Z) = WW(:,R) 520 !C=== 521 522 !C 523 !C +---------------+ 524 !C | {RHO}= {r}{z} | 525 !C +---------------+ 526 !C=== 527 RHO0= 0.0 528 529 do i= 1, N 530 RHO0= RHO0 + WW(i,R)*WW(i,Z) 531 enddo 532 533 534 RHO=RHO0 535 536 !C=== 537 538 !C 539 !C +-----------------------------+ 540 !C | {p} = {z} if ITER=1 | 541 !C | BETA= RHO / RHO1 otherwise | 542 !C +-----------------------------+ 543 !C=== 544 if ( ITER.eq.1 ) then 545 do i= 1, N 546 WW(i,P)= WW(i,Z) 547 enddo 548 else 549 BETA= RHO / RHO1 550 do i= 1, N 551 WW(i,P)= WW(i,Z) + BETA*WW(i,P) 552 enddo 553 endif 554 !C=== 555 556 !C 557 !C +-------------+ 558 !C | {q}= [A]{p} | 559 !C +-------------+ 560 !C=== 561 562 563 !C 564 !C-- BEGIN calculation 565 566 do j= 1, N 567 WVAL= D(j) * WW(j,P) 568 569 isU= INU(j-1) + 1 570 ieU= INU(j ) 571 572 do i= isU, ieU 573 inod= IAU(i)+ZERO_ORIGIN 574 WVAL= WVAL + AU(i) * WW(inod,P) 575 enddo 576 577 isL= INL(j-1) + 1 578 ieL= INL(j ) 579 580 do i= isL, ieL 581 inod= IAL(i)+ZERO_ORIGIN 582 WVAL= WVAL + AL(i) * WW(inod,P) 583 enddo 584 WW(j,Q)= WVAL 585 enddo 586 !C=== 587 588 !C 589 !C +---------------------+ 590 !C | ALPHA= RHO / {p}{q} | 591 !C +---------------------+ 592 !C=== 593 C10= 0.d0 594 595 do i= 1, N 596 C10= C10 + WW(i,P)*WW(i,Q) 597 enddo 598 599 C1=C10 600 601 ALPHA= RHO / C1 602 !C=== 603 604 !C 605 !C +----------------------+ 606 !C | {x}= {x} + ALPHA*{p} | 607 !C | {r}= {r} - ALPHA*{q} | 608 !C +----------------------+ 609 !C=== 610 DNRM20= 0.d0 611 612 X1= 0.0d0 613 X2= 0.0d0 614 615 do i= 1, N 616 X(i) = X (i) + ALPHA * WW(i,P) 617 WW(i,R)= WW(i,R) - ALPHA * WW(i,Q) 618 enddo 619 620 DNRM20 = 0.0 621 do i= 1, N 622 DNRM20= DNRM20 + WW(i,R)**2 623 enddo 624 625 626 DNRM2=DNRM20 627 RESID= sqrt(DNRM2/BNRM2) 628 629 630!!$#ifdef PRINT_REZ 631 if (my_rank.eq.0) write (*, 1000) ITER, RESID 6321000 format (i5, 1pe16.6) 633!!$#endif 634 if ( RESID.le.TOL ) then 635 exit 636 end if 637 if ( ITER .eq.MAXIT ) ERROR= -100 638 639 RHO1 = RHO 640 end do 641 !C=== 642 643 644 call clear_matrix(LEVEL_NUM) 645 646 647 !C count non zeros of hierarchy of matrixes 648 DEALLOCATE (WW) 649 DEALLOCATE (Temp) 650 651 END SUBROUTINE AMGCG 652 653 subroutine clear_matrix(LEVEL_NUM) 654 use data_structure_for_AMG 655 656 IMPLICIT NONE 657 INTEGER(kind=kint) :: i, LEVEL_NUM 658 659 NULLIFY( HIERARCHICAL_DATA(1) % INU, & 660 & HIERARCHICAL_DATA(1) % INL, & 661 & HIERARCHICAL_DATA(1) % IAU, & 662 & HIERARCHICAL_DATA(1) % IAL, & 663 & HIERARCHICAL_DATA(1) % AU, & 664 & HIERARCHICAL_DATA(1) % AL, & 665 & HIERARCHICAL_DATA(1) % D, & 666 & HIERARCHICAL_DATA(1) % X, & 667 & HIERARCHICAL_DATA(1) % B ) 668 669 DO i=2,LEVEL_NUM 670 deallocate(HIERARCHICAL_DATA(i) % INU, & 671 & HIERARCHICAL_DATA(i) % INL, & 672 & HIERARCHICAL_DATA(i) % IAU, & 673 & HIERARCHICAL_DATA(i) % IAL, & 674 & HIERARCHICAL_DATA(i) % AU, & 675 & HIERARCHICAL_DATA(i) % AL, & 676 & HIERARCHICAL_DATA(i) % D, & 677 & HIERARCHICAL_DATA(i) % X, & 678 & HIERARCHICAL_DATA(i) % B, & 679 & HIERARCHICAL_DATA(i) % R % IN, & 680 & HIERARCHICAL_DATA(i) % R % CN, & 681 & HIERARCHICAL_DATA(i) % R % V, & 682 & HIERARCHICAL_DATA(i) % P % IN, & 683 & HIERARCHICAL_DATA(i) % P % CN, & 684 & HIERARCHICAL_DATA(i) % P % V ) 685 deallocate(HIERARCHICAL_DATA(i) % P, & 686 & HIERARCHICAL_DATA(i) % R ) 687 END DO 688 689 690 end subroutine clear_matrix 691 692END MODULE solver_AMGCG 693 694 695