Class | Algebra |
In: |
algebra.f90
|
代数演算を主に行うモジュール 固有値計算, 行列演算, ベクトル演算を主に行う.
Subroutine : | |||
x(:,:) : | real, intent(in)
| ||
y(size(x,1),size(x,2)) : | real, intent(in)
| ||
dis(size(x,1),size(x,2)) : | real, intent(inout)
|
2 次元ベクトルの絶対値を計算するルーチン
subroutine abst_2d(x,y,dis) ! 2 次元ベクトルの絶対値を計算するルーチン implicit none real, intent(in) :: x(:,:) ! x 方向のベクトル成分 real, intent(in) :: y(size(x,1),size(x,2)) ! y 方向のベクトル成分 real, intent(inout) :: dis(size(x,1),size(x,2)) ! 各点での絶対値ベクトル integer :: i, j, nx, ny nx=size(x,1) ny=size(x,2) do j=1,ny do i=1,nx dis(i,j)=sqrt(x(i,j)**2+y(i,j)**2) end do end do end subroutine abst_2d
Subroutine : | |||
x(:,:,:) : | real, intent(in)
| ||
y(size(x,1),size(x,2),size(x,3)) : | real, intent(in)
| ||
z(size(x,1),size(x,2),size(x,3)) : | real, intent(in)
| ||
dis(size(x,1),size(x,2),size(x,3)) : | real, intent(inout)
|
3 次元ベクトルの絶対値を計算するルーチン 配列数を調整することにより, 2 次元での計算も可能.
subroutine abst_3d(x,y,z,dis) ! 3 次元ベクトルの絶対値を計算するルーチン ! 配列数を調整することにより, 2 次元での計算も可能. implicit none real, intent(in) :: x(:,:,:) ! x 方向のベクトル成分 real, intent(in) :: y(size(x,1),size(x,2),size(x,3)) ! y 方向のベクトル成分 real, intent(in) :: z(size(x,1),size(x,2),size(x,3)) ! z 方向のベクトル成分 real, intent(inout) :: dis(size(x,1),size(x,2),size(x,3)) ! 各点での絶対値ベクトル integer :: i, j, k, nx, ny, nz nx=size(x,1) ny=size(x,2) nz=size(x,3) !$omp parallel default(shared) !$omp do schedule(runtime) private(i,j,k) do k=1,nz do j=1,ny do i=1,nx dis(i,j,k)=sqrt(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2) end do end do end do !$omp end do !$omp end parallel end subroutine abst_3d
Subroutine : | |||
xp : | real, intent(in)
| ||
yp : | real, intent(in)
| ||
zp : | real, intent(in)
| ||
x(:) : | real, intent(in)
| ||
y(:) : | real, intent(in)
| ||
z(:) : | real, intent(in)
| ||
rad(size(x),size(y),size(z)) : | real, intent(inout)
|
ある位置からの距離を計算するルーチン 配列数を調整することにより, 2 次元での計算も可能.
subroutine calc_radius(xp,yp,zp,x,y,z,rad) ! ある位置からの距離を計算するルーチン ! 配列数を調整することにより, 2 次元での計算も可能. implicit none real, intent(in) :: xp ! 中心位置座標の x 成分 real, intent(in) :: yp ! 中心位置座標の y 成分 real, intent(in) :: zp ! 中心位置座標の z 成分 real, intent(in) :: x(:) ! x 方向の位置座標 real, intent(in) :: y(:) ! y 方向の位置座標 real, intent(in) :: z(:) ! z 方向の位置座標 real, intent(inout) :: rad(size(x),size(y),size(z)) ! 距離 integer :: i, j, k, nx, ny, nz nx=size(x) ny=size(y) nz=size(z) !$omp parallel default(shared) !$omp do schedule(runtime) private(i,j,k) do k=1,nz do j=1,ny do i=1,nx rad(i,j,k)=sqrt((x(i)-xp)**2+(y(j)-yp)**2+(z(k)-zp)**2) end do end do end do !$omp end do !$omp end parallel end subroutine calc_radius
Subroutine : | |||
x(:,:) : | real, intent(in)
| ||
y(size(x,1),size(x,2)) : | real, intent(in)
| ||
u(size(x,1),size(x,2)) : | real, intent(in)
| ||
v(size(x,1),size(x,2)) : | real, intent(in)
| ||
dot(size(x,1),size(x,2)) : | real, intent(inout)
| ||
undeff : | real, intent(in), optional
|
2ベクトルの内積計算ルーチン dot_prod の 2 次元版
subroutine dot_prod_2d(x,y,u,v,dot,undeff) ! 2ベクトルの内積計算ルーチン ! dot_prod の 2 次元版 implicit none real, intent(in) :: x(:,:) ! x 方向のベクトル成分 real, intent(in) :: y(size(x,1),size(x,2)) ! y 方向のベクトル成分 real, intent(in) :: u(size(x,1),size(x,2)) ! x 方向のベクトル成分 real, intent(in) :: v(size(x,1),size(x,2)) ! y 方向のベクトル成分 real, intent(inout) :: dot(size(x,1),size(x,2)) ! 内積 real, intent(in), optional :: undeff ! 未定義値 integer :: i, j, nx, ny nx=size(x,1) ny=size(x,2) if(present(undeff))then do j=1,ny do i=1,nx if(x(i,j)==undeff.or.u(i,j)==undeff.or. y(i,j)==undeff.or.v(i,j)==undeff)then dot(i,j)=undeff else dot(i,j)=x(i,j)*u(i,j)+y(i,j)*v(i,j) end if end do end do else do j=1,ny do i=1,nx dot(i,j)=x(i,j)*u(i,j)+y(i,j)*v(i,j) end do end do end if end subroutine dot_prod_2d
Subroutine : | |||
x(:,:,:) : | real, intent(in)
| ||
y(size(x,1),size(x,2),size(x,3)) : | real, intent(in)
| ||
z(size(x,1),size(x,2),size(x,3)) : | real, intent(in)
| ||
u(size(x,1),size(x,2),size(x,3)) : | real, intent(in)
| ||
v(size(x,1),size(x,2),size(x,3)) : | real, intent(in)
| ||
w(size(x,1),size(x,2),size(x,3)) : | real, intent(in)
| ||
dot(size(x,1),size(x,2),size(x,3)) : | real, intent(inout)
| ||
undeff : | real, intent(in), optional |
2ベクトルの内積計算ルーチン 配列を工夫すると, 2 次元での内積を計算することも可能
subroutine dot_prod_3d(x,y,z,u,v,w,dot,undeff) ! 2ベクトルの内積計算ルーチン ! 配列を工夫すると, 2 次元での内積を計算することも可能 implicit none real, intent(in) :: x(:,:,:) ! x 方向のベクトル成分 real, intent(in) :: y(size(x,1),size(x,2),size(x,3)) ! y 方向のベクトル成分 real, intent(in) :: z(size(x,1),size(x,2),size(x,3)) ! z 方向のベクトル成分 real, intent(in) :: u(size(x,1),size(x,2),size(x,3)) ! x 方向のベクトル成分 real, intent(in) :: v(size(x,1),size(x,2),size(x,3)) ! y 方向のベクトル成分 real, intent(in) :: w(size(x,1),size(x,2),size(x,3)) ! z 方向のベクトル成分 real, intent(inout) :: dot(size(x,1),size(x,2),size(x,3)) ! 内積 real, intent(in), optional :: undeff integer :: i, j, k, nx, ny, nz nx=size(x,1) ny=size(x,2) nz=size(x,3) if(present(undeff))then !$omp parallel default(shared) !$omp do schedule(dynamic) private(i,j,k) do k=1,nz do j=1,ny do i=1,nx if(x(i,j,k)==undeff.or.u(i,j,k)==undeff.or.y(i,j,k)==undeff.or. v(i,j,k)==undeff.or.z(i,j,k)==undeff.or.w(i,j,k)==undeff)then dot(i,j,k)=undeff else dot(i,j,k)=x(i,j,k)*u(i,j,k)+y(i,j,k)*v(i,j,k)+z(i,j,k)*w(i,j,k) end if end do end do end do !$omp end do !$omp end parallel else !$omp parallel default(shared) !$omp do schedule(runtime) private(i,j,k) do k=1,nz do j=1,ny do i=1,nx dot(i,j,k)=x(i,j,k)*u(i,j,k)+y(i,j,k)*v(i,j,k)+z(i,j,k)*w(i,j,k) end do end do end do !$omp end do !$omp end parallel end if end subroutine dot_prod_3d
Subroutine : | |||
x(:) : | real, intent(in)
| ||
y(size(x)) : | real, intent(in)
| ||
bot : | real, intent(in)
| ||
top : | real, intent(in)
| ||
res : | real, intent(inout)
| ||
undeff : | real, intent(in), optional |
1 次元台形積分 不等間隔でも計算可能であるが, 精度は保証しない.
subroutine rectangle_int( x, y, bot, top, res, undeff ) ! 1 次元台形積分 ! 不等間隔でも計算可能であるが, 精度は保証しない. implicit none real, intent(in) :: bot ! 積分区間左端 real, intent(in) :: top ! 積分区間右端 real, intent(in) :: x(:) ! 積分変数 real, intent(in) :: y(size(x)) ! 非積分関数 real, intent(inout) :: res ! 台形積分の積分値 real, intent(in), optional :: undeff integer :: i, nx, i_bot, i_top real :: y_bot, y_top nx=size(x) res=0.0 !-- bot < top でなければエラー if(bot>top)then write(*,*) "#### ERROR (algebra:rectangle_int) ####" write(*,*) "integrated interval must be bot < top. STOP" stop end if !-- 以下で積分域を決定 !-- 実際には, 積分境界に最近する点 if(x(1)>=bot)then if(x(1)>bot)then write(*,*) "#### WARNING ####" write(*,*) "there is NOT the bot in the x(i)." ! このときは, i_bot=1としておいても, x(1)=bot なので, 余分計算はゼロとなり影響はない. end if i_bot=1 end if if(x(nx)<=top)then if(x(nx)<top)then write(*,*) "#### WARNING ####" write(*,*) "there is NOT the top in the x(i)." ! このとき, i_top=nx としても, x(nx)=top なので, 余分計算はゼロとなる. end if i_top=nx end if do i=2,nx-1 if(x(i)>=bot)then ! i_bot は bot - top 内の最も bot に近づく格子 i_bot=i exit end if end do do i=nx-1,2,-1 if(x(i)<=top)then ! i_top は bot - top 内の最も top に近づく格子 i_top=i exit end if end do !-- 以下で格子に当てはまらない部分を内挿で補完 if(present(undeff))then if(i_bot/=1)then if(y(i_bot)/=undeff.and.y(i_bot-1)/=undeff.and.x(i_bot)/=x(i_bot-1))then y_bot=y(i_bot-1) +((y(i_bot)-y(i_bot-1))/(x(i_bot)-x(i_bot-1)))*(bot-x(i_bot-1)) else y_bot=-y(i_bot) ! 最後の積分でゼロにするための措置 end if else y_bot=-y(i_bot) ! 最後の積分でゼロにするための措置 end if if(i_top/=nx)then if(y(i_top)/=undeff.and.y(i_top+1)/=undeff.and.x(i_top+1)/=x(i_top))then y_top=y(i_top) +((y(i_top+1)-y(i_top))/(x(i_top+1)-x(i_top)))*(x(i_top+1)-top) else y_top=-y(i_top) end if else y_top=-y(i_top) ! 最後の積分でゼロにするための措置 end if else if(i_bot/=1)then y_bot=y(i_bot-1) +((y(i_bot)-y(i_bot-1))/(x(i_bot)-x(i_bot-1)))*(bot-x(i_bot-1)) else y_bot=-y(i_bot) ! 最後の積分でゼロにするための措置 end if if(i_top/=nx)then y_top=y(i_top) +((y(i_top+1)-y(i_top))/(x(i_top+1)-x(i_top)))*(x(i_top+1)-top) else y_top=-y(i_top) ! 最後の積分でゼロにするための措置 end if end if if(i_bot<i_top)then ! 積分区間内に格子が 2 つ以上あるとき if(present(undeff))then do i=i_bot,i_top-1 if(y(i+1)/=undeff.and.y(i)/=undeff)then if(i==i_bot)then res=res+0.5*(x(i)-bot)*(y(i)+y_bot) +0.5*(x(i+1)-x(i))*(y(i+1)+y(i)) ! 下端の余りと通常の短冊分 else res=res+0.5*(x(i+1)-x(i))*(y(i+1)+y(i)) ! 通常の短冊 end if end if end do res=res+0.5*(top-x(i_top))*(y_top+y(i_top)) ! 上端の余りのみ else do i=i_bot,i_top-1 if(i==i_bot)then res=res+0.5*(x(i)-bot)*(y(i)+y_bot) +0.5*(x(i+1)-x(i))*(y(i+1)+y(i)) ! 下端の余りと通常の短冊分 else res=res+0.5*(x(i+1)-x(i))*(y(i+1)+y(i)) ! 通常の短冊 end if end do res=res+0.5*(top-x(i_top))*(y_top+y(i_top)) ! 上端の余りのみ end if else if(present(undeff))then if(y(i_bot)/=undeff)then res=res+0.5*(x(i)-bot)*(y(i)+y_bot) +0.5*(top-x(i))*(y_top+y(i)) else res=0.5*(top-bot)*(y_top+y_bot) end if else res=res+0.5*(x(i_bot)-bot)*(y(i_bot)+y_bot) +0.5*(top-x(i_bot))*(y_top+y(i_bot)) end if end if end subroutine rectangle_int
Subroutine : | |||
x(:,:) : | real, intent(in)
| ||
y(size(x,1),size(x,2)) : | real, intent(in)
| ||
u(size(x,1),size(x,2)) : | real, intent(in)
| ||
v(size(x,1),size(x,2)) : | real, intent(in)
| ||
vecz(size(x,1),size(x,2)) : | real, intent(inout)
| ||
undeff : | real, intent(in), optional |
2ベクトルの外積計算ルーチン vec_prod の 2 次元版 外積なので, 2 次元で計算すると, 外積の 1 成分しか現れないので, intent(inout) の引数が 1 つしかないことに注意.
subroutine vec_prod_2d(x,y,u,v,vecz,undeff) ! 2ベクトルの外積計算ルーチン ! vec_prod の 2 次元版 ! 外積なので, 2 次元で計算すると, 外積の 1 成分しか現れないので, ! intent(inout) の引数が 1 つしかないことに注意. implicit none real, intent(in) :: x(:,:) ! x 方向のベクトル成分 real, intent(in) :: y(size(x,1),size(x,2)) ! y 方向のベクトル成分 real, intent(in) :: u(size(x,1),size(x,2)) ! x 方向のベクトル成分 real, intent(in) :: v(size(x,1),size(x,2)) ! y 方向のベクトル成分 real, intent(inout) :: vecz(size(x,1),size(x,2)) ! 外積の直交成分 real, intent(in), optional :: undeff integer :: i, j, nx, ny nx=size(x,1) ny=size(x,2) if(present(undeff))then do j=1,ny do i=1,nx if(x(i,j)==undeff.or.u(i,j)==undeff.or. y(i,j)==undeff.or.v(i,j)==undeff)then vecz(i,j)=undeff else vecz(i,j)=x(i,j)*v(i,j)-y(i,j)*u(i,j) end if end do end do else do j=1,ny do i=1,nx vecz(i,j)=x(i,j)*v(i,j)-y(i,j)*u(i,j) end do end do end if end subroutine vec_prod_2d
Subroutine : | |||
x(:,:,:) : | real, intent(in)
| ||
y(size(x,1),size(x,2),size(x,3)) : | real, intent(in)
| ||
z(size(x,1),size(x,2),size(x,3)) : | real, intent(in)
| ||
u(size(x,1),size(x,2),size(x,3)) : | real, intent(in)
| ||
v(size(x,1),size(x,2),size(x,3)) : | real, intent(in)
| ||
w(size(x,1),size(x,2),size(x,3)) : | real, intent(in)
| ||
vecx(size(x,1),size(x,2),size(x,3)) : | real, intent(inout)
| ||
vecy(size(x,1),size(x,2),size(x,3)) : | real, intent(inout)
| ||
vecz(size(x,1),size(x,2),size(x,3)) : | real, intent(inout)
| ||
undeff : | real, intent(in), optional |
2ベクトルの外積計算ルーチン
subroutine vec_prod_3d(x,y,z,u,v,w,vecx,vecy,vecz,undeff) ! 2ベクトルの外積計算ルーチン implicit none real, intent(in) :: x(:,:,:) ! x 方向のベクトル成分 real, intent(in) :: y(size(x,1),size(x,2),size(x,3)) ! y 方向のベクトル成分 real, intent(in) :: z(size(x,1),size(x,2),size(x,3)) ! z 方向のベクトル成分 real, intent(in) :: u(size(x,1),size(x,2),size(x,3)) ! x 方向のベクトル成分 real, intent(in) :: v(size(x,1),size(x,2),size(x,3)) ! y 方向のベクトル成分 real, intent(in) :: w(size(x,1),size(x,2),size(x,3)) ! z 方向のベクトル成分 real, intent(inout) :: vecx(size(x,1),size(x,2),size(x,3)) ! 外積の x 成分 real, intent(inout) :: vecy(size(x,1),size(x,2),size(x,3)) ! 外積の y 成分 real, intent(inout) :: vecz(size(x,1),size(x,2),size(x,3)) ! 外積の z 成分 real, intent(in), optional :: undeff integer :: i, j, k, nx, ny, nz nx=size(x,1) ny=size(x,2) nz=size(x,3) if(present(undeff))then !$omp parallel default(shared) !$omp do schedule(dynamic) private(i,j,k) do k=1,nz do j=1,ny do i=1,nx if(x(i,j,k)==undeff.or.u(i,j,k)==undeff.or.y(i,j,k)==undeff.or. v(i,j,k)==undeff.or.z(i,j,k)==undeff.or.w(i,j,k)==undeff)then vecx(i,j,k)=undeff vecy(i,j,k)=undeff vecz(i,j,k)=undeff else vecx(i,j,k)=y(i,j,k)*w(i,j,k)-z(i,j,k)*v(i,j,k) vecy(i,j,k)=z(i,j,k)*u(i,j,k)-x(i,j,k)*w(i,j,k) vecz(i,j,k)=x(i,j,k)*v(i,j,k)-y(i,j,k)*u(i,j,k) end if end do end do end do !$omp end do !$omp end parallel else !$omp parallel default(shared) !$omp do schedule(runtime) private(i,j,k) do k=1,nz do j=1,ny do i=1,nx vecx(i,j,k)=y(i,j,k)*w(i,j,k)-z(i,j,k)*v(i,j,k) vecy(i,j,k)=z(i,j,k)*u(i,j,k)-x(i,j,k)*w(i,j,k) vecz(i,j,k)=x(i,j,k)*v(i,j,k)-y(i,j,k)*u(i,j,k) end do end do end do !$omp end do !$omp end parallel end if end subroutine vec_prod_3d