subroutine rowmbs(nx,ng,t0,te,h0,x,v,
     &  rtol,atol,itol,
     &  ffunc,mfunc,gxfunc,gfunc,
     &   work,lwork,iwork,liwork,idid)
C*************************************************************
C
C  ROWMBS4 implements a 4th order rosenbrock method
C  for the solution of mechanical multibody systems   
C     
C      x'=v               T
C  m(x)v'=f(t,x.v)-(dg/dx)  lambda
C      0 =g(x)
C      
C  in index 3 formulation.
C  (Note: only the applied forces may depend on t)     
C
C  AUTHOR: J. Wensch, MLU HALLE, FB MATHEMATIK UND INFORMATIK, 
C          INSTITUT FUER NUMERISCHE MATHEMATIK
C          wensch@mail.mathematik.uni-halle.de
C
C  version of oct 2 1997
C
C*************************************************************  
C
C  input parameters:
C  
C  nx: dimension of position coordinates x
C
C  ng: dimension of constraints g(x)
C
C  t0: left endpoint of integration interval
C
C  te: right endpoint of integration interval    
C      
C  h0: initial stepsize, if h0=0.d0 the code puts
C      h0=(te-e0)*rtol**0.5d0
C      
C  x,v: initial values in position and velocitiy
C      
C  rtol: relative tolerance    
C
C  atol: absolute tolerance    
C      
C  ffunc:  Name (external) of user-supplied subroutine computing the 
c          value f(t,x,v): 
c                subroutine ffunc(n,t,x,v,f)
c                real*8 t,x(n),v(n),f(n)
c                integer n
c                f(1)= ...
c            
C  mfunc:  Name (external) of user-supplied subroutine computing the 
c          value m(x): 
c                subroutine mfunc(n,ld,t,x,amat)
c                real*8 t,x(n),amat(ld,n)
c                integer n,ldamat
c                amat(1,1)= ...
c            
C  gxfunc: Name (external) of user-supplied subroutine computing the 
c          value dg/dx: 
c                subroutine gxfunc(m,n,ld,t,x,gx)
c                real*8 t,x(n),gx(ld,n)
c                integer n,m,ld
c                gx(1,1)= ...
c            
C  gfunc:  Name (external) of user-supplied subroutine computing the 
c          value g(x): 
c                subroutine gfunc(m,n,t,x,gx)
c                real*8 t,x(n),g(m)
c                integer n,m
c                g(1)= ...
c            
C  work: real working array of dimension lwork    
C        real*8 work(lwork)
C      
C  LWORK: length of array work, hast to be at least 
C         13*nx+17*ng+(nx+ng)**2+nx**2+1
C
C  iwork: integer working array of dimension liwork
C         integer iwork(liwork)
C                           
C  LIWORK: length of array iwork, has to be at least nx+ng+4
C        
C**************************************************************
C      
C  output parameters:
C
C  T: endpoint of integration      
C        
C  X,V: values of position and volocity at T
C      
C  IDID =  1: integration succesfull
C       = -1: stop at t, h<hmin
C      
C  iwork(1): number of executed steps
C  iwork(2): number of accepted steps
C  iwork(3): number of rejected steps  
C
C  further statistical parameters:
C    lu-decompositions: iwork(2)
C    function calls (f,g) : iwork(1)*8     
C    function calls (m,gx): iwork(1)*7+iwork(2)     
C      
C**************************************************************
C
C  further input parameters     
C
C  work(1) : minimum stepsize, if iwork(1)<=0      
C            the code puts hmin=1.d-6
C      
C  iwork(4): maximum number of accepted steps,
C            if iwork(4)<=0 the code puts      
C            maxsteps=100000
C
C      
      
      
      
      
      
      
      
      IMPLICIT real*8 (A-H,O-Z)
      DIMENSION X(NX),V(NX)
      DIMENSION WORK(LWORK),IWORK(LIWORK)  
      EXTERNAL ffunc,mfunc,gxfunc,gfunc
C
C initialization
C
      nxg=ng+nx
C      
      if (work(1).le.0) then 
         hmin=1.e-6
      else
         hmin=work(1)
      end if   
C
      if (iwork(1).le.0) then 
         maxsteps=100000
      else
         maxsteps=iwork(1)
      end if      
C check input parameters      
      nwork=13*nx+17*ng+(nx+ng)**2+nx**2+1
      niwork=nx+ng+4
      
      
C prepare entries in work, iwork
      iyxi=2
      iyvi=iyxi+nx
      iywli=iyvi+nx
      ivs=iywli+nxg
      igvi=ivs+8*nxg
      iomat=igvi+8*ng
      iam=iomat+nxg*nxg
      iip=4
C Call the main integrator      
      CALL ROW(t0,te,nx,ng,nxg,h0,rtol,atol,itol,x,v,
     & hmin,
     & ffunc,mfunc,gxfunc,gfunc,
     & iwork(2),iwork(3),maxsteps, 
     & work(iyxi),work(iyvi),work(iywli),
     & work(ivs),work(igvi),work(iomat),work(iam),
     &  iwork(iip),idid)
      iwork(1)=iwork(2)+iwork(3)
      return 
      end
C      
C main integrator
C      
      subroutine row(t,tend,n,m,npm,h,rtol,atol,itol,yx,yv,
     & hmin,
     & ffunc,mfunc,gxfunc,gfunc,
     & naccept, nreject,maxsteps, 
     & yxi,yvi,ywli,vs,gvi,omat,am,ip,idid) 
       IMPLICIT real*8(a-h,o-z)
        DIMENSION YXI(N),YVI(N),YWLI(NPM)
        DIMENSION YX(N),YV(N)
        DIMENSION VS(NPM,8),GVI(M,8)
        DIMENSION AX(8,8),AV(8,8),AW(8,8),BX(8),BV(8),BXE(8)
        DIMENSION GAMMA(8,8),GAM(8),CC(8),GAMM(8)
        DIMENSION OMAT(NPM,NPM),IP(NPM)
        DIMENSION AM(N,N)
        logical laststep,repeated
        EXTERNAL FFUNC,MFUNC,GXFUNC,GFUNC
C
C ***************************************************
C
C  set the coefficients
C
C ***************************************************        
      CALL SETCOEFF(AX,AV,AW,BX,BV,CC,GAMMA,GAM,GAMM,bxe)
C ***************************************************
C
C     Main integration loop
C
C *******************************************************
      NP1=N+1
      naccept=0
      nreject=0
      Laststep=.false.
      repeated=.false.
 100  Continue
      if (t+h.lt.tend) goto 110
      h=tend-t
      laststep=.true.
 110  continue
C *******************************************************
C
C  basic integration step
C
C *******************************************************
      h2=h*h     
C *******************************************************
C
C  Computation and decomposition of optimization matrix
C  Computation of g_p*v for initial projection          
C
C *******************************************************          
      if (repeated) goto 120
      CALL MFUNC(N,NPM,T,YX,OMAT(1,1))
      CALL GXFUNC(M,N,NPM,T,YX,OMAT(NP1,1))
      do i=np1,npm 
        do j=1,n 
          omat(j,i)=omat(i,j)
        end do
        do j=np1,npm 
          omat(j,i)=0.0D0
        end do
      end do
C
      do i=1,npm 
        ywli(i)=0.0D0
      end do
      do j=1,n 
        do i=np1,npm 
          ywli(i)=ywli(i)-OMAT(i,j)*yv(j)
        end do
      end do 
C     
      Call DEC (NPM, NPM, OMAT, IP, IER)
C **********************
C      
C initial projection of v 
C
C **********************      
      CALL SOL(npm,npm,OMAT,ywli,ip)
      Do k=1,n
         YV(k)=YV(k)+ywli(k)
      end do
C ******************************************************
C      
C first stage
C      
C ******************************************************
 120  CALL FFUNC(N,T,YX,YV,VS(1,1))
      CALL GFUNC(M,N,T,YX,GVI(1,1))
      gamih2=-GAM(1)/h2
      do i=1,m        
        GVI(i,1)=gamih2*GVI(i,1)
        VS(i+n,1)=GVI(i,1)
      end do
      goto 130
C     
 130  CALL SOL(npm,npm,omat,vs(1,1),ip)
C ******************************************************
C      
C stages k>1
C      
C ******************************************************
      do i=2,8 
C ******************************************************
C     
C compute YXI,YVI,YWLI 
C      
C ******************************************************
        hci=h*cc(i)
        ti=t+hci
        do k=1,n 
          YXI(k)=hci*YV(k)+YX(k)
          YVI(k)=YV(k)
        end do  
        
        do j=1,i-1 
          axh2=ax(i,j)*h2
          avh=av(i,j)*h
          do k=1,n 
            YXI(k)=YXI(k)+axh2*VS(k,j)
            YVI(k)=YVI(k)+avh*VS(k,j)
          end do
        end do  
        
        do k=1,npm
          YWLI(k)=0.0D0
        end do
       
        do j=1,i-1 
          awij=aw(i,j)
          if (awij.NE.0.0D0) Then
            do k=1,npm
              YWLI(k)=YWLI(k)+awij*VS(k,j)
            end do
          end if 
        end do  
C ******************************************************
C       
C compute righthandside      
C      
C ******************************************************
        Call ffunc(N,Ti,YXI,YVI,VS(1,i))  
        CALL GFUNC(M,N,Ti,YXI,GVI(1,i))
        
        gamih2=-GAM(i)/h2
        do k=1,m        
          GVI(k,i)=gamih2*GVI(k,i)
        end do
        do j=1,i-1 
          gammaij=-gamma(i,j)
          
          do k=1,m 
            GVI(k,i)=GVI(k,i)+gammaij*GVI(k,j)
          end do
        end do
        do k=1,m 
          VS(n+k,i)=GVI(k,i)
        end do
C
        CALL MFUNC(N,N,Ti,YXI,AM)
        do j=1,n 
          wj=YWLI(j)
          do k=1,n 
            VS(k,i)=VS(k,i)-AM(k,j)*wj
          end do
        end do  
C    
        CALL GXFUNC(M,N,N,Ti,YXI,AM)
        do k=1,n 
          do j=1,m 
            VS(k,i)=VS(k,i)-AM(j,k)*YWLI(n+j)
          end do
        end do              
C ******************************************************
C
C   solve linear equation
C     
C ******************************************************
        CALL SOL(npm,npm,omat,VS(1,i),ip)
C
      end do
C****************************************
C
C  error estimation
C
C****************************************
      err=0.0d0
      do k=1,n
         y=0.0d0
         do j=1,8
            y=y+bxe(j)*vs(k,j)
         end do
         err=err+(y/(rtol*DABS(yx(k))+atol))**2.0
      end do
      err=h2*err**0.5
      if (err.gt.1) goto 300
C****************************************
C     accept step      
C****************************************
        t=t+h
        naccept=naccept+1
        repeated=.false.
        do k=1,n 
          YX(k)=h*YV(k)+YX(k)
        end do  
        
        do j=1,i-1 
          axh2=bx(j)*h2
          avh=bv(j)*h
          do k=1,n 
            YX(k)=YX(k)+axh2*VS(k,j)
            YV(k)=YV(k)+avh*VS(k,j)
          end do
       end do 
       if (.not.laststep) goto 310
       idid=0
       goto 400
 310   h=h*DMIN1(4.d0,.8d0/err**0.25)
       repeated=.false.
       goto 100
 300   continue
C****************************************
C step rejected
C****************************************
       nreject=nreject+1
       h=h*DMAX1(0.25D0,.8d0/err**0.25)
       laststep=.false.
       if (h.ge.hmin) goto 100
C stepsize too small       
       idid=-1
 400   continue
       return
       end
C
C      
C
      SUBROUTINE SETCOEFF(ax,av,aw,bx,bv,cc,gamma,gam,gamm,bxe)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION AX(8,8),AV(8,8),AW(8,8),BX(8),BV(8),BXE(8)
      DIMENSION GAMMA(8,8),GAM(8),GAMM(8),CC(8)
      DO I=1,8 
        DO J=1,8 
          AX(I,J)=0.0D0
          AV(I,J)=0.0D0
          AW(I,J)=0.0D0
          GAMMA(I,J)=0.0D0
        END DO
      END DO
      DO I=1,8 
        BX(I)=0.0D0
        BV(I)=0.0D0
        GAM(I)=0.0D0
        CC(I)=0.0D0
      END DO  

      AX( 2,1)=-5.308779335409982147453433753981e-01
      AX( 3,1)=5.000000000000000000000000000000e-01
      AX( 3,2)=1.000000000000000000000000000000e+00
      AX( 4,1)=1.249999999999999722444243843711e-01
      AX( 4,2)=3.299920321433769121455270578736e-01
      AX( 4,3)=-1.999800803584422109748786056116e-02
      AX( 5,1)=5.555555555555555247160270937457e-02
      AX( 5,2)=3.703703703703698640881114556578e-02
      AX( 5,3)=-6.172839506172798397509726697763e-03
      AX( 5,4)=2.469135802469129073455356149225e-02
      AX( 6,1)=3.472222222222222098864108374983e-01
      AX( 6,2)=2.314814814814823429323098480381e-01
      AX( 6,3)=7.716049382716026749928772687781e-02
      AX( 6,4)=3.858024691358084845571596588343e-02
      AX( 6,5)=-7.277815340696820560484336937179e-03
      AX( 7,1)=2.222222222222221821308352218693e-01
      AX( 7,2)=1.245901640109368935238620679229e-01
      AX( 7,3)=3.647035009329596677574159002688e-02
      AX( 7,4)=2.582473191217309407829993972427e-02
      AX( 7,5)=9.689526482936797191491962166765e-02
      AX( 7,6)=9.700346409435248173913635127974e-03
      AX( 8,1)=1.388888888888889332207110527406e-02
      AX( 8,2)=5.995895437615072001702465342987e-03
      AX( 8,3)=-1.454737842264335750952497505750e-03
      AX( 8,4)=4.452685561071953283807101087177e-03
      AX( 8,5)=1.343942753612087544212716494485e-02
      AX( 8,6)=1.156783645220162658195928173654e-03
      AX( 8,7)=-6.089372263294507738040639566179e-04
      AV( 2,1)=0.000000000000000000000000000000e+00
      AV( 3,1)=1.000000000000000000000000000000e+00
      AV( 3,2)=2.000000000000000000000000000000e+00
      AV( 4,1)=5.000000000000000000000000000000e-01
      AV( 4,2)=1.000000000000000000000000000000e+00
      AV( 4,3)=0.000000000000000000000000000000e+00
      AV( 5,1)=3.333333333333333148296162562474e-01
      AV( 5,2)=2.222222222222239029765233908620e-01
      AV( 5,3)=0.0D0
      AV( 5,4)=1.111111111111140470342206754140e-01
      AV( 6,1)=8.333333333333334813630699500209e-01
      AV( 6,2)=5.555555555555595770300669755670e-01
      AV( 6,3)=4.166666666666627993897975557047e-01
      AV( 6,4)=-1.388888888888832329193689929525e-01
      AV( 6,5)=-3.126173841186789814639013229680e-02
      AV( 7,1)=6.666666666666666296592325124948e-01
      AV( 7,2)=3.453405421989679258132355244015e-01
      AV( 7,3)=2.717741733449563490410127997166e-01
      AV( 7,4)=-9.910390224546974935471155276900e-02
      AV( 7,5)=4.186455269440948478987252201478e-01
      AV( 7,6)=4.080748915988272051968976938952e-02
      AV( 8,1)=1.666666666666664631257788187213e-01
      AV( 8,2)=-6.469966412196534455425478427060e-02
      AV( 8,3)=6.012760983875156672873174557026e-02
      AV( 8,4)=-9.247744189972938178012640264569e-02
      AV( 8,5)=6.737037296502148553400957098347e-01
      AV( 8,6)=1.127384916928332453389671741206e-01
      AV( 8,7)=1.314121067853177560191824113645e-01
      Aw( 2,1)=1.000000000000000000000000000000e+00
      Aw( 3,1)=1.000000000000000000000000000000e+00
      Aw( 3,2)=2.000000000000000000000000000000e+00
      Aw( 4,1)=1.000000000000000000000000000000e+00
      Aw( 4,2)=2.000000000000000000000000000000e+00
      Aw( 4,3)=0.000000000000000000000000000000e+00
      Aw( 5,1)=9.999999999999998889776975374843e-01
      Aw( 5,2)=6.666666666666566376520108860859e-01
      Aw( 5,3)=3.333333333333265979803172740503e-01
      Aw( 5,4)=0.0D0
      Aw( 6,1)=1.000000000000000000000000000000e+00
      Aw( 6,2)=6.666666666667380169997159100603e-01
      Aw( 6,3)=1.333333333333313053259416847141e+00
      Aw( 6,4)=-9.999999999999599209488110318489e-01
      Aw( 6,5)=-4.583276988957349362685533833428e-01
      Aw( 7,1)=9.999999999999997779553950749687e-01
      Aw( 7,2)=6.069985560983433003912068670616e-01
      Aw( 7,3)=1.029834055284167382282589642273e+00
      Aw( 7,4)=-7.263347772349871833696965950367e-01
      Aw( 7,5)=-8.250372705052208197051832883062e-03
      Aw( 7,6)=2.456922199871738529686204799418e-02
      Aw( 8,1)=1.000000000000000444089209850063e+00
      Aw( 8,2)=-5.074961563384361440398606646340e-01
      Aw( 8,3)=5.870814115027163104798546555685e-01
      Aw( 8,4)=-8.408294896719942235208122838230e-01
      Aw( 8,5)=-3.268819345212530258493188739521e-01
      Aw( 8,6)=2.541115605913969677231989408028e+00
      Aw( 8,7)=6.702017446050071214358467841521e+00
      GAMMA( 2,1)=1.530877933540998103723040912882e+00
      GAMMA( 3,1)=-2.000000000000000000000000000000e+00
      GAMMA( 3,2)=-4.000000000000000000000000000000e+00
      GAMMA( 4,1)=-1.999999999999999555910790149937e+00
      GAMMA( 4,2)=-5.919808771441045891492649388965e+00
      GAMMA( 4,3)=4.799521928602613063397086534678e-01
      GAMMA( 5,1)=2.666666666666277052399891545065e+00
      GAMMA( 5,2)=0.0D0
      GAMMA( 5,3)=0.0D0
      GAMMA( 5,4)=0.0D0
      GAMMA( 6,1)=-4.595238095237503372914034116548e+00
      GAMMA( 6,2)=-2.428571428572706025761362980120e+00
      GAMMA( 6,3)=1.214285714286230000169553022715e+00
      GAMMA( 6,4)=-2.428571428572512402865868352819e+00
      GAMMA( 6,5)=1.097241701278391801110956293996e+01
      GAMMA( 7,1)=9.058406218647874652560858521610e+00
      GAMMA( 7,2)=1.770853160187233765743286539873e-01
      GAMMA( 7,3)=-8.854265800933359964464131053319e-02
      GAMMA( 7,4)=1.770853160186926511521221527801e-01
      GAMMA( 7,5)=-1.140478687116817990698791618343e+00
      GAMMA( 7,6)=2.341000400905052236666392673214e-01
      GAMMA( 8,1)=3.145939777369665790729413856752e+00
      GAMMA( 8,2)=-1.428225684103834547400424526131e+00
      GAMMA( 8,3)=7.141128420520610475819012208376e-01
      GAMMA( 8,4)=-1.428225684104026838028289603244e+00
      GAMMA( 8,5)=3.418059458027589148088054571417e+00
      GAMMA( 8,6)=3.178771000370391508482725839713e+00
      GAMMA( 8,7)=5.681066068363909593585958646145e+00
      BX( 1)=4.999999999999993338661852249061e-01
      BX( 2)=2.582124245480629243232328917657e-01
      BX( 3)=2.042271210592968844199646127890e-01
      BX( 4)=-7.512090878526295201211837593291e-02
      BX( 5)=2.839215849034626426572458512965e-01
      BX( 6)=4.438074542425791274569135680395e-02
      BX( 7)=4.380403559510566424695099385644e-02
      BX( 8)=0.000000000000000000000000000000e+00
      BV( 1)=9.999999999999988897769753748435e-01
      BV( 2)=2.653343739036009285570116844610e-01
      BV( 3)=8.673328130482576003856820534565e-01
      BV( 4)=-7.346656260964796736345761019038e-01
      BV( 5)=-9.644493992316871544545620054123e-02
      BV( 6)=1.184466738971105437272512972413e+00
      BV( 7)=2.400672482016675601812494278420e+00
      BV( 8)=3.333333333333314274504743934813e-01
      GAM( 1)=-1.232061470943563286084554420086e+01
      GAM( 2)=1.000000000000000000000000000000e+00
      GAM( 3)=6.000000000000000000000000000000e+00
      GAM( 4)=2.400000000000000000000000000000e+01
      GAM( 5)=-5.399999999999097610725584672764e+01
      GAM( 6)=4.937142857142866603226138977334e+01
      GAM( 7)=-1.027302607822642102064492064528e+02
      GAM( 8)=-1.492764286952148609088908415288e+02
      CC( 1)=0.000000000000000000000000000000e+00
      CC( 2)=1.000000000000000000000000000000e+00
      CC( 3)=1.000000000000000000000000000000e+00
      CC( 4)=5.000000000000000000000000000000e-01
      CC( 5)=3.333333333333333148296162562474e-01
      CC( 6)=8.333333333333333703407674875052e-01
      CC( 7)=6.666666666666666296592325124948e-01
      CC( 8)=1.666666666666666574148081281237e-01
      BXE( 1)=0.0d0 
      BXE( 2)=  -5.5776949463803970e-01
      BXE( 3)=  2.7888474731896425e-01
      BXE( 4)= -5.5776949463795765e-01
      BXE( 5)=  3.0661177537300327e+00
      BXE( 6)=  -2.6888429388934698e-01
      BXE( 7)= -1.3324290350002405e+00
      BXE( 8)= -1.0570106296385264e-01

      RETURN
      END