Community
cancel
Showing results for 
Search instead for 
Did you mean: 
nichols__john
Beginner
86 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

 

0 Kudos
11 Replies
nichols__john
Beginner
86 Views

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

 

Steven_L_Intel1
Employee
86 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.

nichols__john
Beginner
86 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

FortranFan
Honored Contributor I
86 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.

jimdempseyatthecove
Black Belt
86 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

nichols__john
Beginner
86 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

nichols__john
Beginner
86 Views

How do I find out what causes this error

Thanks

John

jimdempseyatthecove
Black Belt
86 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

farzin_s_
Beginner
86 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'
          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

jimdempseyatthecove
Black Belt
86 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

mecej4
Black Belt
86 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.