くだすれFORTRAN（超初心者用）その4

1 名前：デフォルトの名無しさん [2009/01/24(土) 18:32:01 ] このスレッドは、他のスレッドでは書き込めない超低レベル、 もしくは質問者自身何が何だか分からない質問を勇気を持って書き込むスレッドです。 FORTRAN使いが優しくコメントを返しますが、 お礼はFORTRANの布教と初心者の救済と次期Fortran2008規格でのCOMEFROM文採用をお願いします。 ●注意事項 ・質問する前にGoogle等の検索サイトで検索しましょう。 ・回答者にわかりやすい様に、質問内容はできる限り詳しく書きましょう。 ・エラーの場合は起きた状況、環境（OS・コンパイラ・バージョン）、エラーメッセージも詳しく書きましょう。 ●前スレ くだすれFORTRAN（超初心者用）その3 pc11.2ch.net/test/read.cgi/tech/1196384126/ くだすれFORTRAN（超初心者用）その2 pc11.2ch.net/test/read.cgi/tech/1164121236/ くだすれFORTRAN（超初心者用） pc8.2ch.net/test/read.cgi/tech/1138063703/ ●関連スレ FORTRAN �W pc11.2ch.net/test/read.cgi/tech/1163319215/ 292 名前：デフォルトの名無しさん mailto:sage [2009/06/16(火) 00:33:41 ] >>291 そう、それをもう一度頼む！(><) MSのFortran部隊はどこ行っちまったんだ？ 293 名前：デフォルトの名無しさん mailto:sage [2009/06/16(火) 00:49:49 ] >>292 Netscapeがブラウザー戦争を仕掛けたせいで、VC++部隊以外が全部IE開発に 回されたという噂が当時あった。 まぁMSの歴史をたどると、BASICの次に作って売ったのはFORTRANだから MS Visual FORTRANは原点回帰といえなくもないが、VBすら捨ててしまうMS様は 容赦ないのだｗ でも、高い金出してコンパイラ律儀に買ってくれるのはFORTRANユーザーくらい？ 294 名前：デフォルトの名無しさん mailto:sage [2009/06/16(火) 00:55:45 ] >>292 MS → DEC → Compaq → Intel 295 名前：デフォルトの名無しさん mailto:sage [2009/06/16(火) 13:53:00 ] >>294 MSからDECに流れたとはしらなんだ IntelとMSがもっとみっちり連携してくれればいいんだけど、またWintelって非難が起こりそうだな 296 名前：arrayの地獄 mailto:sage [2009/06/16(火) 22:24:17 ] 時間がかかって仕方ないのだが何とかなりませんか？ ifort 9.1です。 real(8)::x(25,25),y(50,50) x=なんたら y=0 y(1:25,1:25)=x(1:25,1:25) y=cshift(y,12,1) y=cshift(y,12,2) 特にy=xのところ。 なんでこんなんに時間食う？？ （メモリ書き込みに難儀してるのはわかるが、 y全体がキャッシュ上に乗るから、キャッシュで整えてから書けば時間かからんはず・・） 297 名前：デフォルトの名無しさん mailto:sage [2009/06/16(火) 22:27:38 ] >>296 それだけじゃ分からんっつの 時間食うってどれくらい？ 298 名前：arrayの地獄 mailto:sage [2009/06/16(火) 22:46:24 ] ９２ｘ４９回で数秒のオーダー ＠xeon5160？ 1core使用 なんかおかしいな・・もうちょっと調べます。 299 名前：デフォルトの名無しさん mailto:sage [2009/06/16(火) 22:49:48 ] ド初心者ですが質問です。コンパイルできたファイルを実行すると fmt: end of file apparent state: unit 10 named inp-2.dat last format: (F9.5) lately reading sequential formatted external IO となるのですが何が問題なんでしょうか？プログラムは外部データ（inp-2.dat）から数値を読み最大値を求めるものです 数値は実数です 300 名前：デフォルトの名無しさん mailto:sage [2009/06/16(火) 22:56:05 ] >>296 >y(1:25,1:25)=x(1:25,1:25) なぜ素直にｙ＝ｘとしない。　 多分 ALLOCATE（TMP(26,26)) TMP=x ｙ＝TMP DEALLOCATE（TMP) のように展開されるから遅いんだろ。 301 名前：デフォルトの名無しさん mailto:sage [2009/06/16(火) 22:59:20 ] >>300 TMP(25,25)だったｗ　0:25と勘違いしたよ。1:25とか俺的にはあり得ないｗ >>299 データ数が足りないように見受けられる。 F9.5に当たるところでデータが足りないか、変数の数が合わないかしてるのでは？ 302 名前：296です。 mailto:sage [2009/06/16(火) 23:06:23 ] [hoge@xeon prog20]$ cat test.f90 program test implicit none integer,parameter::nx=25,ny=50 real(8)::x(nx,nx),y(ny,ny) integer::i x=reshape((/1:25*25/),shape(x)) do i=1,100000 call sub1(x,y,nx,ny) call sub2(x,y,nx,ny) end do end program test subroutine sub1(x,y,nx,ny) implicit none real(8)::x(nx,nx),y(ny,ny) integer::nx,ny y=0 y(1:nx,1:nx)=x(1:nx,1:nx) y=cshift(y,12,1) y=cshift(y,12,2) end subroutine subroutine sub2(x,y,nx,ny) implicit none real(8)::x(nx,nx),y(nx,nx) !nyでなくnx integer::nx,ny y=0 y(1:nx,1:nx)=x(1:nx,1:nx) y=cshift(y,12,1) y=cshift(y,12,2) end subroutine　続く 303 名前：デフォルトの名無しさん mailto:sage [2009/06/16(火) 23:07:10 ] >>301 データ数ですか。それは入力元のファイルの値に問題ありということですか？プログラムは IMPLICIT REAL(A-H,O-Z) DIMENSION X(100) OPEN(UNIT=10,FILE='inp-2.dat',STATUS='old',FORM='FORMATTED') OPEN(UNIT=11,FILE='list-2.dat',STATUS='UNKNOWN',FORM='FORMATTED') READ(10,100)N DO 20 I=1,N READ(10,200)X(I) 20 CONTINUE IMAX=0.0 DO 30 I=1,N IF(X(I).GT.IMAX)MAXN=I IF(X(I).GT.IMAX)IMAX=X(I) 30 CONTINUE WRITE(11,300)'MAX=',IMAX,'NO.=',MAXN 100 FORMAT(A5) 200 FORMAT(F9.5) 300 FORMAT(A4,F10.5,A4,I3) CLOSE(10) CLOSE(11) STOP END となってます。サンプルプログラムを改変しただけなのでおかしいところが多いと思います。入力ファイルはとりあえず実数を5つ縦に並べてあります 304 名前：296です。 mailto:sage [2009/06/16(火) 23:07:11 ] [hoge@xeon prog20]$ ifort test.f90 -p [hoge@xeon prog20]$ time a.out real 0m2.359s user 0m2.357s sys 0m0.002s [hoge@xeon prog20]$ gprof -b Each sample counts as 0.01 seconds. % cumulative self self total time seconds seconds calls s/call s/call name 71.83 1.53 1.53 100000 0.00 0.00 sub1_ 28.17 2.13 0.60 100000 0.00 0.00 sub2_ 0.00 2.13 0.00 1 0.00 2.13 MAIN__ 305 名前：296です。 mailto:sage [2009/06/16(火) 23:13:31 ] みにくてすみませんでした。 CPUはcore2の6600 2.4Gでcashe＝4Mでした。 yとxの寸法が違うので、y=xにできないはずです。 逆フーリエ変換にかけるとき外を0にしておきたいので、 こんなになっています。 実プログラムでは逆フーリエ変換にかかる時間と 同じくらいの時間がかかって大損しています。 306 名前：296です。 mailto:sage [2009/06/16(火) 23:18:31 ] まとめると、 y(1:25,1:25)=x(1:25,1:25) を100000回実行すると、 real(8)::x(25,25),y(50,50) のとき　=>1.53秒かかる。 real(8)::x(25,25),y(25,25) のとき　=>0.60秒かかる。 2.5倍もかからんでええやろ ということです。 307 名前：296です。 mailto:sage [2009/06/16(火) 23:27:05 ] >>299 私もできる範囲で手伝い。 100 FORMAT(A5) が気になります。 実際に実行させたらNにどんな値が入るでしょうか？ 思っている値かどうかwrite文で確かめるのが吉かと思います。 308 名前：デフォルトの名無しさん mailto:sage [2009/06/16(火) 23:32:44 ] >>307 このままだと実行すらできないんですよね・・・　 READ（〜）N　　のあとに N=5　　としてみましたが実行エラーの文は変わりませんでした 309 名前：デフォルトの名無しさん mailto:sage [2009/06/16(火) 23:40:27 ] >>296 corei7のifort 11ですが $ cat test.f90 integer,parameter::n=25 real(8)::x(n,n),y(2*n,2*n) x=0 y=1 !y=x y(1:n,1:n)=x(1:n,1:n) end $ time ./a.out real 0m0.001s user 0m0.000s sys 0m0.001s n=25000だと real 0m1.974s user 0m0.754s sys 0m1.217s 何か環境おかしいんじゃないんでしょうか？他のホストで試したどうっすかね 310 名前：296です。 mailto:sage [2009/06/16(火) 23:44:15 ] >>299 適当にinp-2.datに 5 1 12 123 1234 12345 123456 と入れてます。 read(〜)Nのあとにwrite(*,*)Nとすれば、ここまでは走りませんか？ これで走らせると、私の場合はN= 538976288と、とんでもないことになってました。 数字を数として読み込みたいのに、文字として読み込んだらまずいですよね、やっぱり。 A5はI5ぐらいかと思いますが、どうでしょう？ 311 名前：296です。 mailto:sage [2009/06/16(火) 23:49:47 ] >>309 情報ありがとうです。 corei7＋ifort 11ですか・・ なるほど。別のマシンで試してみます。 312 名前：デフォルトの名無しさん mailto:sage [2009/06/16(火) 23:54:28 ] >>310 inp-2.datに小数点以下１桁を含んだ実数を適当に5つ入れて、writeしてみたところN＝1073741824となってました。 フォーマットをi5にすると実数を受け取れないのでF5.5にしてみたのですが、Nは上記の通りです。何が起こってるんですかね・・・？ 313 名前：296です。 mailto:sage [2009/06/16(火) 23:55:32 ] >>309 すんません。メインルーチンで do i=1,100000 と回してますので、timeの結果が直接比較できないです。 すみませんが追加でテストお願いできませんか？ 314 名前：デフォルトの名無しさん mailto:sage [2009/06/17(水) 00:14:43 ] >>313 いやメインルーチンとかじゃなくて >>309でテストすりゃいいじゃん 315 名前：296です。 mailto:sage [2009/06/17(水) 00:18:42 ] >>299 Nはデータの個数ですよね。実数じゃまずいのでは？ inp-2.datは データの個数 データ１ データ２ ： というような構造じゃありませんか？ 少なくとも最初のデータだけは整数かな。 それとIMAXを実数宣言しなきゃいけないようです。 316 名前：デフォルトの名無しさん mailto:sage [2009/06/17(水) 00:27:14 ] >>315 datファイルにデータの個数を入れてなかったんですね・・・やっとわかりました！ありがとうございます。実行できました プログラム自体はうまくいってないんですが・・・ｗ 実行結果がMAX＝０．０００・・・、No.＝２となっているんでうまくいってないですね 317 名前：296です。 mailto:sage [2009/06/17(水) 00:33:29 ] >>314 お世話になってます。 確かにこちらでも結果0秒になりました。 いやしかし、たぶん0.000,0060秒（6マイクロ秒）くらいかかっているのが、 timeコマンドの分解能不足で測れてないだけの様な気がしてますです。 実際は100000ｘ100000回とかもっとぐらい実行するので、何とかしたいです。 >>316 IMAXの方もよろしく。 318 名前：309 mailto:sage [2009/06/17(水) 00:39:12 ] >>317 問題点を整理（できればテストコードを用意）して頂ければ協力できますよ 他のマシンもいくつかあります 319 名前：デフォルトの名無しさん mailto:sage [2009/06/17(水) 00:53:40 ] >>317 実数宣言とは何でしょうか？imax=0.0ではダメなんですよね？ 320 名前：300 mailto:sage [2009/06/17(水) 01:02:12 ] >>305 >yとxの寸法が違うので、y=xにできないはずです。 すまんこｗ　ようつべ見ながらだったので、よく見なかったｗｗ こんなものに時間がかかるも糞もないだろと思ったが、 １０＾５＊１０＾５回も繰り返すなら気になるかｗ 部分配列でメモリーが連続していないので、コンパイラのバージョンが古いと 余計なテンポラリを暗黙にとっている可能性が高いと思う。 DO i = 1, 25 y(1:25, i) = x(:, i) END DO と、でもすれば早くなるべ。最適化でアンローリングとかもあろうし。 321 名前：296です。 mailto:sage [2009/06/17(水) 01:07:49 ] >>318 ありがとうございます！ テストコードを302に書き込みました。これでお願いできたらありがたいです！ 問題点のほうは、まとめると、 302のプログラムで y(1:25,1:25)=x(1:25,1:25) を100000回実行すると、 real(8)::x(25,25),y(50,50) のとき　=>1.53秒かかる。 real(8)::x(25,25),y(25,25) のとき　=>0.60秒かかる。 「2.5倍もかからんでええやろ」 ということです。 おわかりいただけたでしょうか？ ifort10.1を試みてみました。XEON5160です。 % time a10.out　　　　　：ifort10.1 1.288u 0.000s 0:01.31 97.7% % time a09.out　　　　　：ifort9.1 1.916u 0.004s 0:01.96 97.4% だいぶ差があるなあ・・というところです。 322 名前：296です。 mailto:sage [2009/06/17(水) 01:16:14 ] >>319 real imax を宣言文のところに追加です。 FORTRANの慣習ではI〜Nから始まる変数名は整数変数に使うことになっているので、 imaxのままだと動くけどちょっと減点かな？ imax->xmaxにするほうがいいかも。この場合real xmaxは不要。 >>320 ご指摘ありがとうですです。 さっそく試してみますです！ 323 名前：デフォルトの名無しさん mailto:sage [2009/06/17(水) 01:21:31 ] >>322 うまくいきました！　なんとなくfortranがわかりはじめてきました！ ありがとうございます 324 名前：デフォルトの名無しさん mailto:sage [2009/06/17(水) 01:30:04 ] 初心者は呪文のようにimplicit noneを宣言しておきましょう。 325 名前：309 mailto:sage [2009/06/17(水) 01:39:38 ] >>321 こんな感じでした。時間はreal timeです Corei7 920、ifort 11.0 real(8)::x(25,25),y(50,50)・・・0m0.799s real(8)::x(25,25),y(25,25)・・・0m0.255s Xeon5472、ifort 10.1 real(8)::x(25,25),y(50,50)・・・0m0.922s real(8)::x(25,25),y(25,25)・・・0m0.280s Opteron2218、pgf90 6.2 real(8)::x(25,25),y(50,50)・・・0m9.821s real(8)::x(25,25),y(25,25)・・・0m4.372s >>300さんの指摘通りってことですかね Opteronが糞遅いのが気になりますが・・・（；^ω^） 326 名前：296です。 mailto:sage [2009/06/17(水) 01:56:29 ] >>309様（325） ありがとうございます！！お手数掛けました。 COREi7万歳というところですか。 320様のご指摘も含めていろいろ検討したいのですが、 今晩はこの辺で失礼します。 （いろいろありがとうございました。また明日頑張ります。） 327 名前：デフォルトの名無しさん mailto:sage [2009/06/17(水) 03:05:28 ] corei7はかなり評判いいよね opteronを完全に引き離した 328 名前：デフォルトの名無しさん mailto:sage [2009/06/17(水) 03:11:03 ] >>325 OpteronというよりPGIのコンパイラのせいのような気がする。 オプションがどうなってるのかわからんし。 最新版は8.0だし。 でもたしかにi7はなかなかいい。 329 名前：309 mailto:sage [2009/06/17(水) 13:56:34 ] >>328 ＞OpteronというよりPGIのコンパイラのせいのような気がする。 その通りでしたｗ >>325のコンパイルは全てオプション指定なしでして、調べたところ ifortはデフォルトで-O2が、pgf90は-O1がかかるようです。 -O2を付けた場合↓ real(8)::x(25,25),y(50,50)・・・0m1.306s real(8)::x(25,25),y(25,25)・・・0m0.784s 330 名前：デフォルトの名無しさん mailto:sage [2009/06/17(水) 20:36:22 ] >>301 ＞1:25とか俺的にはあり得ないｗ Cみたいに0から開始ってこと？ Fortranだと1から始めるのが一般的だと思うけど 331 名前：デフォルトの名無しさん mailto:sage [2009/06/17(水) 23:16:26 ] >>330 いや、：だけで済ますとか、それすら書かないと言うこと。 332 名前：デフォルトの名無しさん mailto:sage [2009/06/18(木) 18:02:57 ] >>331 そういうことね で、実際どれがいいんだろうか。例えば↓の３つでは �@ x(1:n) = y(1:n) �A x(:) = y(:) �B x = y 自分は�Aが良いと思ってる。�@は冗長だし、�Bは簡潔だけど配列であることが自明じゃないので。 問題はどれが一番速いかなんだけど・・・ ↓のサイトでは「並列最適化コンパイラの動作を助け、バグを少なくするため」�@が推奨されてるし ttp://www.mri-jma.go.jp/Project/mrinpd/coderule.html でも、これって>>300氏の言ってることとは反対だよね。 333 名前：デフォルトの名無しさん [2009/06/18(木) 19:19:25 ] 例えば、配列 data1 から data4 があって、次々にsubroutine に引き渡したい場合ですが、 call hoge(data1) call hoge(data2) call hoge(data3) call hoge(data4) しないで、 do i=1,4 call hoge( ) end do の様にまとめる事って加能ですか？　もし可能ならやり方教えてください。 334 名前：デフォルトの名無しさん [2009/06/18(木) 19:22:10 ] data1 から data4 は既に７次元配列なので、それらをまとめる事はできないです。 335 名前：デフォルトの名無しさん mailto:sage [2009/06/18(木) 19:22:23 ] >>333 普通は配列をdata(n,4)と2次元にして do i=1,4 call hoge(data(:,i)) enddo 336 名前：デフォルトの名無しさん mailto:sage [2009/06/18(木) 19:25:49 ] >>334 7次元配列ってすごいなｗ Fortranの配列って7が限界だったっけ？それなら無理じゃないかな あきらめれ 337 名前：デフォルトの名無しさん mailto:sage [2009/06/18(木) 23:01:17 ] >>332 新しいコンパイラでは、どれを書いても最適化で同じになるらしいが、 初期のFortran90コンパイラでは、n：nnみたいな指定子が入ると、 テンポラリをとることが多かったみたい。 それは多分一般的には、n：nn：mみたいな不連続の番地をとる可能性があるためだと思う。 >>333 ポインターの配列を、TYPEによって定義しておいて、それに配列を割り付けておいて、 それを引数に置けばいいんじゃないかな？ まぁポインターに割り当てるときに、だらだらとdata1~4を並べないと駄目なので ゴミの位置が移動するだけだがｗ 338 名前：337 mailto:sage [2009/06/19(金) 01:27:21 ] MODULE m_billy IMPLICIT NONE INTEGER, PARAMETER :: nmax = 2 TYPE :: t_billy INTEGER, POINTER :: mm(:, :, :, :, :, :, :) END TYPE t_billy CONTAINS SUBROUTINE anal(mm) INTEGER, INTENT(OUT) :: mm(:, :, :, :, :, :, :) INTEGER :: i1, i2, i3, i4, i5, i6, i7 DO i1 = 1, nmax DO i2 = 1, nmax DO i3 = 1, nmax DO i4 = 1, nmax DO i5 = 1, nmax DO i6 = 1, nmax DO i7 = 1, nmax mm(i1, i2, i3, i4, i5, i6, i7) = i1 + 10 * (i2 + 10 * (i3 + 10 * (i4 + 10 * (i5 + 10 * (i6 + 10 * i7))))) END DO END DO END DO END DO END DO END DO END DO RETURN END SUBROUTINE anal END MODULE m_billy 339 名前：337 mailto:sage [2009/06/19(金) 01:28:03 ] !==================================================== PROGRAM herrington USE m_billy IMPLICIT NONE INTEGER, TARGET :: data1(nmax, nmax, nmax, nmax, nmax, nmax, nmax), data2(nmax, nmax, nmax, nmax, nmax, nmax, nmax) INTEGER, TARGET :: data3(nmax, nmax, nmax, nmax, nmax, nmax, nmax), data4(nmax, nmax, nmax, nmax, nmax, nmax, nmax) TYPE (t_billy) :: kazuya(4) INTEGER :: i kazuya(1)%mm => data1 kazuya(2)%mm => data2 kazuya(3)%mm => data3 kazuya(4)%mm => data4 DO i = 1, 4 CALL anal( kazuya(i)%mm ) END DO print *, data1 STOP END PROGRAM herrington 340 名前：333 [2009/06/19(金) 01:58:08 ] サンクス。でも結構大掛かりだな。。。 real*8 data1(n,n,n,n,n,n,n), data2(n,n,n,n,n,n,n),....... char*1 cnum char*5 target do i = 1, 4 write( cnum, '(i1)' ) i target = "data" // cnum call hoge( target ) end do みたいにできんもんかね。。。 target と data$ は型が違うがな。 341 名前：337 mailto:sage [2009/06/19(金) 02:24:46 ] >>340 77の掟破りでよければ、 real*8 data1(n,n,n,n,n,n,n), data2(n,n,n,n,n,n,n), data3(n,n,n,n,n,n,n), data4(n,n,n,n,n,n,n) real*8 tmp(n**7, 4) equivalence (tmp(1, 1), data1), (tmp(1, 2), data2), (tmp(1, 3), data3), (tmp(1, 4), data4) do i = 1, 4 call hoge( tmp(i) ) end do こんな感じでEUIVALENCEでいけるかな？ 342 名前：デフォルトの名無しさん mailto:sage [2009/06/19(金) 03:00:32 ] >>341 call hoge( tmp(1, i) ) の過ち。 ハルヒ新OP糞杉ｗｗｗ 343 名前：デフォルトの名無しさん mailto:sage [2009/06/19(金) 10:22:30 ] >>334 本当に7次元いるの？　変数設計を再検討した方がいいんじゃない？ 座標3次元＋時間で、4次元はざらにあるが。 344 名前：デフォルトの名無しさん mailto:sage [2009/06/19(金) 14:06:33 ] >>343の言うとおり。 変数設計は、アルゴリズム的な自然さ（人間にとっての見やすさ）と速度の 両方の観点から決めるべき。 345 名前：デフォルトの名無しさん mailto:sage [2009/06/20(土) 20:01:44 ] 7次元配列の数値計算ってどういう分野で出てくるの？ 346 名前：デフォルトの名無しさん mailto:sage [2009/06/20(土) 21:51:11 ] 流れぶった切ってすいません。 n×n 次元行列　M と　n 次元ベクトルa の積 b = Ma を計算するプログラムを 作っていたのですが、詰まってしまいました。 どこが、おかしいか教えていただけないでしょうか？ program integer ::i,k,n integer,parameter ::limit=5 real ::M(1:limit,1:limit),a(1:limit),b(1:limit) do print *,'nの値は？' read *,n if ((n<=limit).and.(n>=1)) then exit end if end do ! print *,'行列Mの成分' do i=1,n print *,i,'行目?' read *,(M(i,k),k=1,n) end do print *,'ベクトルaの成分' do k=1,n read *,a(k) end do 347 名前：346 mailto:sage [2009/06/20(土) 21:52:41 ] !>>346の続き do i=1,n do k=1,n b(0)=0 b(k)=b(k)+M(i,k)*a(k) end do end do print *,'b:' do k=1,n print *,b(k) end do stop end program よろしくお願いします 348 名前：346 mailto:sage [2009/06/20(土) 21:57:18 ] たびたびすいません。 上のプログラムはfortran90を使っています。 コンパイルする際は、ifort (プログラム名）.f90 -o (プログラム名） としています。 349 名前：デフォルトの名無しさん mailto:sage [2009/06/20(土) 21:59:30 ] >>346 全く見ていないが、b(0)=0　がおかしい。 b（k）=0.0 でかつ一個DOLOOPの外。 350 名前：デフォルトの名無しさん mailto:sage [2009/06/21(日) 10:45:18 ] ループのところ do i=1, n b(i) = 0 do k=1,n b(i) = b(i) +M(i, k) * a(k) end do end do 直接関係ないけど、Mの行と列を入れかえたほうがいいかも。 351 名前：346 mailto:sage [2009/06/21(日) 18:17:49 ] >>350さん ありがとうございました。 うまくいきました！ 352 名前：デフォルトの名無しさん mailto:sage [2009/06/23(火) 11:08:29 ] PROGRAM GAUSS_SOLVER IMPLICIT NONE REAL,DIMENSION(100,100)::A !行列A（最大100元まで解ける） REAL,DIMENSION(100)::B,X !右辺Bと解Xのベクトル INTEGER::N !元数 INTEGER::I,J,K !制御変数 WRITE(*,*)'このプログラムは連立一次方程式AX=Bのい解Xを求めます' WRITE(*,*)'AandB.TXTの入力フォーマットは次のようです' WRITE(*,*)'たとえば,３元の連立方程式の場合' WRITE(*,*)'3 (元数を与えてください)' WRITE(*,*)'A11 A12 A13 B1 (データは空白やタブで別けてください)' WRITE(*,*)'A21 A22 A23 B2 (AIJとBIには実数を与えてください)' WRITE(*,*)'A31 A32 A33 B3' WRITE(*,*)'解XはX.TXTに出力されます' OPEN(1,FILE='AandB.TXT') READ(1,*) N DO I=1,N READ(1,*) (A(I,J),J=1,N),B(I) END DO CLOSE(1) DO I=1,N-1 DO J=I+1,N DO K=I+1,N A(J,K)=A(J,K)-A(J,I)*A(I,K)/A(I,I) END DO B(J)=B(J)-B(I)*A(J,I)/A(I,I) END DO END DO 353 名前：デフォルトの名無しさん mailto:sage [2009/06/23(火) 11:10:04 ] DO I=N,1,-1 X(I)=B(I)/A(I,I) DO J=I+1,N X(I)=X(I)-A(I,J)*X(J)/A(I,I) END DO END DO OPEN(1,FILE='X.TXT') DO I=1,N WRITE(1,*) 'X(',I,')=',X(I) END DO CLOSE(1) STOP END PROGRAM GAUSS_SOLVER 354 名前：デフォルトの名無しさん mailto:sage [2009/06/23(火) 11:20:09 ] このプログラムをもとにN次行列の逆行列をもとめるプログラムを作るのですが どうしたらいいかわかりません。教えてくださると助かります。。 355 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 00:59:41 ] ガウス・ジョルダン消去法でぐぐれ ちなみに逆行列を求める時は元の配列に逆行列を 入れる事が出来て効率が良い 356 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:19:20 ] subroutine gauss_jordan(a,b,flag) implicit none integer,parameter::N=3 double precision ::a(N,N),b(N) integer ipv,i,j double precision :: inv_pivot,temp double precision :: big integer pivot_row,row(N) integer flag do ipv=1,N !ｺﾇﾂ酖ﾍﾃｵｺ・ big=0.0 do i=ipv,N if(abs(a(ipv,i)) > big) then big = abs(a(ipv,i)) pivot_row = i endif enddo if(abs(big) < 0.0001) then write(6,*) "A is singular matrix" flag = 0 return endif 357 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:20:34 ] row(ipv) = pivot_row !ｹﾔ､ﾎﾆ�､・ﾘ､ｨ if(ipv .ne. pivot_row) then do i=1,N temp = a(i,ipv) a(i,ipv) = a(i,pivot_row) a(i,pivot_row) = temp enddo temp = b(ipv) b(ipv) = b(pivot_row) b(pivot_row) = temp endif ! ﾂﾐｳﾑﾀｮﾊｬ=1(･ﾔ･ﾜ･ﾃ･ﾈｹﾔ､ﾎｽ靉�) inv_pivot = 1.0/a(ipv,ipv) a(ipv,ipv) = 1.0 358 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:21:30 ] do j=1,N a(j,ipv) = a(j,ipv) * inv_pivot enddo b(ipv) = b(ipv) * inv_pivot ! ･ﾔ･ﾜ･ﾃ･ﾈﾎ・0(･ﾔ･ﾜ･ﾃ･ﾈｹﾔｰﾊｳｰ､ﾎｽ靉�) do i=1,N if(i .ne. ipv) then temp = a(ipv,i) a(ipv,i) = 0.0 do j=1,N a(j,i) = a(j,i) - temp*a(j,ipv) enddo b(i) = b(i) - temp*b(ipv) endif enddo 359 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:23:23 ] enddo ! ﾎ�､ﾎﾆ�､・ﾘ､ｨ(ｵﾕｹﾔﾎ・ do j=N,1,-1 if(j .ne. row(j)) then do i=1,N temp = a(j,i) a(j,i) = a(row(j),i) a(row(j),i) = temp enddo endif enddo flag = 1 return end subroutine !### TEST MAIN ###### ! program MAIN ! implicit none ! integer,parameter::N=3 ! double precision :: a(N,N),b(N) ! integer :: flag 360 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:25:41 ] ! a(1,1) = 0.0 ! a(2,1) = 1.0 ! a(3,1) = 0.0 ! a(1,2) = 1.0 ! a(2,2) = 0.0 ! a(3,2) = 0.0 ! a(1,3) = 0.0 ! a(2,3) = 0.0 ! a(3,3) = 1.0 ! ! b(1) = 3.0 ! b(2) = 4.0 ! b(3) = 5.0 ! call gauss_jordan(a,b,flag) ! write(6,*) b(1),b(2),b(3) ! stop ! end program 361 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:27:01 ] UTF-8か何かで貼ったのか？ 文字化けしてて日本語部分読めないぞ 362 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:27:21 ] これのことですか？長々とすみません 363 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:29:24 ] ダウンロードしたものをノートパッドで開いてコピペしたんですが。。 364 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:33:53 ] で、わかったっしょ？ガウスジョルダン消去法を使うと逆行列が 求められるという事が もちろん行列式が零だと潰れてしまうが 365 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:40:14 ] 半分くらい・・・（汗 366 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:42:15 ] Cのソースで良いなら nyan11.ciao.jp/B/NA/EWSNA/linkANNA.html のINV1.Cが逆行列を実際にガウスジョルダン消去法で解くプログラムだ これを参考にFORTRANに移植するとよい 367 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:45:03 ] Intel Fortran v11.1 が出たようだな。 まだ完全ではないようだが、ようやくF2003時代がやってきたｗ 368 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:49:38 ] まだ分からないと言われそうなんで一応逃げを打っておく www.sic.shibaura-it.ac.jp/~i020201/programming_a/lecture_2005/lecture_09.pdf ここのPDFをじっくり眺めてくれ さすがに理解できるだろ 369 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 22:58:35 ] 去年も逆行列の質問をしていた奴がいたが、さては留年したかｗ 今年はがんばれ。 質問でうまくおだてて解答ゲットしろｗ 370 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 23:05:27 ] これでも分からないならネットに頼るのはやめて明日本屋に行け そして線形代数の本を買え 掃き出し法という名前で載っているかもしれないがとにかくそれだ 手で一回3×3程度の行列の逆行列を求めてみればおのずと プログラムでやっている事も理解できる 371 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 23:08:37 ] detail.chiebukuro.yahoo.co.jp/qa/question_detail/q1418060584 これとかな いくらでもネット上にも情報はあるんだけどな 372 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 23:25:05 ] 数学的なやり方そのものに疑問はないです・・・ ただプログラムできないです・・・ちなみにその去年の人じゃないよ 373 名前：デフォルトの名無しさん mailto:sage [2009/06/24(水) 23:33:18 ] >>372 おｋ。 疑問点を具体的に書いてくれれば、具体的に答えられるが、 漠然と聞かれると、漠然としか答えられない。 まぁがんばれｗ 374 名前：デフォルトの名無しさん mailto:sage [2009/06/25(木) 04:52:56 ] 最初は誰でも、どこがわからないかわからない状態だよね まあがんばれw 375 名前：デフォルトの名無しさん mailto:sage [2009/06/25(木) 14:49:13 ] >>372 というか、質問の仕方というものがあるだろう・・・。 いきなりソースコースばしばし貼りまくってから質問内容を言うのは、分かりにくいというか それ以前のマナーの問題。せめて要点をまとめろよ。 376 名前：デフォルトの名無しさん mailto:sage [2009/06/25(木) 15:24:15 ] いいやん こっちも動くソースをばしばし貼って、それで終わりにして 質問は一切受け付けませんよという態度を取れば 377 名前：デフォルトの名無しさん mailto:sage [2009/06/25(木) 15:52:01 ] そしたらスレが汚くなるじゃないか。質問者はソースを貼るときは、必要でない限り ・ソース丸ごとではなく、最小限の部位に留める ・変数名を見やすい、もしくは一般的なものに変える ぐらいはやってくれよと。こんなの常識じゃないかと思うんだが。 378 名前：デフォルトの名無しさん mailto:sage [2009/06/25(木) 20:32:26 ] >>377　それができたら初心者じゃないから、まあ、大目にみようよ 379 名前：デフォルトの名無しさん mailto:sage [2009/06/25(木) 22:27:34 ] >>372ですが今日一日教科書じっくり読み込んでだいたいわかりました ぐだぐだとプログラムを貼り付けてすいませんでした。。 あと、教えてくれた皆さんありがとうございます 380 名前：デフォルトの名無しさん mailto:sage [2009/06/25(木) 23:03:19 ] >>379 そうかそれはよかったな 君のした努力は決して無駄にならないと思うよ これから線形代数・微積分学・複素関数論とかいろんな 数学をプログラミングする機会があると思うが、今した努力は 必ず生きてくる 特にFFTのプログラムなんか芸術的だぞ バタフライ演算って言ってな 381 名前：デフォルトの名無しさん mailto:sage [2009/06/25(木) 23:07:12 ] >>380　 ありがとうございます 読んでる内にどんどん楽しくなってきて、一気に今やってるところまで 読み進めれました。 これからもっとがんばります（｀・ω・´） 382 名前：デフォルトの名無しさん mailto:sage [2009/06/25(木) 23:20:35 ] 立派なフォートランナーになってくれ^^ 383 名前：デフォルトの名無しさん mailto:sage [2009/06/25(木) 23:44:58 ] >>377 おまえの常識なんか初心者はしったこっちゃねーんだよｗ 384 名前：デフォルトの名無しさん mailto:age [2009/06/27(土) 15:04:27 ] C IMPLICIT REAL(A-H,O-Z) REAL T(100),XZ(100,2) DATA GRAV,V0,ANGLE/9.8,30.0,45.0/ DATA T0,VINT,TMAX/0.0,0.5,7.5/ PI=ATAN(1.0)*4.0 R=ANGLE*PI/180.0 C DO 10 I=1,50,1 T(I)=REAL(I-1)*VINT IF(T(I).GT.TMAX)THEN STOP END IF XZ(I,1)=V0*COS(R)*T(I) XZ(I,2)=-0.5*GRAV*T(I)**2+V0*SIN(R)*T(I) 10 CONTINUE S=AREA(T,X,Z,T0,TMAX,R,GRAV) CALL OUTPUT(T,X,Z,S) STOP END 385 名前：デフォルトの名無しさん mailto:sage [2009/06/27(土) 15:06:47 ] C FUNCTION AREA (T,X,Z,T0,TMAX,R,GRAV) IMPLICIT REAL(A-H,O-Z) XMAX=V0*0.5*COS(R)*TMAX X0=V0*0.5*COS(R)*T0 DO 20 I=1,16 AREA=0.0 DX=(XMAX-X0)/REAL(I) T=REAL(I-1)*VINT X1=V0*COS(R)*T X2=X1+DX Z1=TAN(R)*X1-0.5*(GRAV/(V0*(COS(R))**2))*X1**2 Z2=TAN(R)*X2-0.5*(GRAV/(V0*(COS(R))**2))*X2**2 AREA=AREA+(Z1+Z2)*DX*0.5 20 CONTINUE RETURN END 386 名前：デフォルトの名無しさん mailto:sage [2009/06/27(土) 15:09:07 ] C SUBROUTINE OUTPUT(T,X,Z,S) IMPLICIT REAL(A-H,O-Z) INTEGER NO(100) REAL T(100),XZ(100,2) C OPEN(16,FILE='menseki2.res') WRITE(16,'(A)')' TIME CX CZ SPACE' DO 30 I=1,16 WRITE(16,'(F5.1,2F10.2)') T(I),(XZ(I,J),J=1,2) WRITE(16,'(F10.2)') S 30 CONTINUE CLOSE(16,STATUS='KEEP') RETURN END 物体の投げあげの時間、X座標とZ座標、台形積分の計算をファイルに出力したいのですが、コンパイルしてもファイルが作られません。 ファイルは時間が１次元配列、X座標とZ座標が２次元配列です。 どう改善すればよいのでしょうか 387 名前：デフォルトの名無しさん mailto:sage [2009/06/27(土) 15:45:54 ] >>386 T(I) が TMAX (=7.5) を超えたところでSTOPするようになってるから。 てゆーかぱっと見、変数の使われ方が変だ。 XZ は宣言されてるが X とか Z が突然出てきてるぞ。 IMPLICIT NONEで動くところまで書き直せ。 388 名前：デフォルトの名無しさん mailto:sage [2009/06/27(土) 17:50:39 ] >>384-386 あんた>>352-354か？ ソースを貼るなとは言わんが、まず「○○で困っています。××はどうすればいいでしょうか？ソースは以下の通りです」とか書いてから貼れよ 389 名前：デフォルトの名無しさん mailto:sage [2009/06/27(土) 17:57:12 ] プログラムの書き型が相当古いから別人だろうな。 390 名前：デフォルトの名無しさん mailto:sage [2009/06/27(土) 22:05:47 ] >>387 それ以前に、コンパイルが通らないんだろ。 AREAの引数のTが、スカラー定義なのに配列を渡している。 書き方は古いが、非常にオーソドックスで正当な書き方をしている。 しかし、実際の論理というか中身はかなり？？ 簡単には直らんｗ 391 名前：１／２ mailto:sage [2009/06/27(土) 22:59:57 ] >>386 IMPLICIT REAL(A-H,O-Z) REAL T(100),XZ(100,2) DATA GRAV,V0,ANGLE/9.8,30.0,45.0/ DATA T0,TINT,TMAX/0.0,0.5,7.5/ PI=ATAN(1.0)*4.0 R=ANGLE*PI/180.0 C C ! IMAX <= 100 no check www IMAX = INT((TMAX - T0) / TINT) + 1 DO 10 I=1,IMAX T(I)=T0 + REAL(I-1)*TINT XZ(I,1)=V0*COS(R)*T(I) XZ(I,2)=-0.5*GRAV*T(I)**2+V0*SIN(R)*T(I) 10 CONTINUE S=AREA(IMAX, XZ) CALL OUTPUT(IMAX, T,XZ,S) STOP END C FUNCTION AREA(IMAX, XZ) IMPLICIT REAL(A-H,O-Z) REAL XZ(100, 2) AREA=0.0 DO 20 I=1, IMAX - 1 AREA=AREA+( XZ(I,2)+XZ(I+1,2) )*( XZ(I+1,1)-XZ(I,1) )*0.5 20 CONTINUE RETURN END 392 名前：２／２ mailto:sage [2009/06/27(土) 23:02:17 ] C SUBROUTINE OUTPUT(IMAX,T,XZ,S) IMPLICIT REAL(A-H,O-Z) INTEGER NO(100) REAL T(100), XZ(100,2) C OPEN(16,FILE='menseki2.res') WRITE(16,'(A)')' TIME CX CZ SPACE' DO 30 I=1,IMAX WRITE(16,'(F5.1,2F10.2)') T(I),(XZ(I,J),J=1,2) 30 CONTINUE WRITE(16,'(A, F10.2)') 'AREA=', S CLOSE(16,STATUS='KEEP') RETURN END イマイチよく分からんが、こうかな？ とりあえず放物線のようなものは描く。 しかし、地面をつきぬけていくので全積分面積は負になる。 この辺は適宜修正してくれ。

---

---

[ 続きを読む ] / [ 携帯版 ]

read.cgi ver5.27 [feat.BBS2 +1.6] / e.0.2 (02/09/03) / eucaly.net products.
担当:undef