C Parallel ICCG Method with Renumbering Process (CRS) Ver. 1.0 C C　 並列化ICCG法 C　 不完全コレスキー分解前処理 C 解ベクトル初期値は０とする C C ICCG METHOD CRS形式による格納 C 長谷川里見、長谷川秀彦、藤野清次訳、「反復法TEMPLATE」p.78　参照 C ISBN4-254-11401-X C C COLIND 列のインデクス　VAL 非ゼロ要素の値 C　　 ROWPTR i行目に対応したCOLINDの先頭の要素番号を表す C NA 配列ROWPTRの次元数　 C　　 NNONZERO 配列COLIND, VALの次元数 C　　 B　右辺ベクトル C 　　N　次元数　 C 　　SOLX　解ベクトル C ERR 打ち切り誤差 C C NPM 最大プロセッサ数 C MAXITE 最大反復数 C C 京都大学大型計算機センター C　　　助手　　　　　岩下　武史　　製作　　 C　　　最終修正日　平成１２年１０月１０日 C　　 C　　　質問などの宛先 take@kudpc.kyoto-u.ac.jp C C SUBROUTINE PICCGRP(ROWPTR,COLIND,VAL,B,N,SOLX,ERR) IMPLICIT REAL*8(A-H,O-Z) INCLUDE 'param.h' INCLUDE 'mpif.h' PARAMETER(NPM=500,MAXITE=5000) INTEGER STAT(MPI_STATUS_SIZE) INTEGER ROWPTR(NA+1),COLIND(NNONZERO) DIMENSION IUHEAD(NA+1),IUCOL(NNONZERO) C REAL*8 VAL(NNONZERO) REAL*8 B(N),SOLX(N) REAL*8 P(NA),Q(NA),PN(NA),R(NA) REAL*8 CGROPP,CGROP REAL*8 ALPHA,ALPHAT,BETA,AR0 REAL*8 DIAG(NA),U(NNONZERO),Z(NA) C DIMENSION INPSR(0:NPM-1),INPER(0:NPM-1) DIMENSION INPSBR(0:NPM-1),INPEBR(0:NPM-1) DIMENSION IRNEW(0:NPM-1),IDNEW(0:NPM-1) DIMENSION IRNEW2(0:NPM-1),IDNEW2(0:NPM-1) C CALL MPI_INIT(IEER) CALL MPI_COMM_RANK(MPI_COMM_WORLD,MYID,IEER) CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NUMPROCS,IEER) C For Symmetric Matrix NB=0 DO 333 I=1,N JMAX=0 DO 334 J=ROWPTR(I),ROWPTR(I+1)-1 JMAX=MAX(JMAX,COLIND(J)) 334 CONTINUE JBAND=JMAX-I NB=MAX(NB,JBAND) 333 CONTINUE C WRITE(6,*) 'BAND WIDTH=',NB NIN=(N-(NUMPROCS-1)*NB)/NUMPROCS IF (NIN.LE.NB) THEN WRITE(6,*) 'INTERIOR PART < INTERFACE PART' WRITE(6,*) 'BAND WIDTH IS TOO LARGE' WRITE(6,*) 'STOP' WRITE(6,*) 'PLEASE SELECT BLOCK TYPE APPROACH' STOP ENDIF C 内部領域、インターフェース領域の設定 DO 100 I=0,NUMPROCS-2 INPSR(I)=(NIN+NB)*I+1 INPER(I)=(NIN+NB)*I+NIN INPSBR(I)=(NIN+NB)*I+NIN+1 INPEBR(I)=(NIN+NB)*(I+1) 100 CONTINUE INPSR(NUMPROCS-1)=(NIN+NB)*(NUMPROCS-1)+1 INPER(NUMPROCS-1)=N INPSBR(NUMPROCS-1)=N+1 INPEBR(NUMPROCS-1)=N DO 1 I=1,N DIAG(I)=0D0 IUHEAD(I)=0 1 CONTINUE IUHEAD(N+1)=0 DO 2 I=1,NNONZERO IUCOL(I)=0 U(I)=0D0 2 CONTINUE KU=0 IUHEAD(INPSR(MYID))=1 DO 155 I=INPSR(MYID),INPSR(MYID)+NB KU=0 DO 156 J=ROWPTR(I),ROWPTR(I+1)-1 JJ=COLIND(J) IF (JJ.GE.INPSR(MYID)-NB.AND.JJ.LE.INPSR(MYID)-1) THEN IUCOL(KU+IUHEAD(I))=JJ U(KU+IUHEAD(I))=VAL(J) KU=KU+1 ELSEIF(JJ.GE.I) THEN IF (JJ.NE.I) THEN IUCOL(KU+IUHEAD(I))=JJ U(KU+IUHEAD(I))=VAL(J) KU=KU+1 ELSE DIAG(I)=VAL(J) ENDIF ENDIF 156 CONTINUE IUHEAD(I+1)=IUHEAD(I)+KU 155 CONTINUE DO 700 I=INPSR(MYID)+NB+1,INPER(MYID) KU=0 DO 710 J=ROWPTR(I),ROWPTR(I+1)-1 JJ=COLIND(J) IF (JJ.EQ.I) THEN DIAG(I)=VAL(J) ELSEIF (JJ.GT.I) THEN IUCOL(KU+IUHEAD(I))=JJ U(KU+IUHEAD(I))=VAL(J) KU=KU+1 ENDIF 710 CONTINUE IUHEAD(I+1)=IUHEAD(I)+KU 700 CONTINUE IF (MYID.NE.NUMPROCS-1) THEN DO 170 I=INPSBR(MYID),INPEBR(MYID) KU=0 DO 171 J=ROWPTR(I),ROWPTR(I+1)-1 JJ=COLIND(J) IF (JJ.GE.I.AND.JJ.LE.INPEBR(MYID)) THEN IF (JJ.NE.I) THEN IUCOL(KU+IUHEAD(I))=JJ U(KU+IUHEAD(I))=VAL(J) KU=KU+1 ELSE DIAG(I)=VAL(J) ENDIF ENDIF 171 CONTINUE IUHEAD(I+1)=IUHEAD(I)+KU 170 CONTINUE DO 150 I=INPEBR(MYID)+1,MIN(INPEBR(MYID)+NB,N) KU=0 DO 151 J=ROWPTR(I),ROWPTR(I+1)-1 JJ=COLIND(J) IF(JJ.GE.INPSBR(MYID).AND.JJ.LE.INPEBR(MYID))THEN IUCOL(KU+IUHEAD(I))=JJ U(KU+IUHEAD(I))=VAL(J) KU=KU+1 ELSEIF(JJ.GE.I) THEN IF (JJ.NE.I) THEN IUCOL(KU+IUHEAD(I))=JJ U(KU+IUHEAD(I))=VAL(J) KU=KU+1 ELSE DIAG(I)=VAL(J) ENDIF ENDIF 151 CONTINUE IUHEAD(I+1)=IUHEAD(I)+KU 150 CONTINUE ENDIF CALL MPI_BARRIER(MPI_COMM_WORLD,IEER) C STOP C GANMA: 加速係数 C GANMA=1.03D0 C DO 145 I=1,N C DIAG(I)=DIAG(I)*GANMA C 145 CONTINUE DO 810 I=INPEBR(MYID)+1,MIN(INPEBR(MYID)+NB,N) DO 820 J=IUHEAD(I),IUHEAD(I+1)-1 JJ=IUCOL(J) DIAG(IUCOL(J))=DIAG(IUCOL(J))-U(J)*U(J)/DIAG(I) DO 830 JP=IUHEAD(I),IUHEAD(I+1)-1 JJP=IUCOL(JP) C 以下のIF文はなくてもよい。 C IF (JJP.NE.JJ) THEN DO 840 JI=IUHEAD(JJ),IUHEAD(JJ+1)-1 IF (IUCOL(JI).EQ.JJP) THEN U(JI)=U(JI)-U(J)*U(JP)/DIAG(I) GO TO 848 ENDIF 840 CONTINUE ENDIF 848 CONTINUE 830 CONTINUE 820 CONTINUE 810 CONTINUE DO 1190 I=INPSR(MYID),INPSR(MYID)+NB DO 1191 J=IUHEAD(I),IUHEAD(I+1)-1 JJ=IUCOL(J) IF (JJ.GE.INPSR(MYID)) THEN DIAG(IUCOL(J))=DIAG(IUCOL(J))-U(J)*U(J)/DIAG(I) DO 1192 JP=IUHEAD(I),IUHEAD(I+1)-1 JJP=IUCOL(JP) C 以下のIF文はなくてもよい。 C IF (JJP.NE.JJ) THEN DO 1193 JI=IUHEAD(JJ),IUHEAD(JJ+1)-1 IF (IUCOL(JI).EQ.JJP) THEN U(JI)=U(JI)-U(J)*U(JP)/DIAG(I) GO TO 1198 ENDIF 1193 CONTINUE ENDIF 1198 CONTINUE 1192 CONTINUE ENDIF 1191 CONTINUE 1190 CONTINUE DO 200 I=INPSR(MYID)+NB+1,INPEBR(MYID) DO 201 J=IUHEAD(I),IUHEAD(I+1)-1 JJ=IUCOL(J) DIAG(IUCOL(J))=DIAG(IUCOL(J))-U(J)*U(J)/DIAG(I) DO 202 JP=IUHEAD(I),IUHEAD(I+1)-1 JJP=IUCOL(JP) IF (JJP.GT.JJ) THEN DO 203 JI=IUHEAD(JJ),IUHEAD(JJ+1)-1 IF (IUCOL(JI).EQ.JJP) THEN U(JI)=U(JI)-U(J)*U(JP)/DIAG(I) GO TO 208 ENDIF 203 CONTINUE ENDIF 208 CONTINUE 202 CONTINUE 201 CONTINUE 200 CONTINUE DO 230 I=INPSR(MYID),INPEBR(MYID) IF (ABS(DIAG(I)).LT.1D-3) THEN WRITE(6,*) 'DIAG ERROR' ,I,DIAG(I) STOP ENDIF 230 CONTINUE C 不完全これスキー分解終了 BVNORM=0D0 DO 260 I=INPSR(MYID),INPEBR(MYID) BVNORM=BVNORM+B(I)*B(I) 260 CONTINUE CALL MPI_ALLREDUCE(BVNORM,BVNO2,1, &MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IEER) BNORM=DSQRT(BVNO2) DO 11 I=1,N SOLX(I)=0D0 11 CONTINUE C DO 20 I=INPSR(MYID),INPEBR(MYID) DO 20 I=1,N AR0=0D0 DO 30 J=ROWPTR(I),ROWPTR(I+1)-1 JJ=COLIND(J) AR0=AR0+VAL(J)*SOLX(JJ) 30 CONTINUE R(I)=B(I)-AR0 20 CONTINUE IUP=MYID-1 IDOWN=MYID+1 IF(MYID.EQ.0) IUP=MPI_PROC_NULL IF(MYID.EQ.NUMPROCS-1) IDOWN=MPI_PROC_NULL CGROP=0D0 DO 1001 ITE=1,MAXITE C DO 1001 ITE=1,1 DO 312 I=INPSR(MYID),INPEBR(MYID) Z(I)=R(I) 312 CONTINUE DO 301 I=INPSR(MYID),INPER(MYID) DO 302 J=IUHEAD(I),IUHEAD(I+1)-1 JJ=IUCOL(J) Z(JJ)=Z(JJ)-Z(I)*U(J)/DIAG(I) 302 CONTINUE 301 CONTINUE INSFST=INPSR(MYID) INSD=NB INSEND=INPEBR(MYID)+1 IF (MYID.EQ.NUMPROCS-1) INSEND=N CALL MPI_SENDRECV(Z(INSFST),INSD,MPI_DOUBLE_PRECISION,IUP &,0,Z(INSEND),INSD,MPI_DOUBLE_PRECISION,IDOWN &,0,MPI_COMM_WORLD,STAT,IEER) IF (MYID.NE.NUMPROCS-1) THEN DO 451 I=INPEBR(MYID)+1,MIN(INPEBR(MYID)+NB,N) DO 452 J=IUHEAD(I),IUHEAD(I+1)-1 JJ=IUCOL(J) C 以下のIF文はいらない。ない場合は関係ないところが引かれる。 IF (JJ.LT.I) THEN Z(JJ)=Z(JJ)-Z(I)*U(J)/DIAG(I) ENDIF 452 CONTINUE 451 CONTINUE DO 351 I=INPSBR(MYID),INPEBR(MYID) DO 352 J=IUHEAD(I),IUHEAD(I+1)-1 JJ=IUCOL(J) Z(JJ)=Z(JJ)-Z(I)*U(J)/DIAG(I) 352 CONTINUE 351 CONTINUE ENDIF C 前進代入終了 IF(MYID.NE.NUMPROCS-1)THEN Z(INPEBR(MYID))=Z(INPEBR(MYID))/DIAG(INPEBR(MYID)) DO 630 I=INPEBR(MYID)-1,INPSBR(MYID),-1 DO 631 J=IUHEAD(I),IUHEAD(I+1)-1 JJ=IUCOL(J) Z(I)=Z(I)-U(J)*Z(JJ) 631 CONTINUE Z(I)=Z(I)/DIAG(I) 630 CONTINUE ENDIF INSFST=INPSBR(MYID) IF (MYID.EQ.NUMPROCS-1) INSFST=N INSD=NB INSEND=INPSR(MYID)-NB IF (MYID.EQ.0) INSEND=1 CALL MPI_SENDRECV(Z(INSFST),INSD,MPI_DOUBLE_PRECISION,IDOWN, &1,Z(INSEND),INSD,MPI_DOUBLE_PRECISION,IUP, &1,MPI_COMM_WORLD,STAT,IEER) DO 305 I=INPER(MYID),INPSR(MYID),-1 DO 306 J=IUHEAD(I),IUHEAD(I+1)-1 JJ=IUCOL(J) Z(I)=Z(I)-U(J)*Z(JJ) 306 CONTINUE 316 Z(I)=Z(I)/DIAG(I) 305 CONTINUE C DO 395 I=1,N C Z(I)=R(I) C 395 CONTINUE CGROPP=CGROP CGROP=0D0 DO 265 I=INPSR(MYID),INPEBR(MYID) CGROP=CGROP+R(I)*Z(I) 265 CONTINUE CALL MPI_ALLREDUCE(CGROP,CGROP2,1, &MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IEER) CGROP=CGROP2 IF (ITE.EQ.1) THEN DO 261 I=INPSR(MYID),INPEBR(MYID) PN(I)=Z(I) 261 CONTINUE ELSE BETA=CGROP/CGROPP DO 262 I=INPSR(MYID),INPEBR(MYID) PN(I)=Z(I)+BETA*PN(I) 262 CONTINUE ENDIF INSFST=INPSR(MYID) INSD=NB INSEND=INPEBR(MYID)+1 IF (MYID.EQ.NUMPROCS-1) INSEND=N C WRITE(6,*) MYID,INSEND C CALL MPI_BARRIER(MPI_COMM_WORLD,IEER) C STOP C STOP CALL MPI_SENDRECV(PN(INSFST),INSD,MPI_DOUBLE_PRECISION, &IUP,2,PN(INSEND),INSD,MPI_DOUBLE_PRECISION, &IDOWN,2,MPI_COMM_WORLD,STAT,IEER) INSFST=INPSBR(MYID) IF (MYID.EQ.NUMPROCS-1) INSFST=N INSD=NB INSEND=INPSR(MYID)-NB IF (MYID.EQ.0) INSEND=1 CALL MPI_SENDRECV(PN(INSFST),INSD,MPI_DOUBLE_PRECISION, &IDOWN,3,PN(INSEND),INSD,MPI_DOUBLE_PRECISION, &IUP,3,MPI_COMM_WORLD,STAT,IEER) DO 50 I=INPSR(MYID),INPEBR(MYID) Q(I)=0D0 DO 51 J=ROWPTR(I),ROWPTR(I+1)-1 JJ=COLIND(J) Q(I)=Q(I)+VAL(J)*PN(JJ) 51 CONTINUE 50 CONTINUE ALPHAT=0D0 DO 1400 I=INPSR(MYID),INPEBR(MYID) ALPHAT=ALPHAT+PN(I)*Q(I) 1400 CONTINUE CALL MPI_ALLREDUCE(ALPHAT,ALPHA2,1, &MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IEER) ALPHA=CGROP/ALPHA2 DO 1500 I=INPSR(MYID),INPEBR(MYID) SOLX(I)=SOLX(I)+ALPHA*PN(I) R(I)=R(I)-ALPHA*Q(I) 1500 CONTINUE RNORM=0D0 DO 1600 I=INPSR(MYID),INPEBR(MYID) RNORM=RNORM+R(I)*R(I) 1600 CONTINUE CALL MPI_ALLREDUCE(RNORM,RNORM2,1, &MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IEER) RNORM=DSQRT(RNORM2) CCCCCCCCCCCCCCCCCCCCCCCCCC C 収束の振る舞いを出力しない場合には下行をコメントアウト IF (MYID.EQ.NUMPROCS-1) THEN WRITE(6,*) ITE,RNORM/BNORM ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCC IF (RNORM/BNORM.LT.ERR) GO TO 400 1001 CONTINUE 400 CONTINUE C 通信用配列の設定 DO 241 I=0,NUMPROCS-2 IRNEW(I)=NIN+NB IDNEW(I)=(NIN+NB)*I 241 CONTINUE IRNEW(NUMPROCS-1)=N-NB*(NUMPROCS-1)-NIN*(NUMPROCS-1) IDNEW(NUMPROCS-1)=(NIN+NB)*(NUMPROCS-1) INSTAT=INPSR(MYID) ITD=INPEBR(MYID)-INPSR(MYID)+1 CALL MPI_ALLGATHERV(SOLX(INSTAT),ITD, &MPI_DOUBLE_PRECISION,SOLX(1),IRNEW(0),IDNEW(0), &MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IEER) CALL MPI_FINALIZE(IEER) RETURN END