Class | Matrix_Calc |
In: |
matrix_calc.f90
|
行列計算を主に行うルーチン集. 固有値計算や逆行列計算, 連立方程式の求解などを行う. 一部にベクトル処理も入る.
Subroutine : | |||
a(size(b),size(b)) : | real, intent(inout)
| ||
b(:) : | real, intent(inout)
| ||
eps : | real, intent(in)
| ||
x(size(b)) : | real, intent(inout)
|
ガウスザイデル法による連立 1 次方程式ソルバ
subroutine Gau_Sei(a, b, eps, x) ! ガウスザイデル法による連立 1 次方程式ソルバ implicit none real, intent(inout) :: b(:) ! ax=b のベクトル real, intent(inout) :: a(size(b),size(b)) ! 係数行列 (第 1 要素が行列の行成分を表す) real, intent(in) :: eps ! 収束条件 real, intent(inout) :: x(size(b)) ! 解く解 integer :: i, j ! イテレーション用添字 real :: xn ! 更新した x(i) のテンプ領域 real :: err, err_max ! 誤差 integer :: nx !-- 初期値は 0.0 からスタートする --- x=0.0 nx=size(b) !-- ピボッティング call Pivot_part( a, b ) !-- 以下, while を使用するため, 1 回目のイテレートは単独で行う --- err_max=0.0 do i=1,nx xn=0.0 do j=1,nx if(j/=i)then xn=xn+a(i,j)*x(j) end if end do xn=(b(i)-xn)/a(i,i) err=errata(x(i),xn,1) if(err_max<=err)then err_max=err end if x(i)=xn end do !-- 以下より, 収束条件を満たすまでループする --- do while(err_max>=eps) err_max=0.0 do i=1,nx xn=0.0 do j=1,nx if(j/=i)then xn=xn+a(i,j)*x(j) end if end do xn=(b(i)-xn)/a(i,i) err=errata(x(i),xn,1) if(err_max<=err)then err_max=err end if x(i)=xn end do end do end subroutine Gau_Sei
Subroutine : | |||
a(:,:) : | real, intent(in)
| ||
b(size(a,1),size(a,2)) : | real, intent(inout)
| ||
p(size(a,1),size(a,2)) : | real, intent(inout), optional
| ||
sym_opt : | logical, intent(in), optional
|
ハウスホルダー変換を用いて, 任意実非対称行列をヘッセンベルグ行列に変換する. ただし, このルーチンでは計算速度を考慮して, 非対角の上三角成分は計算しない. 対称行列の場合, 上三角成分は下三角成分と同じになるので, 上三角成分の値を得たい場合は, b の下三角成分を求めればよい. 参考は川上 (2007) の 4.5.2.
subroutine Householder( a, b, p, sym_opt ) ! ハウスホルダー変換を用いて, 任意実非対称行列をヘッセンベルグ行列に変換する. ! ただし, このルーチンでは計算速度を考慮して, 非対角の上三角成分は計算しない. ! 対称行列の場合, 上三角成分は下三角成分と同じになるので, ! 上三角成分の値を得たい場合は, b の下三角成分を求めればよい. ! 参考は川上 (2007) の 4.5.2. implicit none real, intent(in) :: a(:,:) ! 任意の実非対称行列 real, intent(inout) :: b(size(a,1),size(a,2)) ! ヘッセンベルグ行列 real, intent(inout), optional :: p(size(a,1),size(a,2)) ! 変換に用いた直交行列 ! これを返す場合, 専用に計算処理を行うので, ! 計算時間が遅くなることに注意. logical, intent(in), optional :: sym_opt ! a が対称行列かどうか. ! default = .false. (対称ではない) ! .true. の場合, b の対称成分をそのまま返す. ! 計算速度を考慮して通常は対称成分を計算しない. integer :: i, j, k, n real :: coe, s, t, tmp real, allocatable, dimension(:) :: w, pconv, qconv real, dimension(size(a,1),size(a,2)) :: tmp_p, tmp_pp logical :: sym_link, pconv_flag n=size(a,1) if(present(sym_opt))then sym_link=sym_opt else sym_link=.false. end if if(present(p))then pconv_flag=.true. p=0.0 else pconv_flag=.false. end if if(size(a,1)/=size(a,2))then write(*,*) "*** ERROR (Householder) ***" write(*,*) "Input array must be square matrix. STOP." stop end if b(:,:)=a(:,:) allocate(w(n)) allocate(pconv(n)) allocate(qconv(n)) if(sym_link.eqv..true.)then do j=1,n-2 pconv=0.0 qconv=0.0 w=0.0 tmp_p=0.0 !-- 変換行列の作成 t=sum_sq( b(j+1:n,j) ) if(b(j+1,j)>0.0)then s=sqrt(t) else s=-sqrt(t) end if w(1:j)=0.0 w(j+1)=b(j+1,j)+s coe=1.0/(t+b(j+1,j)*s) ! この係数は川上 (2007) の c ではなく, 2c である. w(j+2:n)=b(j+2:n,j) do k=j+1,n do i=j+1,n pconv(k)=pconv(k)+b(k,i)*w(i) end do pconv(k)=coe*pconv(k) end do do k=j+1,n do i=j+1,n qconv(k)=qconv(k)+pconv(i)*w(i) end do qconv(k)=pconv(k)-0.5*coe*qconv(k)*w(k) end do !-- 変換行列による変換 b(j+1,j)=-s b(j+2:n,j)=0.0 do k=j+1,n do i=j+1,n b(i,k)=b(i,k)-qconv(i)*w(k)-w(i)*qconv(k) end do end do if(pconv_flag.eqv..true.)then do k=j+1,n ! w は j+1 以上の成分のみ非ゼロ. do i=j+1,n tmp_p(i,k)=-coe*w(i)*w(k) end do end do do k=1,n tmp_p(k,k)=1.0+tmp_p(k,k) end do if(j==1)then p(1:n,1:n)=tmp_p(1:n,1:n) else call mat_dot( p, tmp_p, tmp_pp ) p(1:n,1:n)=tmp_pp(1:n,1:n) end if end if end do do k=1,n do j=k,n b(k,j)=b(j,k) end do end do else ! 非対称行列の場合 (まだ対応していない.) do j=1,n-2 pconv=0.0 qconv=0.0 w=0.0 tmp_p=0.0 !-- 変換行列の作成 t=sum_sq( b(j+1:n,j) ) if(b(j+1,j)>0.0)then s=sqrt(t) else s=-sqrt(t) end if w(1:j)=0.0 w(j+1)=b(j+1,j)+s coe=1.0/(t+b(j+1,j)*s) ! この係数は川上 (2007) の c ではなく, 2c である. w(j+2:n)=b(j+2:n,j) do k=j+1,n do i=j+1,n pconv(k)=pconv(k)+b(k,i)*w(i) end do pconv(k)=coe*pconv(k) end do do k=j+1,n do i=j+1,n qconv(k)=qconv(k)+pconv(i)*w(i) end do qconv(k)=pconv(k)-0.5*coe*qconv(k)*w(k) end do !-- 変換行列による変換 b(j+1,j)=-s b(j+2:n,j)=0.0 do k=j+1,n do i=j+1,n b(i,k)=b(i,k)-qconv(i)*w(k)-w(i)*qconv(k) end do end do if(pconv_flag.eqv..true.)then do k=j+1,n ! w は j+1 以上の成分のみ非ゼロ. do i=j+1,n tmp_p(i,k)=-coe*w(i)*w(k) end do end do do k=1,n tmp_p(k,k)=1.0+tmp_p(k,k) end do if(j==1)then p(1:n,1:n)=tmp_p(1:n,1:n) else call mat_dot( p, tmp_p, tmp_pp ) p(1:n,1:n)=tmp_pp(1:n,1:n) end if end if end do end if end subroutine Householder
Subroutine : | |||
a(size(b),size(b)) : | real, intent(inout)
| ||
b(:) : | real, intent(inout)
| ||
eps : | real, intent(in)
| ||
x(size(b)) : | real, intent(inout)
|
ヤコビ法による連立 1 次方程式ソルバ
subroutine Jacobi_algebra(a, b, eps, x) ! ヤコビ法による連立 1 次方程式ソルバ implicit none real, intent(inout) :: b(:) ! ax=b のベクトル real, intent(inout) :: a(size(b),size(b)) ! 係数行列 (第 1 要素が行列の行を表す) real, intent(in) :: eps ! 収束条件 real, intent(inout) :: x(size(b)) ! 解く解 real :: y(size(b)) ! ヤコビ法で使用する一時格納用配列, この配列で一斉更新する.(xn の代わり) integer :: i, j ! イテレーション用添字 real :: err, err_max ! 誤差 integer :: nx nx=size(b) !-- 初期値は 0,0 からスタートする --- x=0.0 y=0.0 !-- ピボッティング call Pivot_part( a, b ) !-- 以下, 実際のソルバ(while を使用するため, 1 回目のイテレートは単独で行う) --- err_max=0.0 do i=1,nx y(i)=0.0 do j=1,nx if(j/=i)then y(i)=y(i)+a(i,j)*x(j) end if end do y(i)=(b(i)-y(i))/a(i,i) err=errata(x(i),y(i),1) if(err_max<=err)then err_max=err end if end do do i=1,nx ! データの一斉更新 x(i)=y(i) end do !-- 以下より, 収束条件を満たすまでループする --- do while(err_max>=eps) err_max=0.0 do i=1,nx y(i)=0.0 do j=1,nx if(j/=i)then y(i)=y(i)+a(i,j)*x(j) end if end do y(i)=(b(i)-y(i))/a(i,i) err=errata(x(i),y(i),1) if(err_max<=err)then err_max=err end if end do do i=1,nx ! データの一斉更新 x(i)=y(i) end do end do end subroutine Jacobi_algebra
Subroutine : | |||
a(:,:) : | real, intent(in)
| ||
lambda(size(a,1)) : | real, intent(inout)
| ||
eps : | real, intent(in), optional
|
Jacobi 法を用いて固有値を計算するルーチン. 本ルーチンでは, 計算対象となる行列は n 次の実対称行列でなければならない.
subroutine Jacobi_eigen( a, lambda, eps ) ! Jacobi 法を用いて固有値を計算するルーチン. ! 本ルーチンでは, 計算対象となる行列は n 次の実対称行列でなければならない. implicit none real, intent(in) :: a(:,:) ! 固有値を計算する実対称行列. ! 第 1 要素が行, 第 2 要素が列成分をそれぞれ表す. real, intent(inout) :: lambda(size(a,1)) ! 各固有値 real, intent(in), optional :: eps ! 反復計算の収束条件 [絶対誤差] ! デフォルトでは, 1.0e-6 integer :: i, j, k, n, m, l real :: tmp_a(size(a,1),size(a,2)), new_a(size(a,1),size(a,2)) real :: error, err_max, tan2, cos2, sin2 n=size(a,1) if(present(eps))then error=eps else error=1.0e-6 end if !-- intent(in) 属性なので, tmp へ入れ込む do j=1,n do i=1,n tmp_a(i,j)=a(i,j) new_a(i,j)=a(i,j) end do end do !-- 実際に計算させる. err_max=eps do while (err_max>=eps) err_max=0.0 do j=1,n-1 ! 各行について順に操作 do i=j+1,n ! 上三角のみ行えば, 対角成分は一意に求められる. if(tmp_a(i,j)/=0.0)then ! 以下で, 井桁更新用の係数を計算 if(tmp_a(i,i)/=tmp_a(j,j))then ! 対角成分の重解を処理 tan2=2.0*tmp_a(i,j)/(tmp_a(i,i)-tmp_a(j,j)) cos2=sqrt(0.5*(1.0+1.0/sqrt(1.0+tan2*tan2))) if(tan2>=0.0)then sin2=sqrt(0.5*(1.0-1.0/sqrt(1.0+tan2*tan2))) else sin2=-sqrt(0.5*(1.0-1.0/sqrt(1.0+tan2*tan2))) end if else ! この場合, tan2=\infty なので, cos2=sqrt(0.5) if(tmp_a(i,j)>=0.0)then ! 係数が正なら, \pi / 4 sin2=sqrt(0.5) else sin2=-sqrt(0.5) end if end if ! 以降で実際に井桁を更新 new_a(i,j)=0.0 ! これら 2 つは先に計算しておく. new_a(j,i)=0.0 ! (あとで同じ計算を 2 回行うのを回避するため) do k=1,n if(k/=i.and.k/=j)then new_a(i,k)=tmp_a(i,k)*cos2+tmp_a(j,k)*sin2 new_a(j,k)=-tmp_a(i,k)*sin2+tmp_a(j,k)*cos2 new_a(k,i)=tmp_a(k,i)*cos2+tmp_a(k,j)*sin2 new_a(k,j)=-tmp_a(k,i)*sin2+tmp_a(k,j)*cos2 else if(k==i)then new_a(i,k)=tmp_a(i,i)*cos2*cos2+tmp_a(j,j)*sin2*sin2 +2.0*tmp_a(i,j)*sin2*cos2 else new_a(j,j)=tmp_a(i,i)*sin2*sin2+tmp_a(j,j)*cos2*cos2 -2.0*tmp_a(i,j)*sin2*cos2 end if end if end do ! 以下で, new_a -> tmp_a に戻し, 同時に最大誤差を求める. do m=1,n do l=1,n error=abs(new_a(l,m)-tmp_a(l,m)) tmp_a(l,m)=new_a(l,m) if(error>err_max)then ! 誤差の最大値を求める. err_max=error end if end do end do end if end do end do end do do i=1,n lambda(i)=tmp_a(i,i) end do end subroutine Jacobi_eigen
Subroutine : | |||
a(size(b),size(b)) : | real, intent(inout)
| ||
b(:) : | real, intent(inout)
| ||
x(size(b)) : | real, intent(inout)
| ||
itermax : | integer, intent(in)
|
— LU 分解を計算するサブルーチン —
subroutine LU_devs( a, b, x, itermax ) !-- LU 分解を計算するサブルーチン --- implicit none real, intent(inout) :: b(:) ! 右辺のベクトル real, intent(inout) :: a(size(b),size(b)) ! 係数行列 (第 1 要素が行を表現) real, intent(inout) :: x(size(b)) ! 解 integer, intent(in) :: itermax ! 反復の回数 real :: d(size(b),size(b)), r(size(b)), y(size(b)) integer :: ip(size(b)) real :: scale(size(b)), dx(size(b)) real :: s, t, pivot, xnorm, dxnorm real :: s1, s2, s3, s4, s5, t0, t1, t3, t4, eps integer :: iter, nmax integer :: p, itemp, i, j, k nmax=size(b) !-- 反復改良での精度の設定 --- t4=6.0 !-- 配列 x(i) の初期化 --- do i=1,nmax x(i)=0.0 end do do i=1,nmax do j=1,nmax d(i,j)=a(i,j) end do ip(i)=i !-- 最大値を計算するループ --- s=d(i,1) do j=2,nmax if(d(i,j).gt.s)then s=d(i,j) end if end do scale(i)=1.0/s end do do k=1,nmax t=d(ip(k),k)*scale(ip(k)) p=k do i=k,nmax t0=d(ip(i),k)*scale(ip(i)) if(t0.gt.t)then t=t0 p=i end if end do !-- ip(p) と ip(k) の入れ替え --- if(p.ne.k)then itemp=ip(p) ip(p)=ip(k) ip(k)=itemp end if pivot=d(ip(k),k) do i=k+1,nmax d(ip(i),k)=d(ip(i),k)/pivot do j=k+1,nmax d(ip(i),j)=d(ip(i),j)-d(ip(i),k)*d(ip(k),j) end do end do if(k.ge.nmax-1)then exit end if end do !-- 前進消去 --- y(1)=b(ip(1)) do i=2,nmax s1=0.0 do j=1,i-1 s1=s1+d(ip(i),j)*y(j) end do y(i)=b(ip(i))-s1 end do !-- 後退代入 --- x(nmax)=y(nmax)/d(ip(nmax),nmax) do i=nmax-1,1,-1 s2=0.0 do j=i+1,nmax s2=s2+d(ip(i),j)*y(j) end do x(i)=(y(i)-s2)/d(ip(i),i) end do t1=x(1) xnorm=x(1) do i=2,nmax if(x(i).gt.t1)then t1=x(i) xnorm=x(i) end if end do !-- 反復改良 --- eps=10**(-t4) ! 標準精度を 10 進数の t4 桁とする do iter=1,itermax if(xnorm==0.0)then exit end if !-- 残差の計算 --- do i=1,nmax s3=0.0 do j=1,nmax s3=s3+a(i,j)*x(j) end do r(i)=b(i)-s3 end do !-- 前進消去 --- y(1)=r(ip(1)) do i=2,nmax s4=0.0 do j=1,i-1 s4=s4+d(ip(i),j)*y(j) end do y(i)=r(ip(i))-s4 end do !-- 後退代入 --- dx(nmax)=y(nmax)/d(ip(nmax),nmax) do i=nmax-1,1,-1 s5=0.0 do j=i+1,nmax s5=s5+d(ip(i),j)*y(j) end do dx(i)=(y(i)-s5)/d(ip(i),i) end do do i=1,nmax x(i)=x(i)+dx(i) end do t3=dx(1) dxnorm=dx(1) do i=1,nmax if(dx(i).gt.t3)then t3=dx(i) dxnorm=dx(i) end if end do if(dxnorm/xnorm.le.eps)then exit end if end do end subroutine LU_devs
Subroutine : | |||
a(:,:) : | real, intent(in)
| ||
lambda(size(a,1)) : | real, intent(inout)
| ||
p(size(a,1),size(a,2)) : | real, intent(inout), optional
| ||
eps : | real, intent(in), optional
| ||
sym_opt : | logical, intent(in), optional
|
QR 分解法を用いて行列 a の全固有値を計算するルーチン. 参考文献は川上 (2007) 4.8. 規格化された固有ベクトルも求めるようにしているが, 値が正確かどうかは未確認. 固有値に関しては例題から確認済み. 計算順序としては,
subroutine QR_method( a, lambda, p, eps, sym_opt ) ! QR 分解法を用いて行列 a の全固有値を計算するルーチン. ! 参考文献は川上 (2007) 4.8. ! 規格化された固有ベクトルも求めるようにしているが, 値が正確かどうかは未確認. ! 固有値に関しては例題から確認済み. ! 計算順序としては, ! 1. Householder 変換 ! 2. 直交行列の計算 ! 3. QR 分解 ! 4. QR 逆積 implicit none real, intent(in) :: a(:,:) ! 固有値を求める行列 [第 1 要素が行] real, intent(inout) :: lambda(size(a,1)) ! a の固有値 real, intent(inout), optional :: p(size(a,1),size(a,2)) ! 固有ベクトル. ! これを返すには, Householder 変換時の処理時間が遅くなる. real, intent(in), optional :: eps ! 反復法の収束条件 default = 1.0e-5 logical, intent(in), optional :: sym_opt ! a が対称行列か. ! .true. = 対称行列, default = .false. integer :: i, j, k, n, counter real :: dia, s, t, err, err_max, tmp1, tmp2, mu real, dimension(size(a,1),size(a,2)) :: a0, q, r, b, tmp_p, tmp_pp logical :: sym_link, pconv_flag n=size(a,1) if(present(sym_opt))then sym_link=sym_opt else sym_link=.false. end if if(present(p))then pconv_flag=.true. p=0.0 else pconv_flag=.false. end if if(present(eps))then err_max=eps else err_max=1.0e-5 end if !-- 1. Householder 変換. if(pconv_flag.eqv..true.)then call Householder( a, a0, p=tmp_p, sym_opt=sym_link ) else call Householder( a, a0, sym_opt=sym_link ) end if counter=n do while(counter>1) err=err_max do while(err_max<=err) q=0.0 ! QR 分解の Q_k に相当 do k=1,n q(k,k)=1.0 ! 最初に単位行列として初期化 end do mu=a0(counter,counter) a0(counter,counter)=0.0 do k=1,counter-1 a0(k,k)=a0(k,k)-mu ! moving original point end do !-- QR 分解を行う. do k=1,counter-1 dia=sqrt(abs(a0(k,k))**2+abs(a0(k+1,k))**2) s=a0(k+1,k)/dia t=a0(k+1,k)/(dia+a0(k,k)) a0(k+1,k)=0.0 a0(k,k)=dia do j=k+1,counter tmp1=a0(k,j) tmp2=a0(k+1,j) a0(k,j)=tmp1+s*(tmp2-t*tmp1) a0(k+1,j)=tmp2-s*(tmp1+t*tmp2) end do do i=1,counter tmp1=q(i,k) tmp2=q(i,k+1) q(i,k)=tmp1+s*(tmp2-t*tmp1) q(i,k+1)=tmp2-s*(tmp1+t*tmp2) end do end do do j=1,counter do i=1,counter r(i,j)=a0(i,j) ! QR 分解で得られた上三角行列 R end do end do !-- QR の逆積を計算する. A=RQ で, R は上三角, Q はヘッセンベルグなので, !-- a_{i,j} = R_{i,k} * Q_{k,j} + R_{i,k+1} * Q_{k+1,j}, (k=i) !-- ただし, Q は転置で計算されるので, Q=qT とすると, Q_{k,j} = q_{j,k}. !-- a_{i,j} = R_{i,k} * q_{j,k} + R_{i,k+1} * q_{j,k+1}, (k=i) call mat_dot( r, q, a0 ) !-- mooving grid do i=1,counter a0(i,i)=a0(i,i)+mu end do !-- 収束判定条件 if(abs(err_max)>abs(a0(counter,counter-1)))then err=abs(a0(counter,counter-1)) end if if(pconv_flag.eqv..true.)then call mat_dot( q, tmp_p, tmp_pp ) tmp_p(1:n,1:n)=tmp_pp(1:n,1:n) end if end do lambda(counter)=a0(counter,counter) counter=counter-1 end do if(pconv_flag.eqv..true.)then p(1:n,1:n)=tmp_p(1:n,1:n) end if lambda(1)=a0(1,1) end subroutine QR_method
Subroutine : | |||
a(size(b),size(b)) : | real, intent(inout)
| ||
b(:) : | real, intent(inout)
| ||
eps : | real, intent(in)
| ||
accel : | real, intent(in)
| ||
x(size(b)) : | real, intent(inout)
|
ガウスザイデル法かつ, SOR で加速による連立 1 次方程式ソルバ
subroutine SOR_Gau_Sei(a, b, eps, accel, x) ! ガウスザイデル法かつ, SOR で加速による連立 1 次方程式ソルバ implicit none real, intent(inout) :: b(:) ! ax=b のベクトル real, intent(inout) :: a(size(b),size(b)) ! 係数行列 (第 1 要素が行列の行を表す) real, intent(in) :: eps ! 収束条件 real, intent(in) :: accel ! 加速係数. ただし, 数学的に accel >= 2 では発散するので, ! この値以上が設定されるとエラーで止める. real, intent(inout) :: x(size(b)) ! 解く解 integer :: i, j ! イテレーション用添字 real :: xn ! 更新した x(i) のテンプ領域 real :: err, err_max ! 誤差 integer :: nx nx=size(b) !-- 加速パラメータの確認 if(accel>=2.0)then write(*,*) "***** ERROR *****" write(*,*) "accel parameter must be less than 2.0. STOP." stop end if !-- 初期値は 0.0 からスタートする --- x=0.0 !-- ピボッティング call Pivot_part( a, b ) !-- 以下, while を使用するため, 1 回目のイテレートは単独で行う --- err_max=0.0 do i=1,nx xn=0.0 do j=1,nx if(j/=i)then xn=xn+a(i,j)*x(j) end if end do xn=(b(i)-xn)/a(i,i) xn=x(i)+accel*(xn-x(i)) err=errata(x(i),xn,1) if(err_max<=err)then err_max=err end if x(i)=xn end do !-- 以下より, 収束条件を満たすまでループする --- do while(err_max>=eps) err_max=0.0 do i=1,nx xn=0.0 do j=1,nx if(j/=i)then xn=xn+a(i,j)*x(j) end if end do xn=(b(i)-xn)/a(i,i) xn=x(i)+accel*(xn-x(i)) err=errata(x(i),xn,1) if(err_max<=err)then err_max=err end if x(i)=xn end do end do end subroutine SOR_Gau_Sei
Subroutine : | |||
a(size(b),size(b)) : | real, intent(inout)
| ||
b(:) : | real, intent(inout)
| ||
eps : | real, intent(in)
| ||
accel : | real, intent(in)
| ||
x(size(b)) : | real, intent(inout)
|
ヤコビ法かつ SOR 加速による連立 1 次方程式ソルバ
subroutine SOR_Jacobi_algebra(a, b, eps, accel, x) ! ヤコビ法かつ SOR 加速による連立 1 次方程式ソルバ implicit none real, intent(inout) :: b(:) ! ax=b のベクトル real, intent(inout) :: a(size(b),size(b)) ! 係数行列 (第 1 要素が行列の行を表す) real, intent(in) :: eps ! 収束条件 real, intent(in) :: accel ! 加速係数. ただし, 数学的に accel >= 2 では発散するので, ! この値以上が設定されるとエラーで止める. real, intent(inout) :: x(size(b)) ! 解く解 real :: y(size(b)) ! ヤコビ法で使用する一時格納用配列, この配列で一斉更新する.(xn の代わり) integer :: i, j ! イテレーション用添字 real :: err, err_max ! 誤差 integer :: nx nx=size(b) !-- 加速パラメータの確認 if(accel>=2.0)then write(*,*) "***** ERROR *****" write(*,*) "accel parameter must be less than 2.0. STOP." stop end if !-- 初期値は 0,0 からスタートする --- x=0.0 y=0.0 !-- ピボッティング call Pivot_part( a, b ) !-- 以下, 実際のソルバ(while を使用するため, 1 回目のイテレートは単独で行う) --- err_max=0.0 do i=1,nx y(i)=0.0 do j=1,nx if(j/=i)then y(i)=y(i)+a(i,j)*x(j) end if end do y(i)=(b(i)-y(i))/a(i,i) err=errata(x(i),y(i),1) if(err_max<=err)then err_max=err end if end do do i=1,nx ! データの一斉更新 x(i)=x(i)+accel*(y(i)-x(i)) end do !-- 以下より, 収束条件を満たすまでループする --- do while(err_max>=eps) err_max=0.0 do i=1,nx y(i)=0.0 do j=1,nx if(j/=i)then y(i)=y(i)+a(i,j)*x(j) end if end do y(i)=(b(i)-y(i))/a(i,i) y(i)=x(i)+accel*(y(i)-x(i)) err=errata(x(i),y(i),1) if(err_max<=err)then err_max=err end if end do do i=1,nx ! データの一斉更新 x(i)=y(i) end do end do end subroutine SOR_Jacobi_algebra
Function : | recursive | ||
res : | integer | ||
a : | integer, dimension(2,2), intent(in)
|
2x2 の正方行列の行列式を計算する関数(整数版)
Alias for determ_2d_i
Function : | recursive | ||
res : | real | ||
a : | real, dimension(2,2), intent(in)
|
2x2 の正方行列の行列式を計算する関数(実数版)
Alias for determ_2d_f
Function : | recursive | ||
res : | real | ||
a : | real, dimension(2,2), intent(in)
|
2x2 の正方行列の行列式を計算する関数(実数版)
recursive function determ_2d_f( a ) result(res) ! 2x2 の正方行列の行列式を計算する関数(実数版) implicit none real, dimension(2,2), intent(in) :: a ! 2x2 の正方行列 real :: res res=a(1,1)*a(2,2)-a(1,2)*a(2,1) return end function
Function : | recursive | ||
res : | integer | ||
a : | integer, dimension(2,2), intent(in)
|
2x2 の正方行列の行列式を計算する関数(整数版)
recursive function determ_2d_i( a ) result(res) ! 2x2 の正方行列の行列式を計算する関数(整数版) implicit none integer, dimension(2,2), intent(in) :: a ! 2x2 の正方行列 integer :: res res=a(1,1)*a(2,2)-a(1,2)*a(2,1) return end function
Subroutine : | |||
a(:,:) : | real, intent(in)
| ||
lambda(size(a,1)) : | real, intent(inout)
| ||
eigen_vec(size(a,1),size(a,2)) : | real, intent(inout)
| ||
eps : | real, intent(in), optional
|
べき乗法を用いて, a についての全ての固有値とそれに伴う固有ベクトルを 計算するルーチン.
subroutine eigen_power_all( a, lambda, eigen_vec, eps ) ! べき乗法を用いて, a についての全ての固有値とそれに伴う固有ベクトルを ! 計算するルーチン. implicit none real, intent(in) :: a(:,:) ! 固有値を求める行列 [第 1 要素が行] real, intent(inout) :: lambda(size(a,1)) ! a の固有値 real, intent(inout) :: eigen_vec(size(a,1),size(a,2)) ! lambda(i) に対応する固有ベクトル, 第一要素が固有値 lambda(i) に対応している. real, intent(in), optional :: eps ! 反復法の収束条件 integer :: n n=size(a,1) end subroutine
Subroutine : | |||
a(size(eigenvec),size(eigenvec)) : | real, intent(in)
| ||
eps : | real, intent(in)
| ||
eigenval : | real, intent(inout)
| ||
eigenvec(:) : | real, intent(inout)
|
べき乗法を用いて行列の最大固有値とその固有値に対応する固有ベクトルを求める.
subroutine eigenvalue_power( a, eps, eigenval, eigenvec ) ! べき乗法を用いて行列の最大固有値とその固有値に対応する固有ベクトルを求める. implicit none real, intent(inout) :: eigenvec(:) ! 固有ベクトル real, intent(in) :: a(size(eigenvec),size(eigenvec)) ! 固有値を求める行列 (第 1 要素が行を表す) real, intent(in) :: eps ! 収束判定条件 real, intent(inout) :: eigenval ! 行列 a の最大固有値 integer :: i, j real, dimension(size(eigenvec)) :: x, y real :: tmp1, tmp2, err_max integer :: nx nx=size(eigenvec) do i=1,nx x(i)=1.0 ! 反復法の初期値として非ゼロのベクトルを定義 end do tmp1=sqrt(vec_dot( x, x )) ! 初期ベクトルのノルムを計算 do i=1,nx x(i)=x(i)/tmp1 ! 初期ベクトルを規格化 end do err_max=eps ! while 文に入れるための処理 !-- 反復法開始 do while(err_max>=eps) err_max=0.0 do i=1,nx y(i)=0.0 ! 配列の初期化 do j=1,nx y(i)=y(i)+a(i,j)*x(j) end do end do tmp1=sqrt(vec_dot( y, y )) ! ベクトルのノルムを計算 do i=1,nx x(i)=y(i)/tmp1 ! x(i) の更新 end do !-- 固有値計算 tmp2=vec_dot( x, y ) ! 上で計算した Ay に y^t (つまり, 固有ベクトルの転置) をかける. err_max=errata( eigenval, tmp2, 1 ) ! 過去の x(i) と更新した y(i) の誤差比較 eigenval=tmp2 end do !-- 反復法終了 do i=1,nx eigenvec(i)=x(i) ! 固有ベクトルの変数へ代入 end do end subroutine eigenvalue_power
Subroutine : | |||
c(size(d),size(d)) : | real, intent(in)
| ||
d(:) : | real, intent(in) | ||
x(size(d)) : | real, intent(inout) |
部分ピボット付きガウスの消去法
subroutine gausss( c, d, x ) ! 部分ピボット付きガウスの消去法 implicit none real, intent(in) :: d(:) real, intent(in) :: c(size(d),size(d)) ! 係数行列 (第 1 要素が行を表す) real, intent(inout) :: x(size(d)) real :: b(size(d)) real :: a(size(d),size(d)) ! 係数行列 (第 1 要素が行を表す) real :: s, pivotb real :: pivot(size(d)+1) integer :: piv, i, j, k, nmax nmax=size(b) do k=1,nmax do j=1,nmax a(j,k)=c(j,k) end do b(k)=d(k) end do !-- 前進消去 --- !-- A(I,J) の前進消去 --- do k=1,nmax-1 !-- PIVOT の選択 --- !-- まず, I 成分の最大値を決定する --- piv=k do i=k+1,nmax if(abs(a(i,k)).gt.abs(a(piv,k)))then piv=i end if end do !-- ここまでで, 最大値が決定される --- !-- 以下で, 最大値をとる成分の行を入れ替える --- do j=k,nmax pivot(j)=a(k,j) a(k,j)=a(piv,j) a(piv,j)=pivot(j) end do pivotb=b(k) b(k)=b(piv) b(piv)=pivotb !-- PIVOT 処理ここまで --- do i=k+1,nmax a(k,i)=a(k,i)/a(k,k) end do b(k)=b(k)/a(k,k) a(k,k)=1.0 do j=k+1,nmax do i=k+1,nmax a(j,i)=a(j,i)-a(k,i)*a(j,k) end do b(j)=b(j)-b(k)*a(j,k) a(j,k)=0.0 end do end do b(nmax)=b(nmax)/a(nmax,nmax) a(nmax,nmax)=1.0 !-- X(I) の後進代入 x(nmax)=b(nmax) do i=nmax-1,1,-1 s=b(i) do j=i+1,nmax s=s-a(i,j)*x(j) end do x(i)=s end do end subroutine gausss
Subroutine : | |||
ax(:,:) : | real, intent(in)
| ||
xx(size(ax,1),size(ax,2)) : | real, intent(inout)
|
ガウスの消去法を拡張して, 逆行列を計算する. 具体的には, 右辺の b に単位ベクトルを 1 つずつ代入し,消去した結果のベクトルを並べて 逆行列にする.
subroutine invert_mat( ax, xx ) ! ガウスの消去法を拡張して, 逆行列を計算する. ! 具体的には, 右辺の b に単位ベクトルを 1 つずつ代入し,消去した結果のベクトルを並べて ! 逆行列にする. implicit none real, intent(in) :: ax(:,:) ! 逆行列を求める行列 real, intent(inout) :: xx(size(ax,1),size(ax,2)) ! 求めた逆行列 integer :: i, j, k real :: c(size(ax,1),size(ax,2)) ! 単位行列 real :: d(size(ax,1),size(ax,2)) ! a(i,j) をガウスルーチンに渡すと, 結果が変化するのでこれに一時退避 integer :: nx nx=size(ax,1) c=0.0 do i=1,nx c(i,i)=1.0 end do do i=1,nx do j=1,nx do k=1,nx d(k,j)=ax(k,j) ! ダミー変数へ代入 end do end do call gausss( d, c(:,i), xx(:,i) ) ! 配列第 2 成分が列を表すので, この順番で OK. end do end subroutine invert_mat
Subroutine : | |
a : | real, dimension(:,:), intent(in) |
b : | real, dimension(:,:), intent(in) |
c : | real, dimension(:,:), intent(inout) |
行列同士の積を計算する (a * b = c).
subroutine mat_dot( a, b, c ) ! 行列同士の積を計算する (a * b = c). implicit none real, dimension(:,:), intent(in) :: a, b real, dimension(:,:), intent(inout) :: c integer :: i, j, k, l, m, n m=size(a,1) n=size(a,2) if(n/=size(b,1).or.m/=size(b,2).or.m/=size(c,1).or.m/=size(c,2))then write(*,*) "ERROR (mat_dot) : when a * b = c, a=mxn -> b=nxm, c=mxm" write(*,*) "STOP." stop end if do j=1,m do i=1,m c(i,j)=0.0 do k=1,n c(i,j)=c(i,j)+a(i,k)*b(k,j) end do end do end do end subroutine mat_dot
Subroutine : | |||
u(:,:) : | real, intent(in)
| ||
v(size(u,1),size(u,2)) : | real, intent(inout)
|
シュミットの直交化法を用いて, nx 次元ベクトルを正規直交化する. 引数の配列は第 1 要素がベクトルの各成分, 第 2 成分がベクトル群の 1 ベクトルを表す. つまり, u(i,j) は j 番目ベクトルの第 i 成分ということを表す. 行列で表現するなら, 縦ベクトルを横に並べた形と等しい. 本ルーチンはシュミットの直交化を真面目に計算しているので, 計算速度は比較的遅いらしい.
subroutine schumit_norm( u, v ) ! シュミットの直交化法を用いて, nx 次元ベクトルを正規直交化する. ! 引数の配列は第 1 要素がベクトルの各成分, 第 2 成分がベクトル群の 1 ベクトルを表す. ! つまり, u(i,j) は j 番目ベクトルの第 i 成分ということを表す. ! 行列で表現するなら, 縦ベクトルを横に並べた形と等しい. ! 本ルーチンはシュミットの直交化を真面目に計算しているので, ! 計算速度は比較的遅いらしい. implicit none real, intent(in) :: u(:,:) ! 直交化する前のベクトル real, intent(inout) :: v(size(u,1),size(u,2)) ! 直交化後のベクトル integer :: i, j, k, nx, ny real :: tmpn(size(u,2)) ! 正規化を行うときの, ノルムの値を格納する. real :: tmps(size(u,1)) ! 総和演算の際の一時格納に使用. nx=size(u,1) ny=size(u,2) tmpn(1)=sqrt(vec_dot( u(:,1), u(:,1) )) do i=1,nx ! 1 本目のベクトルを基底の基準とする. v(i,1)=u(i,1)/tmpn(1) end do do j=2,ny do i=1,nx tmps(i)=0.0 end do do k=1,j-1 do i=1,nx tmps(i)=tmps(i)+vec_dot( u(:,j), v(:,k) )*u(i,j) end do end do do i=1,nx v(i,j)=u(i,j)-tmps(i) end do ! 以下で正規化を行う tmpn(j)=sqrt(vec_dot( v(:,j), v(:,j) )) do i=1,nx v(i,j)=v(i,j)/tmpn(j) end do end do end subroutine schumit_norm
Subroutine : | |||
a : | real, intent(inout), dimension(:,:)
|
行列成分の転置を返す(実数版)
subroutine trans_mat_f( a ) ! 行列成分の転置を返す(実数版) implicit none real, intent(inout), dimension(:,:) :: a ! 入力行列 integer :: i, j, nx real :: tmp nx=size(a,1) do j=1,nx do i=1,nx if(i<j)then tmp=a(j,i) a(j,i)=a(i,j) a(i,j)=tmp end if end do end do end subroutine
Subroutine : | |||
a : | integer, intent(inout), dimension(:,:)
|
行列成分の転置を返す
subroutine trans_mat_i( a ) ! 行列成分の転置を返す implicit none integer, intent(inout), dimension(:,:) :: a ! 入力行列 integer :: i, j, tmp, nx nx=size(a,1) do j=1,nx do i=1,nx if(i<j)then tmp=a(j,i) a(j,i)=a(i,j) a(i,j)=tmp end if end do end do end subroutine