zyj8119

木蟲 (著名寫手)

應助: 65 (初中生)
貴賓: 0.003
金幣: 915.1
散金: 1440
紅花: 35
帖子: 2936
在線: 1329.4小時
蟲號: 664177
注冊: 2008-11-29
性別: GG
專業(yè): 理論和計算化學

[求助] 一個FORTRAN90程序，越界了，希望高手幫助

CODE:

module inv_mat

contains

subroutine inv(A,invA,N)

implicit real*8(a-z)

integer::n

integer::i

real*8::A(n,n),invA(n,n),E(n,n)

E=0

do i=1,n

   E(i,i)=1

end do

call mateq(A,B,X,N,N)

end subroutine inv

subroutine mateq(A,B,X,N,M)

implicit real*8(a-z)

integer::N,M,i

real*8::A(M,N),B(N,M),X(N,M)

real*8::btemp(N),xtemp(N)

do i=1,M

   btemp=B(:,i)

   call elgauss(A,btemp,xtemp,N)

   X(:,i)=xtemp

end do

end subroutine mateq

subroutine elgauss(A,b,x,N)

implicit real*8(a-z)

integer::i,k,N

integer::id_max

real*8::A(N,N),b(N),x(N)

real*8::Aup(N,N),bup(N)

real*8::Ab(N,N+1)

real*8::vtemp1(N+1),vtemp2(N+1)

Ab(1:N,1:N)=A

Ab(:,N+1)=b

do k=1,N-1

   elmax=dabs(Ab(k,k))

   id_max=k

   do i=k+1,n

    if(dabs(Ab(i,k))>elmax)then

        elmax=Ab(i,k)

        id_max=i

        end if

        end do

        vtemp1=Ab(k,:)

        vtemp2=Ab(id_max,:)

        Ab(k,:)=vtemp2

        Ab(id_max,:)=vtemp1

 do i=k+1,N

    temp=Ab(i,k)/Ab(k,k)

        Ab(i,:)=Ab(i,:)-temp*Ab(k,:)

 end do

end do

Aup(:,:)=Ab(1:N,1:N)

bup(:)=Ab(:,N+1)

call uptri(Aup,bup,x,n)

end subroutine elgauss

subroutine uptri(A,b,x,N)

implicit real*8(a-z)

integer::i,j,N

real*8::A(N,N),b(N),x(N)

x(N)=b(N)/A(N,N)

do i=n-1,1,-1

    x(i)=b(i)

        do j=i+1,N

        x(i)=x(i)-a(i,j)*x(j)

        end do

        x(i)=x(i)/A(i,i)

end do

end subroutine uptri

end module inv_mat

module m_bfs

use inv_mat

contains

subroutine solve(x0,N,tol)

implicit real*8(a-z)

integer::i,n,itmax=50

real*8::x1(n),x0(n),y(n),f0(n),f1(n),dx(n)

real*8::v1(n),v2(n),v3(n)

real*8::H0(n,n),H1(n,n),df(n,n)

real*8::m1(n,n),m2(n,n),m3(n,n)

call jac(df,x0)

call inv(df,H0,N)

write(11,101)

write(12,102)

do i=1,itmax

   call func(f0,x0)

   x1=x0-matmul(H0,f0)

   dx=x1-x0

   call func(f1,x1)

   y=f1-f0

   v1=matmul(H0,y)

   t1=vdot(y,v1,n)

   t2=vdot(dx,y,N)

   u=1d0+t1/t2

   m1=vvmat(dx,dx,n)*u

   m2=vvmat(dx,y,n)

   m2=matmul(m2,H0)

   v3=matmul(H0,y)

   m3=vvmat(v3,dx,n)

   t3=vdot(dx,y,n)

   H1=(m1-m2-m3)/t3

   H1=H0+H1

   write(12,103)i,H0

   x0x=1

   H0=H1

   write(11,104)i,x0

   dx2=dsqrt(dx(1)**2+dx(2)**2)

   if(dx2
end do

101  format(/,T5,'BFS方法計算序列',/)

102  format(/,T5,'BFS方法H矩陣序列為:',/)

103  format(T5,'iter=',I4,/,(F16.10/),/)

104  format(I8,3F16.10)

end subroutine solve

function vdot(a,b,N)

integer::i,N

real*8::a(n),b(n),vdot

vdot=0

do i=1,N

   vdot=vdot+a(i)*b(i)

end do

end function vdot

function vvmat(a,b,n)

integer::n,i,j

real*8::a(n),b(n),vvmat(n,n)

do i=1,n

   do j=1,n

     vvmat(i,j)=a(i)*b(j)

   end do

end do

end function  vvmat

subroutine func(f,x)

implicit real*8(a-z)

real*8::x(3),f(3),pi=3.14159263589793

f(1)=3*x(1)-dcos(x(2)*x(3))-1d0/2

f(2)=x(1)**2-81*(x(2)+0.1d0)**2+dsin(x(3))+1.06d0

f(3)=dexp(-x(1)*x(2))+20*x(3)+(10*pi-3d0)/3

end subroutine func

subroutine jac(df,x)

implicit real*8(a-z)

real*8::x(3),df(3,3)

df(1,1)=3d0

df(1,2)=x(3)*dsin(x(2)*x(3))

df(1,3)=x(2)*dsin(x(2)*x(3))

df(2,1)=2*x(1)

df(2,2)=-162*(x(1)+0.1)

df(2,3)=dcos(x(3))

df(3,1)=-x(2)*dexp(-x(1)*x(2))

df(3,2)=-x(1)*dexp(-x(1)*x(2))

df(3,3)=20d0

end subroutine jac

end module m_bfs

program main

use inv_mat

use m_bfs

implicit real*8(a-z)

integer::N=3

real*8::x0(3)

open(unit=11,file='result.txt')

open(unit=12,file='Hmatrix.txt')

x0=(/0.1d0,0.1d0,-0.1d0/)

call solve(x0,n,1d-8)

end program main

回復此樓

» 猜你喜歡

» 本主題相關價值貼推薦，對您同樣有幫助:

請教一個應用read時出現(xiàn)數(shù)組越界的問題，謝謝已經有7人回復
有會看fortran程序的，能幫忙看看這個程序嗎？已經有10人回復
怎么使用fortran程序創(chuàng)建一個立方晶格模型已經有5人回復
如何改變fortran的浮點精度已經有28人回復
求助高手幫忙編一個Fortran的小程序，重金懸賞...... 已經有15人回復
求助計算機高人，編寫個Fortran程序已經有13人回復
請教一個fortran小程序編譯出錯的問題，謝謝已經有9人回復
寫了一個fortran90的小程序，編譯通不過，請大俠幫忙已經有59人回復
【求助】如何寫FORTRAN程序實現(xiàn)求平均最近鄰距離已經有34人回復
【求助】求高手關于Fortran數(shù)據(jù)讀寫的問題�！疽淹杲Y】已經有12人回復
【求助】請問fortran90與95區(qū)別大嗎已經有4人回復
【求助】此FORTRAN90程序沒有錯誤，但是運行不出來。已經有4人回復
【求助】求助個Fortran計算，應該不難，半天搞不定，請高手幫忙已經有11人回復

好好學習，天天向上。

1樓 2011-11-11 23:00:16

已閱回復此樓關注TA 給TA發(fā)消息送TA紅花 TA的回帖

zyj8119

木蟲 (著名寫手)

應助: 65 (初中生)
貴賓: 0.003
金幣: 915.1
散金: 1440
紅花: 35
帖子: 2936
在線: 1329.4小時
蟲號: 664177
注冊: 2008-11-29
性別: GG
專業(yè): 理論和計算化學

引用回帖:

2樓: Originally posted by maomao1210 at 2011-11-11 23:15:23:
大哥；
子程序 inv, X是數(shù)組，你忘記定義了。X（n,n）

謝謝提醒，修改過來了，編譯通過。。

贊一下

回復此樓

好好學習，天天向上。

3樓2011-11-11 23:29:18

已閱回復此樓關注TA 給TA發(fā)消息送TA紅花 TA的回帖

查看全部 6 個回答

maomao1210

金蟲 (正式寫手)

程序強帖: 5
應助: 2 (幼兒園)
金幣: 1432.8
散金: 242
紅花: 16
沙發(fā): 1
帖子: 991
在線: 442.1小時
蟲號: 253215
注冊: 2006-05-20
性別: MM
專業(yè): 考古理論

★
jjdg(金幣+1): 感謝參與 2011-11-12 00:40:40
xzhdty: 2011-11-12 22:51:02
ben_ladeng: 專家考核存檔 2011-11-12 23:23:07

引用回帖:

1樓: Originally posted by zyj8119 at 2011-11-11 23:00:16:
[code]module inv_mat
contains
subroutine inv(A,invA,N)
implicit real*8(a-z)
integer::n
integer::i
real*8::A(n,n),invA(n,n),E(n,n)
E=0
do i=1,n
E(i,i)=1
end do
call mateq(A,B,X,N,N)
e ...

大哥；
子程序 inv, X是數(shù)組，你忘記定義了。X（n,n）

贊一下(3人)

回復此樓

2樓2011-11-11 23:15:23

已閱回復此樓關注TA 給TA發(fā)消息送TA紅花 TA的回帖

snoopyzhao

至尊木蟲 (職業(yè)作家)

程序強帖: 16
應助: 157 (高中生)
貴賓: 0.02
金幣: 18844.7
紅花: 29
帖子: 3803
在線: 1422.4小時
蟲號: 183750
注冊: 2006-02-13
專業(yè): 污染生態(tài)化學

【答案】應助回帖

★
dubo(金幣+1): 多謝應助 2011-11-12 12:26:34
zyj8119(金幣+10): 建議很好，會認真考慮撒。。。 2011-11-12 14:04:34
ben_ladeng: 專家考核存檔 2011-11-12 23:23:16

引用回帖:

3樓: Originally posted by zyj8119 at 2011-11-11 23:29:18:
謝謝提醒，修改過來了，編譯通過。。

從這個錯誤來看，一定要用 implicit none……

另外，盡可能不要用 real*4、real*8 這些編譯器擴展的用法

贊一下(2人)

回復此樓

4樓2011-11-12 12:13:01

已閱回復此樓關注TA 給TA發(fā)消息送TA紅花 TA的回帖

exabyss916

新蟲 (小有名氣)

應助: 1 (幼兒園)
金幣: 177
散金: 65
帖子: 78
在線: 57小時
蟲號: 507882
注冊: 2008-02-20
專業(yè): 理論和計算化學

引用回帖:

4樓: Originally posted by snoopyzhao at 2011-11-12 12:13:01:
從這個錯誤來看，一定要用 implicit none……

另外，盡可能不要用 real*4、real*8 這些編譯器擴展的用法

能不能簡單說一下為什么盡可能不要用 real*4、real*8 這些編譯器擴展的用法？

贊一下

回復此樓

5樓2011-11-19 10:45:59

已閱回復此樓關注TA 給TA發(fā)消息送TA紅花 TA的回帖

最具人氣熱帖推薦 [查看全部]	作者	回/看	最后發(fā)表

[考研] 266分，一志愿電氣工程，本科材料，求材料專業(yè)調劑 +10	哇呼哼呼哼 2026-04-01	11/550	2026-04-02 11:31 by lnilvy
[考研] 計算機265跨調環(huán)境 +5	Yumeng_6 2026-03-27	5/250	2026-04-02 10:54 by guanxin1001
[考研] 324求調劑 +4	想上學求調 2026-04-02	4/200	2026-04-02 10:00 by olim
[考研] 各位老師好，我的一志愿為北京科技大學085601材料專碩 +13	Koxui 2026-03-28	13/650	2026-04-02 09:35 by 噠噠噠呱呱呱
[考研] 求生物學調劑 +8	15172915737 2026-04-01	8/400	2026-04-02 06:49 by ilovexiaobin
[考研] 0856，269分求調劑 +8	有學上就行求求?/a> 2026-03-30	11/550	2026-04-01 22:33 by 2026材料調劑
[考研] 332求調劑 +8	Lyy930824@ 2026-03-29	8/400	2026-04-01 18:40 by 千鈞澄玉宇啊
[考研] 353求調劑 +4	拉鉤不許變 2026-04-01	4/200	2026-04-01 18:10 by 記事本2026
[考研] 0710生物學考研調劑 +3	李多米lee. 2026-03-27	4/200	2026-04-01 16:21 by zzchen2000
[考研] 08工科275求調劑，可跨考。 +5	AaAa7420 2026-03-31	5/250	2026-04-01 15:21 by 159357hjz
[考研] 一志愿鄭大085600，310分求調劑 +6	李瀟可 2026-03-26	6/300	2026-04-01 14:44 by chenqifeng666
[考研] 一志愿中國科學院大學265求調劑 +8	恬淡ye 2026-03-31	9/450	2026-04-01 14:34 by 逆水乘風
[考研] 0855機械初試280求調劑 +3	kazenotori 2026-03-31	3/150	2026-04-01 10:08 by JourneyLucky
[考研] 本2一志愿C9-333分，材料科學與工程，求調劑 +9	升升不降 2026-03-31	9/450	2026-03-31 18:01 by 無際的草原
[考研] 262求調劑 +7	ZZ..000 2026-03-30	8/400	2026-03-31 10:05 by cal0306
[考研] 本科211總分289，08工學真心求調劑 +3	utopiaE 2026-03-30	3/150	2026-03-30 23:42 by ms629
[考研] 342求調劑 +4	加油a李zs 2026-03-26	4/200	2026-03-30 16:39 by 晶體之美
[考研] 一志愿中南大學化學0703總分337求調劑 +6	niko- 2026-03-27	6/300	2026-03-30 10:25 by herarysara
[考研] 081200-314 +3	LILIQQ 2026-03-27	4/200	2026-03-28 09:41 by 保護地球你我做?/a>
[考研] 305求調劑 +5	哇盧卡庫 2026-03-26	5/250	2026-03-27 14:01 by laoshidan