C Parallel Block ICCG Method (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 BICCG(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 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 INPS(0:NPM-1),INPE(0:NPM-1) DIMENSION IRNEW(0:NPM-1),IDNEW(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 各PE担当行数の設定 IT=(N-1)/NUMPROCS+1 DO 22 I=0,NUMPROCS-1 INPS(I)=I*IT+1 INPE(I)=(I+1)*IT 22 CONTINUE INPE(NUMPROCS-1)=N DO 1 I=1,N DIAG(I)=0D0 IUHEAD(I)=0 1 CONTINUE DO 2 I=1,NNONZERO IUCOL(I)=0 U(I)=0D0 2 CONTINUE KU=0 IUHEAD(INPS(MYID))=1 DO 700 I=INPS(MYID),INPE(MYID) C DO 700 I=1,N KU=0 KL=0 DO 710 J=ROWPTR(I),ROWPTR(I+1)-1 JJ=COLIND(J) IF (JJ.GE.INPS(MYID).AND.JJ.LE.INPE(MYID)) THEN IF (JJ.EQ.I) THEN DIAG(I)=VAL(J) ELSEIF (JJ.LT.I) THEN KL=KL+1 ELSEIF (JJ.GT.I) THEN IUCOL(KU+IUHEAD(I))=JJ U(KU+IUHEAD(I))=VAL(J) KU=KU+1 ELSE WRITE(6,*) 'ERROR' STOP ENDIF ENDIF IUHEAD(I+1)=IUHEAD(I)+KU 710 CONTINUE 700 CONTINUE C GANMA: 加速係数 GANMA=1.3D0 DO 145 I=INPS(MYID),INPE(MYID) DIAG(I)=DIAG(I)*GANMA 145 CONTINUE DO 800 I=INPS(MYID),INPE(MYID) 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) 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 C IF (ABS(DIAG(I)).LT.1D-3) THEN C WRITE(6,*) 'DIAG ERROR',I,DIAG(I) C STOP C ENDIF 800 CONTINUE DO 171 I=0,NUMPROCS-1 IRNEW(I)=INPE(I)-INPS(I)+1 IDNEW(I)=INPS(I)-1 171 CONTINUE BVNORM=0D0 DO 260 I=INPS(MYID),INPE(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 DO 20 I=INPS(MYID),INPE(MYID) 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 CGROP=0D0 DO 1001 ITE=1,MAXITE C DO 1001 ITE=1,1 DO 312 I=INPS(MYID),INPE(MYID) Z(I)=R(I) 312 CONTINUE DO 301 I=INPS(MYID),INPE(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 J1=INPE(MYID) Z(J1)=Z(J1)/DIAG(J1) DO 305 I=INPE(MYID)-1,INPS(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 200 I=INPS(MYID),INPE(MYID) CGROP=CGROP+R(I)*Z(I) 200 CONTINUE CALL MPI_ALLREDUCE(CGROP,CGROP2,1, &MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,IEER) CGROP=CGROP2 IF (ITE.EQ.1) THEN DO 201 I=INPS(MYID),INPE(MYID) PN(I)=Z(I) 201 CONTINUE ELSE BETA=CGROP/CGROPP DO 202 I=INPS(MYID),INPE(MYID) PN(I)=Z(I)+BETA*PN(I) 202 CONTINUE ENDIF INSTAT=INPS(MYID) ITD=INPE(MYID)-INPS(MYID)+1 CALL MPI_ALLGATHERV(PN(INSTAT),ITD, &MPI_DOUBLE_PRECISION,PN(1),IRNEW(0),IDNEW(0), &MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IEER) DO 50 I=INPS(MYID),INPE(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=INPS(MYID),INPE(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=INPS(MYID),INPE(MYID) SOLX(I)=SOLX(I)+ALPHA*PN(I) R(I)=R(I)-ALPHA*Q(I) 1500 CONTINUE RNORM=0D0 DO 1600 I=INPS(MYID),INPE(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 INSTAT=INPS(MYID) ITD=INPE(MYID)-INPS(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