PROGRAM EXERCISE36 IMPLICIT NONE ! ヤコビ法で行列を対角化するプログラム INTEGER, PARAMETER :: N = 10 !行列のサイズはここで決めているがFORMAT文の大きさはこのNとリンクしていない INTEGER :: i, j, p, q, iter DOUBLE PRECISION, DIMENSION(1:N,1:N) :: Amat, Amat_new, Umat, UDUmat, Umat_new, UTUmat DOUBLE PRECISION :: alpha, beta, sintheta, costheta, t, maxoffdiag ! 行列Aの初期化 DO i = 1, N Amat(i,i:N) = DBLE(i) Amat(i:N,i) = DBLE(i) END DO ! 行列Aが正しく初期化されたか出力してチェック PRINT *, "Amatrix initialized" DO i = 1, N WRITE(*,'(10ES15.5)') (Amat(i,j),j = 1, N) END DO ! 行列Uの初期化(単位行列) Umat(1:N,1:N) = 0.0d0 DO i = 1, N Umat(i,i) = 1.0d0 END DO ! 反復の回数のカウンタをゼロにセット iter = 0 ! 反復法スタート DO !最大の非対角要素を与えるp,qを決定 !まず非対角要素を格納する変数の値をゼロとして maxoffdiag = 0.0d0 ! (i,j)がすべての非対角要素の値を参照するようにDOループを組む DO i = 1, N-1 DO j = i+1, N ! すでに保存されているmaxoffdiagの値よりも非対角要素Amat(i,j)の絶対値が大きい場合は IF( ABS(Amat(i,j)) .GT. maxoffdiag) THEN p = i; q = j ! p,qにこの行列要素の値をセット maxoffdiag = ABS(Amat(i,j)) ! maxoffdiagの値を更新 END IF END DO END DO ! このDOループが終了した段階で最大の非対角行列要素がmaxoffdiagに格納され、その添字がp,qに格納される。 PRINT *, "iter = ", iter, "maxoffdiag = ", maxoffdiag, "p = ", p, "q = ", q ! 最大の非対角行列要素が十分に小さい場合は対角化が完了。反復を終了 IF( maxoffdiag .LT. 1.0D-10) EXIT ! alphaを計算 alpha = (Amat(q,q) - Amat(p,p)) / (2.0d0*Amat(p,q)) ! tをalphaから計算 t = SIGN(1.0d0,alpha) / ( ABS(alpha) + SQRT(alpha**2 + 1.0d0)) PRINT *, "alpha = ", alpha, "t = ", t ! sin theta, cos theta, betaの値をtから計算 costheta = 1.0d0/ SQRT( 1.0d0 + t**2) sintheta = t*costheta beta = sintheta / ( 1.0d0 + costheta) ! AmatからAmat_newを計算 ! ほとんどの項は変わらないためまずはAmatを代入 Amat_new(1:N,1:N) = Amat(1:N,1:N) ! pp と qq成分 Amat_new(p,p) = Amat(p,p) - t * Amat(p,q) Amat_new(q,q) = Amat(q,q) + t * Amat(p,q) ! ip, iq, pj, qj成分 DO i = 1, N IF (i .EQ. p .OR. i .EQ. q) CYCLE ! iはpやq以外の場合に限る Amat_new(i,p) = Amat(i,p) - sintheta * (Amat(i,q) + beta * Amat(i,p)) Amat_new(i,q) = Amat(i,q) + sintheta * (Amat(i,p) - beta * Amat(i,q)) Amat_new(p,i) = Amat_new(i,p) Amat_new(q,i) = Amat_new(i,q) END DO ! pq, qp成分はGivens回転の結果ゼロになるようにthetaを決めているのでゼロを代入 Amat_new(p,q) = 0.0d0 Amat_new(q,p) = 0.0d0 ! Givens回転をかけて直交行列の積を更新 Umat_new(1:N,1:N) = Umat(1:N,1:N) Umat_new(1:N,p) = Umat(1:N,p) - sintheta * ( Umat(1:N,q) + beta * Umat(1:N,p)) Umat_new(1:N,q) = Umat(1:N,q) + sintheta * ( Umat(1:N,p) - beta * Umat(1:N,q)) ! 次の反復操作のためにAmatとUmatの値を更新する Amat(1:N,1:N) = Amat_new(1:N,1:N) Umat(1:N,1:N) = Umat_new(1:N,1:N) ! この回の反復で得られたAmatの値を出力 DO i = 1, N WRITE(*,'(10ES15.5)') (Amat(i,j), j = 1, N) END DO ! 反復のカウンタを1増やす iter = iter + 1 ! 1000回反復しても収束しなければエラーメッセージを出して終了する IF( iter .GT. 1000) THEN WRITE(0,*) "did not converge after 1000 iterations" STOP END IF END DO ! 反復が終了した後の処理 ! 反復の後にAmatには対角行列が入っている。結果の出力 PRINT *, "Diagonalized matrix D" DO i = 1, N WRITE(*,'(10ES15.5)') (Amat(i,j), j = 1, N) END DO ! 直交行列Uの出力 PRINT *, "Orthogonal matrix U" DO i = 1, N WRITE(*,'(10ES15.5)') (Umat(i,j), j = 1, N) END DO ! Uが直交行列になっているか出力してチェック UTUmat(1:N,1:N) = MATMUL( TRANSPOSE(Umat(1:N,1:N)), Umat(1:N,1:N)) PRINT *, "U^T U matrix" DO i = 1, N WRITE(*,'(10ES15.5)') (UTUmat(i,j), j = 1, N) END DO ! U D U^{-1}がもとの行列Aになっているか出力してチェック UDUmat(1:N,1:N) = MATMUL( Umat(1:N,1:N), MATMUL(Amat(1:N,1:N), TRANSPOSE(Umat(1:N,1:N)))) PRINT *, "U^T D U matrix" DO i = 1, N WRITE(*,'(10ES15.5)') (UDUmat(i,j), j = 1, N) END DO RETURN END PROGRAM EXERCISE36