- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Whilst trouble shooting the program I stumbled across this code
SUBROUTINE HQRWT(N,M,G,E,V,A,B,W,ND,P,Q,XM,INT,ITAPE,NP,LF) c hqrw 2 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * hqrw 3 c subroutine to compute eigenvalues and eigenvectors of a hqrw 4 c symmetric real matrix stored in compact triangular form hqrw 5 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * hqrw 6 c carlos a. felippa, feb. 1967. hqrw 7 c ---------------------------------------------------------------------- c DOUBLE PRECISION / LARGE include 'double.h' c ---------------------------------------------------------------------- c CALLED FROM : dynpr c SUBR. CALLS : stor c ---------------------------------------------------------------------- c hqrw 8 c hqrw 9 dimension g(lf),v(np,m),e(n),a(n),b(n),w(n+1),nd(n),p(n), 1 q(n),xm(n) logical int(n) double precision lambda hqrw 14 if (n.lt.1) go to 1000 precs = 1.0d-15 hqrw 16 base = 2.d0 hqrw 17 ilim = 100 hqrw 18 hov = base**60 hqrw 19 b(1) = 0.d0 hqrw 20 sqrt2 = dsqrt(2.d0) hqrw 21 n1 = n - 1 hqrw 22 do 100 i = 1,n hqrw 23 100 nd(i) = n1*(i-1)-((i-1)*(i-2)/2) hqrw 24 if (n.eq.2) go to 280 hqrw 25 c hqrw 26 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * hqrw 27 c tri-diagonalize matrix r by householder's procedure hqrw 28 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * hqrw 29 c hqrw 30 110 do 250 k = 2,n1 hqrw 31 k1 = k - 1 hqrw 32 kj = k + 1 hqrw 33 ld = nd(k1) hqrw 34 ly = ld + k hqrw 35 y = g(ly) hqrw 36 sum = 0.d0 hqrw 37 do 120 i = kj,n hqrw 38 l = ld + i hqrw 39 120 sum = sum + g(l)**2 hqrw 40 if (sum.eq.0.d0) go to 230 hqrw 41 s = dsqrt(sum+y**2) hqrw 42 b(k) = dsign(s,-y) hqrw 43 w(k) = dsqrt(1.d0+dabs(y)/s) hqrw 46 x = dsign(1.d0/(s*w(k)),y) hqrw 47 do 150 i = k,n hqrw 48 l = ld + i hqrw 49 if (i.gt.k) w(i) = x*g(l) hqrw 50 p(i) = 0.d0 hqrw 51 150 g(l) = w(i) hqrw 52 do 180 i = k,n hqrw 53 y = w(i) hqrw 54 if (y.eq.0.d0) go to 180 hqrw 55 i1 = i + 1 hqrw 56 do 160 j = k,i hqrw 57 l = nd(j) + i hqrw 58 160 p(j) = p(j) + y*g(l) hqrw 59 if (i1.gt.n) go to 180 hqrw 60 do 170 j = i1,n hqrw 61 l = nd(i) + j hqrw 62 170 p(j) = p(j) + y*g(l) hqrw 63 180 continue hqrw 64 190 x = 0.d0 hqrw 65 do 200 j = k,n hqrw 66 200 x = x + w(j)*p(j) hqrw 67 x = 0.5d0*x hqrw 68 do 210 j = k,n hqrw 69 210 p(j) = x*w(j) - p(j) hqrw 70 do 220 j = k,n hqrw 71 lj = nd(j) hqrw 72 do 220 i = j,n hqrw 73 l = lj + i hqrw 74 220 g(l) = g(l) + p(i)*w(j) + p(j)*w(i) hqrw 75 go to 250 hqrw 76 230 g(ly) = sqrt2 hqrw 77 b(k) = -y hqrw 78 do 240 i = kj,n hqrw 79 l = nd(k) + i hqrw 80 240 g(l) = -g(l) hqrw 81 250 continue hqrw 82 280 do 290 i = 1,n hqrw 83 l = nd(i) + i hqrw 84 a(i) = g(l) hqrw 85 290 e(i) = a(i) hqrw 86 b(n) = g(l-1) hqrw 87 c hqrw 88 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * hqrw 89 c get eigenvalues of tridiagonal form by kahan-varah q-r method hqrw 90 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * hqrw 91 c hqrw 92 tol = precs/(10.d0*dble(n)) hqrw 93 bmax = 0.d0 hqrw 94 tmax = 0.d0 hqrw 95 w(n+1) = 0.d0 hqrw 96 do 300 i = 1,n hqrw 97 bmax = dmax1(bmax,dabs(b(i))) hqrw 98 300 tmax = dmax1(bmax,dabs(a(i)),tmax) hqrw 99 scale = 1.d0 hqrw 100 do 310 i = 1,ilim hqrw 101 if (scale*tmax.gt.hov) go to 320 hqrw 102 310 scale = scale*base hqrw 103 320 if (bmax.eq.0.d0) go to 520 hqrw 104 do 330 i = 1,n hqrw 105 e(i) = a(i)*scale hqrw 106 330 w(i) = (b(i)*scale)**2 hqrw 107 delta = tmax*scale*tol hqrw 108 eps = delta**2 hqrw 109 k = n hqrw 110 350 l = k hqrw 111 if (l.le.0) go to 460 hqrw 112 l1 = l - 1 hqrw 113 do 360 i = 1,l hqrw 114 k1 = k hqrw 115 k = k - 1 hqrw 116 if (w(k1).lt.eps) go to 380 hqrw 117 360 continue hqrw 118 380 if (k1.ne.l) go to 400 hqrw 119 w(l) = 0.d0 hqrw 120 go to 350 hqrw 121 400 t = e(l) - e(l1) hqrw 122 x = w(l) hqrw 123 y = 0.5d0*t hqrw 124 s = dsqrt(x) hqrw 125 if (dabs(t).gt.delta) s = (x/y)/(1.d0+dsqrt(1.d0+x/y**2)) hqrw 126 e1 = e(l) + s hqrw 127 e2 = e(l1)- s hqrw 128 if (k1.ne.l1) go to 430 hqrw 129 e(l) = e1 hqrw 130 e(l1) = e2 hqrw 131 w(l1) = 0.d0 hqrw 132 go to 350 hqrw 133 430 lambda = e1 hqrw 134 if (dabs(t).lt.delta.and.dabs(e2).lt.dabs(e1)) lambda = e2 hqrw 135 s = 0.d0 hqrw 136 c = 1.d0 hqrw 137 gg = e(k1)-lambda hqrw 138 go to 450 hqrw 139 440 c = f/t hqrw 140 s = x/t hqrw 141 x = gg hqrw 142 gg = c*(e(k1)-lambda) - s*x hqrw 143 e(k) = (x-gg) + e(k1) hqrw 144 450 if (dabs(gg).lt.delta) gg = gg + dsign(c*delta,gg) hqrw 145 f = gg**2/c hqrw 146 k = k1 hqrw 147 k1 = k + 1 hqrw 148 x = w(k1) hqrw 149 t = x + f hqrw 150 w(k) = s*t hqrw 151 if (k.lt.l) go to 440 hqrw 152 e(k) = gg + lambda hqrw 153 go to 350 hqrw 154 460 do 470 i = 1,n hqrw 155 470 e(i) = e(i)/scale hqrw 156 y = isign(1,m) hqrw 157 do 500 l = 1,n1 hqrw 158 k = n - l hqrw 159 do 500 i = 1,k hqrw 160 if (y*(e(i)-e(i+1)).gt.0.d0) go to 500 hqrw 161 x = e(i) hqrw 162 e(i) = e(i+1) hqrw 163 e(i+1) = x hqrw 164 500 continue hqrw 165 520 if (m.eq.0) go to 1000 hqrw 166 c hqrw 167 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * hqrw 168 c compute eigenvectors by inverse iteration hqrw 169 c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * hqrw 170 c hqrw 171 nvec = iabs(m) hqrw 172 if (nvec.gt.n) nvec = n hqrw 173 f = scale/hov hqrw 174 do 530 i = 1,n hqrw 175 a(i) = a(i)*f hqrw 176 530 b(i) = b(i)*f hqrw 177 sep = 25.d0*tmax*precs hqrw 178 x1 = 0.d0 hqrw 179 x2 = sqrt2 hqrw 180 do 800 nv = 1,nvec hqrw 181 if (nv.gt.1.and.dabs(e(nv)-e(nv-1)).lt.sep) go to 550 hqrw 182 do 540 i = 1,n hqrw 183 540 w(i) = 1.d0 hqrw 184 go to 570 hqrw 185 550 do 560 i = 1,n hqrw 186 x = dmod(x1+x2,2.d0) hqrw 187 x1 = x2 hqrw 188 x2 = x hqrw 189 560 w(i) = x - 1.d0 hqrw 190 570 ev = e(nv)*f hqrw 191 x = a(1) - ev hqrw 192 y = b(2) hqrw 193 j = n1 hqrw 194 do 600 i = 1,n1 hqrw 195 c = a(i+1) - ev hqrw 196 s = b(i+1) hqrw 197 if (dabs(x).ge.dabs(s)) go to 580 hqrw 198 p(i) = s hqrw 199 q(i) = c hqrw 200 int(i) = .true. hqrw 201 z = -x/s hqrw 202 x = y + z*c hqrw 203 if (i.lt.n1) y = z*b(i+2) hqrw 204 go to 600 hqrw 205 580 if (dabs(x).lt.tol) x = tol hqrw 206 p(i) = x hqrw 207 q(i) = y hqrw 208 int(i) = .false. hqrw 209 z = -s/x hqrw 210 x = c + z*y hqrw 211 y = b(i+2) hqrw 212 600 xm(i) = z hqrw 213 if (dabs(x).lt.tol) x = tol hqrw 214 niter = 0 hqrw 215 620 niter = niter + 1 hqrw 216 w(n) = w(n)/x hqrw 217 sum = w(n)**2 hqrw 218 do 640 l = 1,n1 hqrw 219 i = n - l hqrw 220 y = w(i) - q(i)*w(i+1) hqrw 221 if (int(i)) y = y - b(i+2)*w(i+2) hqrw 222 w(i) = y/p(i) hqrw 223 640 sum = sum + w(i)**2 hqrw 224 s = dsqrt(sum) hqrw 225 do 660 i = 1,n hqrw 226 660 w(i) = w(i)/s hqrw 227 if (niter.ge.2) go to 760 hqrw 228 do 700 i = 1,n1 hqrw 229 if (int(i)) go to 680 hqrw 230 w(i+1) = w(i+1) + xm(i)*w(i) hqrw 231 go to 700 hqrw 232 680 y = w(i) hqrw 233 w(i) = w(i+1) hqrw 234 w(i+1) = y + xm(i)*w(i) hqrw 235 700 continue hqrw 236 go to 620 hqrw 237 730 k = j hqrw 238 j = j - 1 hqrw 239 j1 = j - 1 hqrw 240 lj = n1*j1-(j1*(j1-1)/2) hqrw 241 x = 0.d0 hqrw 242 do 740 i = k,n hqrw 243 l = lj + i hqrw 244 740 x = x + g(l)*w(i) hqrw 245 do 750 i = k,n hqrw 246 l = lj + i hqrw 247 750 w(i) = w(i) - x*g(l) hqrw 248 760 if (j.gt.1) go to 730 hqrw 249 if (itape.gt.0) go to 790 hqrw 250 do 780 i = 1,n hqrw 251 780 v(i,nv) = w(i) hqrw 252 go to 800 hqrw 253 790 call dwrite (w,n,itape) hqrw 254 800 continue hqrw 255 do 900 i = 1,n hqrw 256 a(i) = a(i)/f hqrw 257 900 b(i) = b(i)/f hqrw 258 1000 return hqrw 259 end hqrw 260
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The professor is still around and he was from Berkeley from Powell's group,
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
See, that's one of the great things about Fortran. That old code still compiles and (probably) works. That "include" line would not have been from 1967, I'd guess.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Steve:
I have a lot of code from this group - it usually works with minor changes. You are looking at the birth of dynamic analysis in FEM code in the world.
I will take the include -- just not today.
I am using VS 2013 and latest Intel compiler -- is there a part of this system that lets me track the call order and the progress through the program -- this will make debugging a bit easier
John
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
jm-nichols@tamu.edu wrote:
.. I have a lot of code from this group - it usually works with minor changes. .. I will take the include -- just not today.
I am using VS 2013 and latest Intel compiler -- is there a part of this system that lets me track the call order and the progress through the program -- this will make debugging a bit easier ..
See your other thread with some recent comments. The exact same source file you show in the original post here is present in the Drain Building program from year 2010 which appears to have been built successfully using gfortran under mingw in that year if one can trust the date on the build scripts. So what you are showing here is most likely a 2010 version of code rather than from 1967 and it's not that big a deal - there are umpteen subroutines out there in each and every field involving computations that look like this.
And why do you think you need to make any changes at all to this subroutine, even "minor"? If you can get your hands on the cases that were run to validate the 2010 gfortran build, why wouldn't you get the same answers? If there are any issues, why wouldn't they be algorithmic rather than Fortran-related?
"is there a part of this system that lets me track the call order and the progress through the program -- this will make debugging a bit easier .." - you should be able to run any of the profiling tools e.g., gprof on *nix system or Intel VTune Analyzer on Windows.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Line 234 may present an issue:
if (int(i)) go to 680
That statement is formulated as a LOGICAL, yet is using an INTEGER. Even if the compiler would re-interpret "cast" the integer to logical, you may have an issue regarding implementation issues of what is the bit pattern for .TRUE. vs .FALSE..
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Yes the init had thrown a flag when I compiled, I changed it to a integer to match the call - need to change it back.
Thanks
John
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
If not run from the debugger, Re-run the program from the debugger.
When the error occurs, click on the Retry. The debugger should trap, but it will likely trap inside (several calls deep) the Fortran and C Runtime library. Look at the Call Stack window (Debug | Windows | Call Stack). As you read lines down in the Call Stack you may be able to locate the first level inside a Fortran source. If you double click on that, it should open that file, and then you can look around for cause.
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello,
I want to apply IMSL library in PGI visual fortran 2013, but my codes were been written in f90,
after starting, the "This project is out of date" comes up, and also ask "Would you like to built it?".
then after continuing, again this error message is shown,"Visual studio cannot debug because the debug target 'C:\Users\wave\documents\visual studio 2013\projects\PVFProject5\PVFProject5\Win32\Debug\PVFProject5.exe'is missing".
USE MSIMSLMD
EXTERNAL FCN
Double precision v1,v2,v3,k12,k23,k13,MV1,MV2,MV3
Double precision x,xguess,fnorm,errrel,v
Dimension x(1),xguess(1),fnorm,errrel,v
Common/const/k12,k23,k13,MV1,MV2,MV3,v3
Data k12,k23,k13,MV1,MV2,MV3,v3
Data k12,k23,k13,MV1,MV2,MV3/0.9,0.5,1.11,18.0,127.8,30532.0/
ERRREL=1E-10
ITMAX=300
N=1
c
open(unit=1,file='spin1.txt')
2 format(20x,16HVolume Fractions,//,12x,2Hv1,25x,2Hv2,22x,2Hv3)
write(1,2)
V(1)=0.4
C Initial guess
10 write(*,*)'enter v3'
read(*,*)v3
xguess(1)=v(1)
CALL DNEQNF(FCN,ERRREL,N,ITMAX,XGUESS,X,FNORM)
v(1)=x(1)
v2=1-v(1)-v3
write(1,*)v(1),v2,V3
GOTO 10
clse(unit=1)
end
C IMSL required subroutine FCN
Subroutine fcn(x,f,n)
Double precision v1,v2,v3,k12,k23,k13,MV1,MV2,MV3
Double precision f,x,G22,G23,G33
Dimension f(n),x(n)
common/const/k12,k23,k13,MV1,MV2,MV3,v3
v1=x(1)
v2=1-v1-v3
G22=1/V1+MV1/(MV2*V2)-2*K12
G23=1/V1-(K12+K13)+MV1/MV2*K23
G33=1V1+MV1/(MV3*V3)-2*K13
f(1)=G22*G33-G23**2
return
end
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Your compilation likely has an error. Thus no executable is produced, and thus not available to be debugged.
Your code contains "clse(unit=1)". If this is in your program (as opposed to typographical error in posting) then your program will contain a compilation error for this statement.
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
PGI has its own forum at http://www.pgroup.com/userforum/viewforum.php?f=4 . Since you are using PGI fortran, you should ask your questions in that forum.
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page