fortran
Re: fortran
コードを貼るの忘れてました
program main
implicit none
double precision fai(10000),t,dt,x,dx,xmax
dx=0.1d0
xmax=2.0d0
dt=0.001d0
t=0.05d0
call youkaihou(fai,t,dt,x,dx,xmax)
call inkaihou(fai,t,dt,x,dx,xmax)
!call cranknicholson(fai,t,dt,x,dx,xmax)
end program main
subroutine youkaihou(fai,t,dt,x,dx,xmax)!t秒後の温度を求める。
implicit none
double precision fai0(10000),fai(10000),dum(10000),t,dt,x,dx,xmax,pi,r
integer n,m,i,j
pi=3.14159265359d0
x=0.0d0
n=dble(xmax/dx)
do i=1,n
x=x+dx
fai0(i)=sin(pi*x/2.0d0)
enddo
m=dble(t/dt)
r=dt/(dx**2.0d0)
do i=1,m
do j=1,n
if(i==1)then
if(j==1) then
dum(j)=(1.0d0-2.0d0*r)*fai0(j)+r*fai0(j+1)
else if(j==n) then
dum(j)=r*fai0(n-1)+(1.0d0-2.0d0*r)*fai0(n)
else!1<j<n-1のとき
dum(j)=r*fai0(j-1)+(1.0d0-2.0d0*r)*fai0(j)+r*fai0(j)
endif
else !i>2のとき
if(j==1) then
dum(j)=(1.0d0-2.0d0*r)*fai(j)+r*fai(j+1)
else if(j==n) then
dum(j)=r*fai(n-1)+(1.0d0-2.0d0*r)*fai(n)
else!1<j<n-1のとき
dum(j)=r*fai(j-1)+(1.0d0-2.0d0*r)*fai(j)+r*fai(j)
endif
endif
enddo
do j=1,n
fai(j)=dum(j)
enddo
enddo
do i=1,n
write(*,*)fai(i)
enddo
write(*,*)' '
end subroutine
subroutine inkaihou(fai,t,dt,x,dx,xmax)
implicit none
double precision fai0(10000),fai(10000),fai2(10000),dum(10000),t,dt,x,dx,xmax,pi,r,a(10000,10000)
integer n,m,i,j
pi=3.14159265359d0
x=0.0d0
n=dble(xmax/dx)
do i=1,n
x=x+dx
fai0(i)=sin(pi*x/2.0d0)
enddo
r=dt/(dx**2.0d0)
do i=1,n
a(i,i)=1+2*r
if(i<n)then
a(i,i+1)=-r
endif
if(i>1)then
a(i,i-1)=-r
endif
enddo
m=dble(t/dt)
do i=1,m
if(i==1) then
call gauss(a,fai,fai0,n)
else
call gauss(a,fai2,fai,n)
fai=fai2
endif
enddo
do i=1,n
write(*,*)fai(i)
enddo
write(*,*)' '
end subroutine
subroutine cranknicholson(fai,t,dt,x,dx,xmax)
implicit none
double precision fai0(10000),fai(10000),fai2(10000),dum(10000),t,dt,x
double precision dx,xmax,pi,r,a(10000,10000),a2(10000,10000),y(10000)
integer n,m,i,j
!pi=3.14159265359d0
!x=0.0d0
!n=dble(xmax/dx)
!do i=1,n
! x=x+dx
! fai0(i)=sin(pi*x/2)
!enddo
r=dt/(dx**2.0d0)
!a(1,1)=1+r
!write(*,*)r,n
!do i=1,n
! a(i,i)=1+r
! if(i<n)then
! a(i,i+1)=-r/2.0d0
! endif
! if(i>1)then
! a(i,i-1)=-r/2.0d0
! endif
!enddo
!do i=1,n
! a2(i,i)=1-r
! if(i<n)then
! a2(i,i+1)=r/2.0d0
! endif
! if(i>1)then
! a2(i,i-1)=r/2.0d0
! endif
!enddo
m=dble(t/dt)
!do i=1,m
!if(i==1) then
!call ax(a2,fai0,n,y) !{n*nの行列A}*{1*nの行列x}
!call gauss(a,fai,y,n)
!else
!call ax(a2,fai,n,y)
!call gauss(a,fai2,y,n)
!fai=fai2
!endif
!enddo
!do i=1,n
!write(*,*)fai(i)
!enddo
write(*,*)' '
end subroutine
subroutine ax(a,x,n,ans)!{n*nの行列A}*{1*nの行列x}
implicit none
double precision a(10000,10000),x(10000),dum(10000),ans(10000)
integer i,j,n
dum=0
do i=1,n
do j=1,n
dum(i)=dum(i)+x(j)*a(i,j)
enddo
enddo
do i=1,n
ans(i)=dum(i)
enddo
end subroutine
subroutine gauss(a,x,y,n)
implicit none
double precision a(10000,10000),x(10000),y(10000),kari,kari2,dum(10000,10000)
integer n,i,j,k,h,p
do i=1,n
do j=1,n
dum(i,j)=a(i,j)
enddo
enddo
h=0
!前進消去43~70
do i=1,n
do j=i,n
if(dum(j,i)/=0) then
exit
endif
if(j==n) then
write(*,*)"Unique solution does not exist"
h=1
endif
enddo
enddo
if(h==0)then
do i=1,n
call pivotselection(i,dum,y,n)
do p=1,n
enddo
do j=i,n
y(i)=y(i)/dum(i,i)
kari=dum(i,i)
if(j<n)then
kari2=dum(j+1,i)
endif
do k=i,n
dum(j,k)=dum(j,k)/kari
if(j<n)then
dum(j+1,k)=dum(j+1,k)-kari2*dum(i,k)
if(k==i)then
y(j+1)=y(j+1)-kari2*y(i)
endif
endif
enddo
enddo
enddo
endif
!後退代入
x(n)=y(n)/dum(n,n)
do i=1,n-1
x(n-i)=y(n-i)/dum(n-i,n-i)
do j=1,i
x(n-i)=x(n-i)-dum(n-i,n-j+1)*x(n-j+1)/dum(n-i,n-i)
enddo
enddo
end subroutine
subroutine pivotselection(s,a,y,n)
implicit none
double precision a(10000,10000),y(10000),saidai,kari
integer i,j,s,n!jは最大の係数をもつのは何行目かを保存、sは対角線上の成分の添え字、nはn*n行列のn
saidai=abs(a(s,s))
j=s
do i=s,n-1
if(saidai<abs(a(i+1,s))) then
saidai=a(i+1,s)
j=i+1
endif
enddo
do i=s,n
kari=a(s,i)
a(s,i)=a(j,i)
a(j,i)=kari
enddo
kari=y(s)
y(s)=y(j)
y(j)=kari
end subroutine
Re: fortran
FORTRAN だと !
とりあえず、「gnu fortran」とか検索してみたら gfortran 出てきたので、
bash on win10(今の正式名称は、なんでしたっけ ?) にダウンロード & インストール。
mixpp,f90 と名付けて、コンパイル。実行。で、
2.5884155414774054E-002
4.1881457553708118E-002
6.7706545758653772E-002
0.10835679610097523
0.16842825478795440
0.24945327123284969
0.34831009271007529
0.45794307685445584
0.56964442703431817
0.67526932461040134
0.76840252845341972
0.84451120711731886
0.90059751057281789
0.93478948229364023
0.94605606140369880
0.93405493279399832
0.89906162635943188
0.84193224211845663
0.76407209596640713
0.36008358621827546
0.13833412222338087
0.27326200759827185
0.40146130647950079
0.51977539973673348
0.62529123225480154
0.71541135206969553
0.78791872207852620
0.84103353015785864
0.87346242081723113
0.88444257997300213
0.87378595175106533
0.84193158098035081
0.79001341516856338
0.71994039887216388
0.63445766576933715
0.53711208017745915
0.43200590533877664
0.32324514300544255
0.21413376452722446
0.10641095232193158
一応、動いてるようですが。
こうなると、動かない原因は、fortran の規格違い、とかじゃないんですか ?
いや、fortran 全然知らないんで、的外れかもしれませんが。
今回、初めて動かしたもので...
オーバーフローだと、浮動小数点のサイズの問題 ?
とりあえず、「gnu fortran」とか検索してみたら gfortran 出てきたので、
bash on win10(今の正式名称は、なんでしたっけ ?) にダウンロード & インストール。
mixpp,f90 と名付けて、コンパイル。実行。で、
2.5884155414774054E-002
4.1881457553708118E-002
6.7706545758653772E-002
0.10835679610097523
0.16842825478795440
0.24945327123284969
0.34831009271007529
0.45794307685445584
0.56964442703431817
0.67526932461040134
0.76840252845341972
0.84451120711731886
0.90059751057281789
0.93478948229364023
0.94605606140369880
0.93405493279399832
0.89906162635943188
0.84193224211845663
0.76407209596640713
0.36008358621827546
0.13833412222338087
0.27326200759827185
0.40146130647950079
0.51977539973673348
0.62529123225480154
0.71541135206969553
0.78791872207852620
0.84103353015785864
0.87346242081723113
0.88444257997300213
0.87378595175106533
0.84193158098035081
0.79001341516856338
0.71994039887216388
0.63445766576933715
0.53711208017745915
0.43200590533877664
0.32324514300544255
0.21413376452722446
0.10641095232193158
一応、動いてるようですが。
こうなると、動かない原因は、fortran の規格違い、とかじゃないんですか ?
いや、fortran 全然知らないんで、的外れかもしれませんが。
今回、初めて動かしたもので...
オーバーフローだと、浮動小数点のサイズの問題 ?
VTuber:
東上☆海美☆(とうじょう・うみみ)
http://atassyu.php.xdomain.jp/vtuber/index.html
レスがついていないものを優先して、レスするみみ。時々、見当外れなレスしみみ。
中の人:
手提鞄あたッしュ、[MrAtassyu] 手提鞄屋魚有店
http://ameblo.jp/mratassyu/
Pixiv: 666303
Windows, Mac, Linux, Haiku, Raspbery Pi, Jetson Nano, 電子ブロック 持ち。
東上☆海美☆(とうじょう・うみみ)
http://atassyu.php.xdomain.jp/vtuber/index.html
レスがついていないものを優先して、レスするみみ。時々、見当外れなレスしみみ。
中の人:
手提鞄あたッしュ、[MrAtassyu] 手提鞄屋魚有店
http://ameblo.jp/mratassyu/
Pixiv: 666303
Windows, Mac, Linux, Haiku, Raspbery Pi, Jetson Nano, 電子ブロック 持ち。
Re: fortran
説明不足すぎました。
サブルーチンのcranknicholson(fai,t,dt,x,dx,xmax)関数内で、例えばa(1,1)=rを書き足してコンパイルすると、オーバーフローします。つまり2次元配列を書き足すとオーバーフローします。わかる方お願いします。
サブルーチンのcranknicholson(fai,t,dt,x,dx,xmax)関数内で、例えばa(1,1)=rを書き足してコンパイルすると、オーバーフローします。つまり2次元配列を書き足すとオーバーフローします。わかる方お願いします。
Re: fortran
subroutine cranknicholson の中の
double precision a(10000,10000),a2(10000,10000) は、
16億バイトの領域を必要とします。
32ビット環境だと、アドレス空間は 4ギガバイトで、約43億。
ユーザプログラムに与えられるメモリはその半分以下かもしれません。
プログラム中の 10000 をすべて 5000 にしてみたらどうなりますか?
OS やコンパイラが 64ビットに対応しているのでしょうか?
フォーラム(掲示板)ルールにあるように、
環境やエラーメッセージなどの正確な情報をお願いします。
double precision a(10000,10000),a2(10000,10000) は、
16億バイトの領域を必要とします。
32ビット環境だと、アドレス空間は 4ギガバイトで、約43億。
ユーザプログラムに与えられるメモリはその半分以下かもしれません。
プログラム中の 10000 をすべて 5000 にしてみたらどうなりますか?
OS やコンパイラが 64ビットに対応しているのでしょうか?
フォーラム(掲示板)ルールにあるように、
環境やエラーメッセージなどの正確な情報をお願いします。
Re: fortran
inkaihou の中の a(10000,10000) や、gauss の中の dum(10000,10000) もかずま さんが書きました: ↑5年前subroutine cranknicholson の中の
double precision a(10000,10000),a2(10000,10000) は、
16億バイトの領域を必要とします。
それぞれ 8億バイトです。
ax の中の a や、gauss の中の a、 pivotselection の中の a は、
引数なので呼び出し元の領域が参照されるだけで大きな領域は
必要ないでしょう。
Re: fortran
どのようにして解決したのですか?
環境やエラーメッセージなどの正確な情報は提示できないのですか?
OS やコンパイラは 64ビットに対応していないのですか?
フォーラム掲示板ルールは読みましたか?
プログラム中の 10000 をすべて 5000 にしてみたらどうなりましたか?
なぜ、回答者の質問に答えないのですか?
環境やエラーメッセージなどの正確な情報は提示できないのですか?
OS やコンパイラは 64ビットに対応していないのですか?
フォーラム掲示板ルールは読みましたか?
プログラム中の 10000 をすべて 5000 にしてみたらどうなりましたか?
なぜ、回答者の質問に答えないのですか?
Re: fortran
今回の FORTRAN の件で感じた違和感。
よく考えたら、自分の中で、FORTRAN のイメージと COBOL のイメージが、ごっちゃになっていた。
というわけで、Unix/Linux 上の GnuCOBOL で、メガデモを作ることにしました。
一度、C のソースコードになるので、C/C++ とリンク可能。
十進演算ランタイムが圧縮できねえよ。
よく考えたら、自分の中で、FORTRAN のイメージと COBOL のイメージが、ごっちゃになっていた。
というわけで、Unix/Linux 上の GnuCOBOL で、メガデモを作ることにしました。
一度、C のソースコードになるので、C/C++ とリンク可能。
十進演算ランタイムが圧縮できねえよ。
VTuber:
東上☆海美☆(とうじょう・うみみ)
http://atassyu.php.xdomain.jp/vtuber/index.html
レスがついていないものを優先して、レスするみみ。時々、見当外れなレスしみみ。
中の人:
手提鞄あたッしュ、[MrAtassyu] 手提鞄屋魚有店
http://ameblo.jp/mratassyu/
Pixiv: 666303
Windows, Mac, Linux, Haiku, Raspbery Pi, Jetson Nano, 電子ブロック 持ち。
東上☆海美☆(とうじょう・うみみ)
http://atassyu.php.xdomain.jp/vtuber/index.html
レスがついていないものを優先して、レスするみみ。時々、見当外れなレスしみみ。
中の人:
手提鞄あたッしュ、[MrAtassyu] 手提鞄屋魚有店
http://ameblo.jp/mratassyu/
Pixiv: 666303
Windows, Mac, Linux, Haiku, Raspbery Pi, Jetson Nano, 電子ブロック 持ち。