Fail to understand code

[nichols__john]

c ---   initialize
        call izero (l(1),kndid,ksofar)
100     Continue
        call indco2 (nnods,l(kndid),l(kcoord),kdata)

L is the blank common.

izero attempts to set l from kndid to ksofar to all equal zero. kndid is 2 and ksofar is 17

But when izero returns l is only one long at the continue

c **********************************************************************
      SUBROUTINE IZERO (ia,n1,n2)
c **********************************************************************
c     PURPOSE : Zero terms n1 through n2 of an iteger array
c ----------------------------------------------------------------------
c     DOUBLE PRECISION / LARGE
      include 'C:\Users\macne\Documents\Visual Studio 2013\Projects\Pro
     1gram067 - Drain3DX\Ares\Ares\Ares\common\double.h'
c ----------------------------------------------------------------------
c     ARGUMENTS
c     ia = array
c     n1,n2 = range to be zeroed
c ----------------------------------------------------------------------
c     ARGUMENT DECLARATIONS
      dimension ia(n2)
c **********************************************************************
      do 10 n=n1,n2
        ia(n) = 0
   10 continue
c **********************************************************************
      RETURN
      END

Any ideas would be appreciated as this is weird to me

Thanks

John

[mecej4]

It is usually incorrect to assume that excerpts from a big piece of code will suffice to reproduce errors that are being seen in the big code. In this case, it may well be the case that at the place where IZERO is called there is no information concerning the extent of the array L().

But when izero returns l is only one long at the continue

How did you reach this conclusion? Using added WRITE statements? Inspecting the Locals tab of the VS debugger?

[nichols__john]

Dear Mecej4

L is declared as a vector of unit length in the blank common files. I slowly added the subroutines to the program to ensure that it at least compiled in VS 2013. The program uses a lot of common blocks to hold the data. This is similar to the original ULARC and AXISHELL programs from the same research group. In the 1980's I rewrote ULARC  to eliminate common blocks.  It is however 2D and does not do dynamic analysis.

I have a manual, but the code does not match the stated input format so I have to work through the code to make sure the input file is in the correct format. I am up to to the nodes and the node co-ordinates. I used VISUAL Studio Debugger and Watch Locals to look at L when it is being passed to the subroutine.  It has a unit length and the error checker says -- the subroutine is using a longer array to hold the node numbers, which should be stored in L. I am not a UNIX boffin, but it would appear that the UNIX Fortran allows for the fact that in the main block the L is declared as a million units, but this data is not available to the subroutine, so Intel trips an execution error. It sees the unit declaration in the blank common file.

I than changed to the blank common file to match the declaration in the main unit and the problem is not currently evident. The input data file format is not directly evident as the code uses scratch files to store copies of the data and then the program rereads the data lines using different routines. So the input is spread out even for a single data line.  This code also has a new node type being added, but there is limited data about it. Removing the code caused other errors as some of the node constants were declared in the new routine.

I can use a commercial FEM package, but UCB code allows me to play directly with the stiffness matrix and in terms of the research this is a big plus.

In the last week I have been thinking it will be better to fix the code and remove the commons and simplify the input. It is really simple input.  I was hoping that would not be necessary.

So at the moment I am debating the alternatives. Benedetti wants to use Strand7, we both have it - but I am not sure that FEM package will handle this analysis as well as required. One could use ABAQUS or other programs, but for research work they are painful. I am going to spend several months on it so giving up a month to get DRAIN3DX without commons is an attractive idea.

As usual I would appreciate your insights. One could change the solver to PARSISO and introduce Newmark techniques to solve the dynamic equation. Or use UCB Unix.

John

[mecej4]

Codes of that vintage, especially if targeted to mainframes, used static allocation, and local variables were saved. Dummy 1-D array arguments would be declared with a dummy (1) dimension, which in modern terms would be declared with (*). Similarly for higher rank array arguments.

What this implies is that you cannot compile much, if any, of such code with subscript checks turned on. Note that Silverfrost Fortran 95 provides the option /Old_arrays, which skips subscript checking for such variables while doing checking for other arrays.

Intel Fortran is well equipped to handle such old code, but one must understand the genre, and refrain from applying modern checks. I do not think that it is worth the effort to fix the old code just to make it conform to modern expectations.

Had the source code been available, I would have been willing to try to get it to run, but I note that there is a fee to be paid before the source code can be obtained. Perhaps, there are other versions available, with no fee to be paid?

[nichols__john]

Dear Mecej4

This is the original code. 

The input file is how this code is set up -- it is slightly different to the manual.

I would prefer as you say to set the compiler options to run the code. It is good code and not hard to follow - just not what I am used to.

I have never used SAVE and L(*) is not something I normally use.

I will try them later once I hang 3000 milk bottles from a rood\f -- do not  ask

John

[nichols__john]

Dear Mecej4:

I added save and it appears to have run through.

The input file is flawed - I can work to fix that -- there is some interesting code that uses a lot od scratch files -- they must have had no memory.

Thanks

In the blank common is it better to use (*) or the main codes 1000000 as declared.

I turned off the interface checker as suggested earlier.

The basic premise is the dimensions on a timber beam vary with heat (very small amounts) I varies and we get a Brownian motion, along with the Brownian motion from the restraint nodes. (Here I am using BM broadly it is really a thermal movement) A normal FEM would not cope with the idea.

John

[LRaim]

You should make sure to compile without interface checking and without checking the boundary of arrays.

However the code is clear.

In this case the definition of IA as "DIMENSION IA(N2)" in the IZERO the subroutine tells that IA is a vector.

Perfectly equivalent definitions would be: 
DIMENSION IA(1)
DIMENSION IA(N1)
DIMENSION IA(*)
DIMENSION IA(any positive number)

In the call statements the parameter L(1) passes the base address of the blank common and the IZERO routine will clear the storage from kndid to ksofar.

[jimdempseyatthecove]

I assume you may have learned how to program in C. Consider:

struct common_t
{
  int L[1];
  double a,b,c;
  float whatsAmaJig[12345];
} common;
...
by indexing common.L[index] with index .gt. 0 you can reach into locations of the struct that do not belong to L. This assumes you have no runtime checking for array subscripting error checks enabled.

It was a common practice when writing out (to file) or reading back (from file) to use unformatted binary WRITE/READ statements using the array l, miss-sized as a proxy for the common block (or portion of your choosing). This increased the performance of I/O and relieved the programmer from updating the I/O routine whenever the struct layout changed (possibly the length of the I/O would change).

Using named commons is similar but better in light that for the same named common, different program units can effectively map the named common in different way. Essentially this is equivalent to a UNION. Using the same named common, the I/O routine would use the declaration that contains the integer :: L(yourMaxSize), and the routines manipulating the data would use the declaration that contained everything else. This eliminates having the extraneous array L in your data. Note, the L array of one element was often repurposed to hold a record ID or record size, record header, or such that the programmer deemed necessary for the I/O routine.

Jim Dempsey

[mecej4]

John, the .TXT files in the Zip indicate that the package is version 1 of Drain3d, and the set-up is for Microsoft FPS 1 with a DOS Extender. Do you have a Unix-configured version of Drain3d, or do you have version 2 of Drain3d? It would be a lot easier to get a Unix-configured version running on Windows than a 16-bit DOS configured version.

This link mentions a 2010 edition: http://nisee.berkeley.edu/elibrary/getpkg?id=DRAIN3DX .

[andrew_4619]

l(1) is passes the address of the start of l. The next 17 locations in memory (of 4 bytes each assuming 4 bye integers) get set to zero. The could be entirely different variables to l, there are whatever is next in the common block....

[nichols__john]

I will get the latest program tomorrow when the office is open.

John

[nichols__john]

Dear mecej4:

I obtained the 2010 version that compiles on GFORTRAN - batch files included in zip.

I have not used GFORTRAN at all - there are no sample files I can find.

JOhn

[mecej4]

There is not much that we can achieve by building an old program with no user manual and no sample input data.  Nevertheless, it seems to me that there is no need for saving local variables. There are only two source files in the sources/main directory, and from them I see that the bulk of the problem data is stored in blank common, so those variable values always remain in scope. The source code can be compiled using almost any F77 compiler, provided no checking of subscripts, call consistency, etc. is attempted. For my tests, I reduced the size of the blank common to a value sufficient for the sample problem input file.

P.S.: There are some input files listed from p.207 onwards in this report: http://mae.cee.illinois.edu/publications/reports/Report00-05.pdf .

I found the sister program DRAIN2DX at ftp://ftp.ecn.purdue.edu/ayhan/SBE/DRAIN2DX . I found that the IFort compiled EXE and the EXE provided in that distribution give nearly identical outputs. The generated file la3ns.ech contains, in either case:

 ERROR - opening file, iv79elcn.acn    in INAXL.
 ERROR: Record  SAC not found

[John_Campbell]

The file "iv79elcn.acn" is probably a data file for an Elcentro earthquake ground acceleration history. There must be another .zip file of data files and sample input files for this program. I recall a reference earthquake from 1940, so not sure of 79. It may be a smaller earthquake that had a better time history data set and produces better answers.

[nichols__john]

I have been playing with the original program - not the ones I got today.

The problem appears to be that someone commented out all of the file numbers, I have slowly been putting them back in.

One of the manuals notes that there are no sample files, so I am slowly creating one - not that hard as all structures are the same. it is just tedious to check the manual against the code and create/ fix the data file.

The version I have started with includes some early code for compound nodes, but it is not complete. They have used a 1000 lines of code to enter simple data, really a 100 line input program would have sufficed.

One interesting problem that the English Professor has created is needing a moving FFT -- I am thinking of stepping one time step ahead and then repeating --  2000 FFT's per second - anyone seen this done before?  Numerically intensive.

Clearly the Drain coders never read Winston and Horne or Abelson - using a if, goto and exit in one block makes for a challenge to debug.

John

[nichols__john]

I am finally starting to understand the ideas behind the code.  The manuals are poorly documented, but the code is not so it is a long slow but interesting process.

The PhD thesis for the original writer or DRAIN3DX notes that it did not play well with Windows Fortran -- I can see why, but it is not hard now to fix.  One could never use this commercially - but for research with the English Professor it has a number of advantages.

The matrix inversion routine is primitive and I was thinking of making it PARDISO to allow for a decent problem size.  The flexibility to play with the code really out weighs the minor annoyances to get it running.

c **********************************************************************
      SUBROUTINE OPTSOL (a,b,na,neq,jcol,kex,leq,lenk)
c **********************************************************************
C     Equation solver.
c ----------------------------------------------------------------------
c     CALLED FROM : flex,grsol,mode,rest,static,step,update,updats
c ----------------------------------------------------------------------
c     DOUBLE PRECISION / LARGE
      include 'C:\Users\macne\Documents\Visual Studio 2013\Projects\Pro
     1gram067 - Drain3DX\Ares\Ares\Ares\common\double.h'
c ----------------------------------------------------------------------
c     ARGUMENTS
c       INPUT
c        a(lenk) = column compacted matrix.
c        na(neq) = storage location of the diagonal terms in A.
c        neq = no. of equations.
c        jcol = no. of first changed column in matrix A.
c        kex = 1: reduce coeff. matrix only.
c              2: reduce vector b and back substitute.
c              3: 1+2.
c              4: reduce vector b only.
c              5: back substitute only.
c      MODIFY
c        b(neq) = solution vector
c ----------------------------------------------------------------------
c     LABELLED COMMONS
      include 'C:\Users\macne\Documents\Visual Studio 2013\Projects\Pro
     1gram067 - Drain3DX\Ares\Ares\Ares\common\tapes.h'
c ----------------------------------------------------------------------
c     ARGUMENT DECLARATIONS
      dimension a(lenk),b(neq),na(neq)
c **********************************************************************
c -------------------------------------------------------SINGLE EQUATION
      if(neq.eq.1)then
        if(kex.eq.1 .or. kex.eq.3) then
          a(1)=1./a(1)
          if(kex.eq.3) b(1)=b(1)*a(1)
        else if(kex.eq.2 .or. kex.eq.5) then
          b(1)=b(1)*a(1)
        endif
        RETURN
      endif

c --------------------------------------------------------------REDUCE A

      neqq=neq-1

      if(kex.eq.1 .or. kex.eq.3)then
   10   jf=max0(jcol,2)
        j1=jf+1
        il=jf-1
        jcl=jcol-1
        najp=na(il)
        do 140 j=jf,neq
        naj=na(j)
        if=j1-naj+najp
        if (if.ge.j) go to 130
        if1=max0(if+1,jcol)
        jk=naj-j
        if (if1.gt.il) go to 80
        jia=jk+if1
        i1=if1+1
        kl=if1-1
        naip=na(kl)
        do 70 i=if1,il
        nai=na(i)
        ik=nai-i
        ii=i1-nai+naip
        if (ii.ge.i) go to 60
        kf=max0(ii,if)
        krl=min0(kl,leq)
        if(kf.gt.krl)go to 60
        jka=jk+kf
        ika=ik+kf
        aa=a(jia)
        if (kf.ge.jcol) go to 30
        jrl=min0(krl,jcl)
        do 20 k=kf,jrl
        nak=na(k)
        aa=aa-a(jka)*a(ika)*a(nak)
        jka=jka+1
   20   ika=ika+1
        if (jcol.gt.krl) go to 50
        kf=jcol

   30   do 40 k=kf,krl
        aa=aa-a(jka)*a(ika)
        jka=jka+1
   40   ika=ika+1
   50   a(jia)=aa
   60   jia=jia+1
        i1=i1+1
        kl=kl+1
   70   naip=nai

   80   kf=if
        jka=jk+if
        aa=a(naj)
        krl=min0(il,leq)
        if(kf.gt.krl)go to 130
        if (if.ge.jcol) go to 100
        jrl=min0(krl,jcl)
        do 90 k=if,jrl
        nai=na(k)
        aa=aa-a(nai)*a(jka)**2
   90   jka=jka+1
        if (jcol.gt.krl)go to 120
        kf=jcol

  100   do 110 k=kf,krl
        nai=na(k)
        cc=a(jka)/a(nai)
        aa=aa-a(jka)*cc
        a(jka)=cc
  110   jka=jka+1
  120   a(naj)=aa
  130   il=il+1
        j1=j1+1
  140   najp=naj
      endif

c -------------------------------------------REDUCE B AND BACKSUBSTITUTE

      if(kex.eq.2 .or. kex.eq.3 .or. kex .eq.4)then
        do 160 n=1,neqq
        if (b(n).ne.0.) go to 170
  160   continue
        n=neqq
  170   n1=n+1
        i1=n1+1
        kl=n
        naip=na(n)
        do 200 i=n1,neq
        nai=na(i)
        ii=i1-nai+naip
        if (ii.ge.i) go to 190
        kf=max0(ii,n)
        krl=min0(leq,kl)
        if(kf.gt.krl)go to 190
        ik=nai-i
        ika=ik+kf
        bb=b(i)
        do 180 k=kf,krl
        bb=bb-a(ika)*b(k)
  180   ika=ika+1
        b(i)=bb
  190   i1=i1+1
        kl=kl+1
  200   naip=nai
        if(n.le.leq)then
          do 210 i=n,leq
          nai=na(i)
  210     b(i)=b(i)/a(nai)
        endif
      endif

c -------------------------------------------------------BACK SUBSTITUTE

      if(kex.eq.2 .or. kex.eq.3 .or. kex.eq.5)then
        j=neq
        if(j.ge.2)then
          j1=j+1
          naj=na(j)
          do 240 i=1,neqq
          najp=na(j-1)
          ii=j1-naj+najp
          if (ii.ge.j) go to 230
          jk=naj-j
          kf=ii
          kl=min0(j-1,leq)
          if( kf.gt.kl)go to 230
          jka=jk+kf
          bb=b(j)
          do 220 k=kf,kl
          b(k)=b(k)-a(jka)*bb
  220     jka=jka+1
  230     j1=j1-1
          j=j-1
  240     naj=najp
        endif
      endif
c **********************************************************************
      RETURN
      END

[mecej4]

Unless you know that OPTSOL is the bottleneck, it seems premature to talk of replacing calls to it with Pardiso calls. So far, I have not seen a single input file that can be run with either of the Drain3D source packages that you posted, so we have zero information on run timing. Secondly, Optsol calls do not seem plug-compatible with Pardiso calls. I do not know what type of sparse matrix storage is expected by Optsol.

Given this state of affairs, I thought that it might be useful to ask the same kind of questions regarding Drain2D. Using the source code in ftp://ftp.ecn.purdue.edu/ayhan/SBE/DRAIN2DX/Drain2DX-2010.zip, I built the program. I used the input file in the same zip file, and found the missing acceleration data in ftp://ftp.ecn.purdue.edu/ayhan/SBE/DRAIN2DX/DRAIN-2DX-UIUC.zip, thanks to John Campbell's hint in #14. The program, built using IFORT 2016 Update 1, ran with no errors. With /Od /Zi, the run spent 34 percent of the time in Optsol. The source code of Optsol in Drain2D is identical to the one in #16. Even if you could replace Optsol with a solver that took zero time, the overall run time would only be improved by a third. If we guess that these findings will carry over to Drain3D, I think that you should postpone trying to do anything to "improve" the source code: It ain't broke yet!

[nichols__john]

Dear mecej4:

I enclose an input file that halts at nodal loads - I am up to that point.

*START 
  TIMBER           0 0 0 0   f          Session One
*NODECOORDS                                                                                   
C        1       0.0       0.0       0.0
C        5      10.0      10.0      10.0
L        1         5         1   0
*NODETYPES
S        1         1         5 
*RESTRAINTS
S   111000 1
S   111000 5  
*MASSES
S 111       1.0         1         5                              9.81       0.0  
*ELEMENTGROUP
   1    0    0            0.1           
   1 
   1 210000000       0.5     0.011      200.0    200.0    1       1.0 
   1         1         4         1    1
*PARAMETERS
OS
       1   -1   -1    0    1
*GRAV                                     STATIC ANALYSIS
N     PAT1      1.0
N     PAT2      1.0
*REST

 I would love to use DRAIN2 - but the data coming from the English Professor is 3D and I want to match her data. DRAIN 2 has been used in a lot of places, so it has been sanity checked. I am not sure, but it appears that DRAIN3d has had limited use. This is not unusual with SEMM programs, as far as I am aware I am the only person in the world using ULARC and AXISHELL. Not that that says much.

I enclose the manuals I have in PDF form.  It appears that the author tried as many Fortran versions as possible and noted problems with some, including Windows

Interestingly I have been stopped for a while whilst I sorted out the truss lengths. The coord array was placed in L, but it did not come out as required. Rather than waste a lot of time trying to understand stuff I do not understand - I merely put the coord array in the necessary calls and the error code went away.

As usual your observations are good and I will leave optsol alone. I am sure if I look and I did - optsol is in ULARC from 1972 but not from AXISHELL from 1968. Here is the 1972 version with variables declared now.  The two are not quite the same- but I could make a good case that ULARC was used to create the latest optsol.

it is not broken - it has worked for more than 40 years so I will leave it alone.

John

SUBROUTINE OPTSOL (A,B,NA,NEQ,LIMS,LIMF,KKK)
      
      Implicit None
C
      REAL A(1),B(1)
      INTEGER NA(1)
      INTEGER NEQQ,J1,IL,NAJP,J,NAJ,IF,IF1,JK,JIA,NAIP,I1,KL,I
      INTEGER NAI,IK,II,KF,JKA,IKA,K
      Integer NEQ, KK,LIMS, LIMF, MAX0, KKK
      REAL AA
      REAL CC
      INTEGER N,N1
      REAL BB
C
      NEQQ=NEQ-1
      GO TO (10,70), KKK
C     TRIANGULARIZE MATRIX A
   10 J1=3
      IL=1
      NAJP=NA(1)
      DO 60 J=2,NEQ
      NAJ=NA(J)
      IF=J1-NAJ+NAJP
      IF1=IF+1
      JK=NAJ-J
      IF (IF1.GT.IL) GO TO 40
      JIA=2+NAJP
      NAIP=NA(IF)
      I1=IF1+1
      KL=IF
      DO 30 I=IF1,IL
      NAI=NA(I)
      IK=NAI-I
      II=I1-NAI+NAIP
      KF=MAX0(II,IF)
      JKA=JK+KF
      IKA=IK+KF
      AA=A(JIA)
      DO 20 K=KF,KL
      AA=AA-A(JKA)*A(IKA)
      JKA=JKA+1
   20 IKA=IKA+1
      A(JIA)=AA
      JIA=JIA+1
      I1=I1+1
      KL=KL+1
   30 NAIP=NAI
C
   40 JKA=JK+IF
      AA=A(NAJ)
      DO 50 K=IF,IL
      NAI=NA(K)
      CC=A(JKA)/A(NAI)
      AA=AA-A(JKA)*CC
      A(JKA)=CC
   50 JKA=JKA+1
      A(NAJ)=AA
      IL=IL+1
      J1=J1+1
   60 NAJP=NAJ
C
      GO TO 150
C     REDUCE VECTOR B AND BACK SUBSTITUTE
   70 DO 80 N=LIMS,NEQQ
      IF (B(N).NE.0.) GO TO 90
   80 CONTINUE
      N=NEQQ
   90 N1=N+1
      I1=N1+1
      KL=N
      NAIP=NA(N)
      DO 110 I=N1,NEQ
      NAI=NA(I)
      II=I1-NAI+NAIP
      KF=MAX0(II,N)
      IK=NAI-I
      IKA=IK+KF
      BB=B(I)
      DO 100 K=KF,KL
      BB=BB-A(IKA)*B(K)
  100 IKA=IKA+1
      B(I)=BB
      I1=I1+1
      KL=KL+1
  110 NAIP=NAI
      DO 120 I=N,NEQ
      NAI=NA(I)
  120 B(I)=B(I)/A(NAI)
C
      J=NEQ
      J1=J+1
      KL=NEQQ
      NAJ=NA(NEQ)
      DO 140 I=LIMF,NEQQ
      NAJP=NA(J-1)
      II=J1-NAJ+NAJP
      JK=NAJ-J
      KF=MAX0(II,LIMF)
      JKA=JK+KF
      BB=B(J)
      DO 130 K=KF,KL
      B(K)=B(K)-A(JKA)*BB
  130 JKA=JKA+1
      J1=J1-1
      KL=KL-1
      J=J-1
  140 NAJ=NAJP
C
  150 RETURN
      END

[nichols__john]

dear mecej4

I got the program to run with a simple input file, the results are wrong -- but it is at least working.  I can now create a real job and sort out the math.

Good fun.

I found a PhD thesis from the UIUC that used DRAIN3DX and he got similar results to ABAQUS - so that is good.

John

[John_Campbell]

It is not a good idea to be changing the program to insert a different equation solver, before you have the original program working correctly. While OPTSOL will consume a significant proportion of run time, these are not typically large equation models. It may be written in old form Fortran, but it still works well.

Identify the cause of the conversion problem and get the program working first. This is not due to a slow equation solver or program bugs but how you are managing the code. Didn't it come with a make file ?

You should find that for the sample problems of 20 years ago, they will run very quickly.

Convergence speed is more affected by other parameters, which you can review with a working program.

I am not sure of your earlier comments of using FFT, I presume for adaptive iterations; it could have some merit, but I would hope you have a sound operating basis with Drain3D first. You might be onto an interesting idea or not.  What advice has your English professor provided ?

[nichols__john]

your points are well made.

I will concentrate on getting the code running.

John