| 查看: 434 | 回复: 2 | ||
songjunann铜虫 (小有名气)
|
[求助]
求助一个调用程序的问题，谢谢
|
|
如下程序，能否在主程序里面这么调用： Call slvasi(a,b,b,idiag, neq, 1,ierr) 我觉得，有两个b很有问题，但是别人给我的一个大程序就是这么干的，我也不知道如何在大程序里面改，先问问这种调用有问题吗（如没有，我也就放过去了，呵呵）？ 非常感谢 subroutine slvasi (a,b,x,mdiag,neq,key,ierr) cc-begin c----------------------------------------------------------------------- c c call slvasi (a,b,x,mdiag,neq,key,ierr) c c usage : in-core non-symmetric system c ===== c c c purpose c ------- c c to solve a linear system of equations [a] [x] = , where a c is total stiffness matrix c cc-doc c input arguments c --------------- c c a = vector containing total stiffness matrix. c b = load vector of length neq (see note 1). c mdiag = vector of length neq which contains the position c of the diagonal element in the total stiffness c matrix.(outputted by diag) c neq = number of equations. c key = operation control key : c c key = 1 perform lu factorisation of a and solve ax = b. c key = 2 perform lu factorisation of a only. c key = 3 reduction of load vector b and back c substitution to obtain solution x (see note 2). c output arguments c ---------------- c c x = solution vector of length neq (see note 1). c ierr = error flag (see note 5). c c reference c --------- c c LSD / FEM Version 2. c A Library for Software Development in the Finite Element Method c Technical Report 85-3, Department of Computer Science c The Hebrew University of Jerusalem. c c notes c ----- c c 1. the vectors b and x can occupy the same storage space i.e. on c output b will contain the solution. c c 2. if key = 3 , slvasi must have previously been called with c key = 1 or 2. slvasi can be called repeatedly with key = 3, c providing a, mdiag and length are not changed. c c 3. slvasi uses the library-supplied function dot(x,y,n) which c calculates the inner (dot,scalar) product of the vectors x,y c of length n. c c 4. the subroutine slvasi can be used independently of the package c presented here. the sole requirement of slvasi is the vectors c a,b,mdiag and length written in the described format. c c 5. ierr = 0 execution successfull c = n a zero pivot was encountered on row n during lu c factorisation. factorisation aborted. c c----------------------------------------------------------------------- cc-end c implicit real*8 (a-h,o-z) dimension a(*),b(*),x(*),mdiag(*) c c ierr = 0 if (key .eq. 3) go to 50 c c loop on rows of matrix c mdi = 0 do 40 i = 1,neq mdim = mdi mdi = mdiag(i) lengli= (mdi - mdim - 1)/2 lengui= mdi - 1 - mdim - lengli lowli = i - lengli lowui = i - lengui loci0 = mdi - lengui - i loc0i = mdi - i im = i - 1 ifj = min0 (lowui,lowli) ilj = im if (ifj .gt. ilj) go to 30 jj = ifj - 1 mdj = 0 if (ifj .gt. 1) mdj = mdiag(jj) c c loop on columns. c do 20 j = ifj,ilj mdjm = mdj mdj = mdiag(j) lenglj= (mdj - mdjm - 1)/2 lenguj= mdj - 1 - mdjm - lenglj lowlj = j - lenglj lowuj = j - lenguj if (lowli .gt. j) go to 10 c c compute l(i,j) c k = max0 (lowli,lowuj) lockj = mdj - j + k locij = loci0 + j locik = loci0 + k a(locij)= (a(locij) - dot(a(locik),a(lockj),j-k))/a(mdj) c c compute u(j,i) c 10 if (lowui .gt. j) go to 20 k = max0 (lowui,lowlj) locji = loc0i + j locki = loc0i + k locjk = mdj - lenguj - j + k a(locji) = a(locji) - dot(a(locjk),a(locki),j-k) 20 continue c c compute u(i,i) c 30 k = max0 (lowli,lowui) locik = loci0 + k locki = loc0i + k a(mdi) = a(mdi) - dot(a(locik),a(locki),i-k) if (a(mdi) .eq. 0.) go to 100 if (key .eq. 2) go to 40 c c reduce b if key = 1 c locik = loci0 + lowli x(i) = b(i) - dot (a(locik),x(lowli),i-lowli) 40 continue if (key .eq. 2) return go to 70 c c reduce b if key = 3 c 50 mdi = 0 do 60 i = 1,neq mdim = mdi mdi = mdiag(i) lengli= (mdi - mdim - 1)/2 lowli = i - lengli locik = mdim + 1 x(i) = b(i) - dot(a(locik),x(lowli),i-lowli) 60 continue c c back substitution c 70 ii = neq mdim = mdiag(neq) im = neq do 90 ll = 1,neq i = im im = im - 1 mdi = mdim ii = ii - 1 mdim = 0 if (i .ne. 1) mdim = mdiag(ii) lengli = (mdi - mdim - 1)/2 lengui = mdi - 1 - mdim - lengli lowui = i - lengui x(i) = x(i)/a(mdi) if (lowui .gt. im) go to 90 locij = mdi - lengui do 80 j = lowui,im x(j) = x(j) - a(locij)*x(i) locij= locij + 1 80 continue 90 continue return c c zero pivot c 100 ierr = i merf = 100*ierr call erfac (4hslva,merf) end |
回复此楼» 猜你喜欢
导师想让我从独立一作变成了共一第一
已经有9人回复
-
博士读完未来一定会好吗
已经有23人回复
-
到新单位后，换了新的研究方向，没有团队，持续积累2区以上论文，能申请到面上吗
已经有11人回复
-
读博
已经有4人回复
-
JMPT 期刊投稿流程
已经有4人回复
-
心脉受损
已经有5人回复
-
Springer期刊投稿求助
已经有4人回复
-
小论文投稿
已经有3人回复
-
申请2026年博士
已经有6人回复
» 本主题相关价值贴推荐，对您同样有帮助:
【求助】matlab 中，几个m文件调用
已经有14人回复- 【求助】mac的编译环境xcode的使用问题
已经有5人回复
- 【求助】请教一个关于解方程的问题
已经有8人回复
- 【求助】求指点，在Fortran里面，怎样实现“数组维度可调”的数组?
已经有5人回复
- 【求助】VC++MFC编程，对话框调用单文档，还要处理一些画图的东西
已经有16人回复
- 【求助】有很多错误的程序，找不出来了
已经有18人回复
- 【求助】求积分程序中被积函数问题
已经有15人回复
锐利的碎片
木虫 (正式写手)
star watcher
- 应助: 136 (高中生)
- 金币: 3637.1
- 散金: 252
- 红花: 22
- 帖子: 988
- 在线: 1224.9小时
- 虫号: 961933
- 注册: 2010-03-05
- 专业: 凝聚态物性 II ：电子结构
2楼2011-09-18 10:59:29
snoopyzhao
至尊木虫 (职业作家)
- 程序强帖: 16
- 应助: 157 (高中生)
- 贵宾: 0.02
- 金币: 18844.7
- 红花: 29
- 帖子: 3803
- 在线: 1422.4小时
- 虫号: 183750
- 注册: 2006-02-13
- 专业: 污染生态化学
【答案】应助回帖
★
dubo(金币+1): 欢迎讨论 2011-09-18 14:44:03
songjunann(金币+10): 恩，谢谢，我也注意到这个了，只是不确定 2011-09-18 14:57:30
dubo(金币+1): 欢迎讨论 2011-09-18 14:44:03
songjunann(金币+10): 恩，谢谢，我也注意到这个了，只是不确定 2011-09-18 14:57:30
|
这样做没有什么问题，如果这么调用后，返回值将使得 b 发生变化，从程序来看，这不会有什么问题。但是，如果你想 b 在调用这个子程后不要发生变化，就不能那么调用…… 看程序中，有一句话： c notes c ----- c c 1. the vectors b and x can occupy the same storage space i.e. on c output b will contain the solution. |
3楼2011-09-18 13:31:12