fortran

フォーラム(掲示板)ルール
フォーラム(掲示板)ルールはこちら  ※コードを貼り付ける場合は [code][/code] で囲って下さい。詳しくはこちら
太郎

fortran

#1

投稿記事 by 太郎 » 5年前

関数cranknicholson(fai,t,dt,x,dx,xmax)の中身でどこが間違っているか探すためにとりあえず、その関数の中身をすべてコメント文にして、上から元に戻して間違いを探すと、a()を入力するとオーバーフローします。なぜですか。お願します。{{{このプログラムは陽解法、陰解法、クランク・ニコルソン法で解を求めるプログラムです。クランク・ニコルソン法の関数を追加したらオーバーフローしました。}}}

太郎

Re: fortran

#2

投稿記事 by 太郎 » 5年前

コードを貼るの忘れてました

コード:

 
 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
 

アバター
あたっしゅ
記事: 663
登録日時: 13年前
住所: 東京23区
連絡を取る:

Re: fortran

#3

投稿記事 by あたっしゅ » 5年前

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, 電子ブロック 持ち。

太郎

Re: fortran

#4

投稿記事 by 太郎 » 5年前

説明不足すぎました。
サブルーチンのcranknicholson(fai,t,dt,x,dx,xmax)関数内で、例えばa(1,1)=rを書き足してコンパイルすると、オーバーフローします。つまり2次元配列を書き足すとオーバーフローします。わかる方お願いします。

かずま

Re: fortran

#5

投稿記事 by かずま » 5年前

subroutine cranknicholson の中の
double precision a(10000,10000),a2(10000,10000) は、
16億バイトの領域を必要とします。

32ビット環境だと、アドレス空間は 4ギガバイトで、約43億。
ユーザプログラムに与えられるメモリはその半分以下かもしれません。
プログラム中の 10000 をすべて 5000 にしてみたらどうなりますか?

OS やコンパイラが 64ビットに対応しているのでしょうか?
フォーラム(掲示板)ルールにあるように、
環境やエラーメッセージなどの正確な情報をお願いします。

かずま

Re: fortran

#6

投稿記事 by かずま » 5年前

かずま さんが書きました:
5年前
subroutine cranknicholson の中の
double precision a(10000,10000),a2(10000,10000) は、
16億バイトの領域を必要とします。
inkaihou の中の a(10000,10000) や、gauss の中の dum(10000,10000) も
それぞれ 8億バイトです。

ax の中の a や、gauss の中の a、 pivotselection の中の a は、
引数なので呼び出し元の領域が参照されるだけで大きな領域は
必要ないでしょう。

太郎

Re: fortran

#7

投稿記事 by 太郎 » 5年前

その通りでした。
ありがとうございます。

かずま

Re: fortran

#8

投稿記事 by かずま » 5年前

どのようにして解決したのですか?
環境やエラーメッセージなどの正確な情報は提示できないのですか?
OS やコンパイラは 64ビットに対応していないのですか?
フォーラム掲示板ルールは読みましたか?
プログラム中の 10000 をすべて 5000 にしてみたらどうなりましたか?
なぜ、回答者の質問に答えないのですか?

アバター
あたっしゅ
記事: 663
登録日時: 13年前
住所: 東京23区
連絡を取る:

Re: fortran

#9

投稿記事 by あたっしゅ » 5年前

今回の FORTRAN の件で感じた違和感。
よく考えたら、自分の中で、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, 電子ブロック 持ち。

返信

“C言語何でも質問掲示板” へ戻る