Community
cancel
Showing results for
Did you mean:
Beginner
34 Views

## Old Code

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

```

Tags (1)
11 Replies
Beginner
34 Views

The professor is still around and he was from Berkeley from Powell's group,

Employee
34 Views

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.

Beginner
34 Views

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

Honored Contributor I
34 Views

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.

Black Belt
34 Views

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

Beginner
34 Views

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

Beginner
34 Views

How do I find out what causes this error

Thanks

John

Black Belt
34 Views

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

Beginner
34 Views

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'
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

Black Belt
34 Views

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

Black Belt
34 Views

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.