UMACR4.F - - Broyden Quasi-Newton Example

The umacr[0-9].f files in the user directory allow one to create user macro commands. To set up a user macro command one needs to edit the file umacr4.f (for example) to include the code one wishes to have executed when the user mesh command is called from inside FEAP. The interface can be programed in fortran or it can be programed in C. In the example that follows a fortran interface is utilzed. The example shows how one can interface an unsymmetric Broyden quasi-Newton method with FEAP.

The main steps in creating a user macro command whether in C or Fortran are

Set the name of the macro command to whatever 4-letter name you wish it to have when FEAP calls it.
Set up the source code you wish to have executed when the command is called from FEAP

The macro command name is set in the common variable uct in the common block umac1. When uct is equal to 'mac4' (for umacr4.f or .c) it should be set to whatever you wish for the name of your mesh command.

       A example umacr4.f is given below. The name of the user macro is set to broyden. The macro command itself is used to perform a given number of quasi-newton updates to the solution vector without form-factoring an additional tangent matrix. The command takes as options the keyword iter and the number of iteration to be performed. The sequence of calls from the Macro prompt would then be:
    Macro> utan,,1
    Macro> broy,iter,10
if one were interested in ten broyden iterations.

As these routines use the user memory allocation system. The needed ualloc file is also contained below.

Source file:

c$Id: umacr4.f,v 1.1 1999/03/21 07:53:57 sanjay Exp $
subroutine umacr4(lct,ctl,prt)

c * * F E A P * * A Finite Element Analysis Program

c-----[--.----+----.----+----.-----------------------------------------]
c      Purpose: User interface for the Broyden Quasi-Newton Method
c
c      Call order: > utan,,1
c                  > broy,iter,15

c      Inputs:
c         lct       - Command character parameters
c         ctl(3)    - Command numerical parameters
c         prt       - Flag, output if true

c      Outputs:
c         N.B. Users are responsible for command actions. See
c               programmers manual for example.
c-----[--.----+----.----+----.-----------------------------------------]

implicit none

      include 'cdata.h'
      include 'iofile.h'
      include 'pointer.h'
      include 'umac1.h'
      include 'upointer.h'
      include 'comblk.h'

      logical   pcomp,prt,tlog,ualloc
      character lct*15
      integer   numvec,vecs,i
      real*8    ctl(3)

save

c Set command word

      if(pcomp(uct,'mac4',4)) then
        uct = 'broy'
        numvec = 0
        return

elseif(pcomp(lct,'iter',4)) then

vecs = min(max(int(ctl(1)),2),20) ! allow for a maximum of 20 vectors and min of one iteration (2 vecs)

        if(numvec .lt. vecs) then ! Re-allocate to get more memory if necessary
          numvec = vecs
          tlog = ualloc(1,'VTIL',neq*numvec,2)
          if(.not.tlog) write(*,*) 'UALLOC ERROR - VTIL'
          tlog = ualloc(2,'DELX',neq*numvec,2)
          if(.not.tlog) write(*,*) 'UALLOC ERROR - DELX'
        endif

call pzero(hr(up(1)),neq*numvec)
call pzero(hr(up(2)),neq*numvec)

c       Grab delta-x from previous utang,,1 solve
        do i = 1,neq
         hr(up(2) + i - 1) = hr(np(26) + i - 1)
        end do

call broyden(hr(up(1)),hr(up(2)),vecs)

endif

end

subroutine broyden(vtil,delx,numvec)

implicit none

      include 'cdata.h'
      include 'fdata.h'
      include 'hdatam.h'
      include 'ndata.h'
      include 'part0.h'
      include 'pointer.h'
      include 'print.h'
      include 'prlod.h'
      include 'comblk.h'

      logical fa,tr
      integer numvec,iter,k,pu,pr,fp(5),j
      real*8 w,z,energy,rnorm,dot
      real*8 vtil(neq,numvec),delx(neq,numvec)
      real*8 v(neq)

save

      pu = np(40)
      pr = np(26)
      fa = .false.
      tr = .true.

c Iterate for the updates
do iter = 2,numvec

c      Allow for history updates
       hflgu = tr
       h3flgu = tr

c     Get external load vector
       call pload(mr(np(31)),hr(pu),hr(np(30)),hr(pr),prop,tr,fa)
c      Get total negative residual in hr(pr)
       call formfe(pu,pr,pr,pr,fa,tr,fa,fa,6,1,numel,1)

rnorm = sqrt(dot(hr(pr),hr(pr),neq))
write(*,*) 'Iter = ',iter,' Residual = ',rnorm

c      Set flags for u-symm solve and solve with result found in hr(pr)
       fp(1) = na
       fp(2) = nau
       fp(3) = nal
       fp(4) = pr
       fp(5) = np(20+npart)
       call psolve(fp,fa,tr,tr,prnt)

       do k=1,neq
        v(k) = -hr(pr+k-1)
       end do

if(iter.eq.2) then

c       Broyden update computation first one
         w = dot(v,delx(1,iter-1),neq)
         z = dot(delx(1,iter-1),delx(1,iter-1),neq)
         z = -1.d0/(w+z)
         do k = 1,neq
           vtil(k,iter-1) = v(k)*z
           delx(k,iter)   = -(v(k) + w*vtil(k,iter-1))
         end do

else

c        Broyden update computation generic form iter > 2
         do j = 1,iter-2
          w = dot(v,delx(1,j),neq)
          do k = 1,neq
            v(k) = v(k) + vtil(k,j)*w
          end do
         end do

         w = dot(v,delx(1,iter-1),neq)
         z = dot(delx(1,iter-1),delx(1,iter-1),neq)
         z = -1.d0/(w+z)
         do k = 1,neq
           vtil(k,iter-1) = v(k)*z
           delx(k,iter)   = -(v(k) + w*vtil(k,iter-1))
         end do

endif

c      Actual update of solution vector
       call update(mr(np(31)),hr(np(30)),hr(pu),hr(np(42)),
     &             delx(1,iter),fl(9),2)

end do

end

c User allocation routine that goes with broyden macro
c$Id: ualloc.f,v 1.1 1999/03/21 07:53:57 sanjay Exp $
logical function ualloc(num,name,length,precis)

c * * F E A P * * A Finite Element Analysis Program

c-----[--.----+----.----+----.-----------------------------------------]
c      Purpose: Define, delete, or resize a user dictionary entry.
c               Pointer defined for integer (single) and real
c               (double precision) arrays.

c N.B. Currently limited to 200 names by dimension of
c common blocks 'allotd','allotn','pointer','upointer'

c Inputs:

c         num        - Entry number for array (see below)
c         name       - Name of array          (see below)
c         length     - Length of array defined: =0 for delete
c         precis     - Precision of array: 1 = integers; 2 = reals

c Output:

c up(num) - Pointer to first word of array in blank common
c-----[--.----+----.----+----.-----------------------------------------]

implicit none

include 'allotd.h'

      logical   usetmem
      character name*(*)
      integer   i, num,length,precis

c Storage definitions for UALLOC variables

integer list
parameter (list = 2)

character names(list)*5

save

c Define and store names

data (names(i),i=1,list)/

& 'VTIL','DELX'/

c Short description of variables

c 'VTIL', ! 1: v-tilde array in broyden method
c 'DELX', ! 2: array of old incremental updates

c Do memory management operations

ualloc = usetmem(list,names,num,name,length,precis)

end