text:
      subroutine foo
        print *, 'FAIL'
      end

ext: f

max_stars_repo_name: envp/fortran-src

max_issues_repo_name: envp/fortran-src

max_forks_repo_name: envp/fortran-src


text:
SUBROUTINE BMG3_SymStd_SETUP_TAUSCH(&
     NOG,&
     DimX, DimY, Dimz, DimXfine, DimYfine, DimZfine,&
     NProcI, NProcJ, NProcK &
     ) BIND(C, NAME='BMG3_SymStd_SETUP_Tausch')

  use ModInterface
  IMPLICIT NONE

  integer(C_INT), VALUE :: NOG, NProcI, NProcJ, NProcK
  INTEGER(len_t) :: DimX(NProcI,NOG), DimY(NProcJ,NOG), DimZ(NProcK,NOG)
  INTEGER(len_t) :: DimXfine(NProcI), DimYfine(NProcJ), DimZfine(NProcK)

  integer(len_t) :: I, J, K, iGs, jGs, kGs
  integer N

  DO I=1, NProcI
     DimX(I,NOG) = DimXfine(I)
  END DO

  DO J=1, NProcJ
     DimY(J,NOG) = DimYfine(J)
  END DO

  do k=1,nprock
     dimz(k,nog) = dimzfine(k)
  enddo

  DO N=NOG-1, 1, -1

     iGs = 1

     DO I=1, NProcI

        if (mod(iGs,2).eq.1) then
           DimX(I,N) = (DimX(I,N+1)+1)/2
        else
           IF (mod(DimX(I,N+1),2).eq.1) THEN
              DimX(I,N) = (DimX(I,N+1)-1)/2
           ELSE
              DimX(I,N) = (DimX(I,N+1)+1)/2
           ENDIF
        endif

        iGs = iGs + DimX(I,N+1)

     END DO

     jGs = 1

     DO J=1, NProcJ

        if (mod(jGs,2).eq.1) then
           DimY(J,N) = (DimY(J,N+1)+1)/2
        else
           IF (mod(DimY(J,N+1),2).eq.1) THEN
              DimY(J,N) = (DimY(J,N+1)-1)/2
           ELSE
              DimY(J,N) = (DimY(J,N+1)+1)/2
           ENDIF
        endif

        jGs = jGs + DimY(J,N+1)

     END DO

     kgs = 1

     do k=1,nprock

        if (mod(kgs,2) .eq. 1) then
           dimz(k,n) = (dimz(k,n+1)+1) / 2
        else
           if (mod(dimz(k,n+1),2) .eq. 1) then
              dimz(k,n) = (dimz(k,n+1)-1)/2
           else
              dimz(k,n) = (dimz(k,n+1)+1)/2
           endif
        endif

        kgs = kgs + dimz(k,n+1)
     enddo

  END DO

end SUBROUTINE BMG3_SymStd_SETUP_TAUSCH

max_stars_repo_name: cedar-framework/cedar

max_issues_repo_name: cedar-framework/cedar

max_forks_repo_name: cedar-framework/cedar


text:
Penny is a stinky dog who likes to hog the bed.

She likes eating gross things and then licking your face.

20121127 16:37:07 nbsp You are just stunning, Penny! Users/PeteB

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program problem52
  implicit none
  integer*8 :: i

  do i = 1,huge(i)
     if (predicate(i)) then
        exit
     end if
  end do

  print *, i

contains

  pure function score(n)
    implicit none
    integer*8, intent(in)             :: n
    integer*8                         :: temp,score
    integer, parameter, dimension(10) :: primes = [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 ]

    score=1
    temp=n
    do while (temp > 0)
       score=score*primes(mod(temp,10)+1)
       temp=temp/10
    end do

  end function score

  pure function predicate(n)
    implicit none
    integer*8, intent(in) :: n
    integer*8             :: temp1
    logical, dimension(5) :: temp2
    logical               :: predicate

    temp1=score(n)
    temp2=(/ (score(n*i) == temp1, i = 2, 6) /)
    predicate=all(temp2)

  end function predicate

end program problem52

max_stars_repo_name: jamesmcclain/Dioskouroi

max_issues_repo_name: jamesmcclain/Dioskouroi

max_forks_repo_name: jamesmcclain/Dioskouroi


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine nonlinmpc(co,vold,ipompc,nodempc,coefmpc,labmpc,
     &  nmpc,ikboun,ilboun,nboun,xbounact,aux,iaux,maxlenmpc,ikmpc,
     &  ilmpc,icascade,kon,ipkon,lakon,ne,reltime,newstep,xboun,fmpc,
     &  iit,idiscon,ncont,trab,ntrans,ithermal,mi)
!
!     updates the coefficients in nonlinear MPC's
!
      implicit none
!
c      logical isochoric
!
      character*8 lakon(*)
      character*20 labmpc(*),label
!
      integer ipompc(*),nodempc(3,*),irefnode,irotnode,idir,
     &  nmpc,index,ii,inode,node,id,ikboun(*),ilboun(*),nboun,
     &  i,j,k,idof,na,nb,nc,np,i1,i2,i3,iaux(*),maxlenmpc,n,
     &  l,m,lmax,mmax,ikmpc(*),ilmpc(*),icascade,neigh(7,8),
     &  mpc,kon(*),ipkon(*),indexe,ne,idofrem,idofins,nmpc0,nmpc01,
     &  newstep,iit,idiscon,ncont,iexpnode,indexexp,nmpcdif,ntrans,
     &  nodei,noded,lathyp(3,6),inum,ndir,number,ithermal,mi(*),
     &  newknot,indexexp1,indexexp2,indexexp3,idim
!
      real*8 co(3,*),coefmpc(*),vold(0:mi(2),*),c(3,3),dc(3,3,3),ww,
     &  e(3,3,3),d(3,3),w(3),f(3,3),c1,c2,c3,c4,c5,c6,xbounact(*),
     &  xboun(*),fmpc(*),expan,dd,a11,a12,a13,a21,a22,a23,a31,a32,a33,
     &  b11,b12,b13,b21,b22,b23,b31,b32,b33,aux(*),const,e1(3),e2(3),
     &  ddmax,a(3,3),b(3,3),xj,xi,et,ze,xlag(3,20),xeul(3,20),t1(3),
     &  coloc(3,8),reltime,csab(7),trab(7,*),pd(3),pi(3),e11(3,3),
     &  ad(3,3),ai(3,3),e22(3,3),e12(3,3),ru(3,3),ru1(3,3),
     &  ru2(3,3),ru3(3,3),u1(3,3),u2(3,3),u3(3,3),dcu(3,3,3),u(3,3),
     &  xi1,xi2,xi3,dco,dsi,dco2,dsi2
!
      data d /1.,0.,0.,0.,1.,0.,0.,0.,1./
      data e /0.,0.,0.,0.,0.,-1.,0.,1.,0.,
     &        0.,0.,1.,0.,0.,0.,-1.,0.,0.,
     &        0.,-1.,0.,1.,0.,0.,0.,0.,0./
      data neigh /1,9,2,12,4,17,5,2,9,1,10,3,18,6,
     &            3,11,4,10,2,19,7,4,11,3,12,1,20,8,
     &            5,13,6,16,8,17,1,6,13,5,14,7,18,2,
     &            7,15,8,14,6,19,3,8,15,7,16,5,20,4/
      data coloc /-1.,-1.,-1.,1.,-1.,-1.,1.,1.,-1.,-1.,1.,-1.,
     &            -1.,-1.,1.,1.,-1.,1.,1.,1.,1.,-1.,1.,1./
!
!     latin hypercube positions in a 3 x 3 matrix
!
      data lathyp /1,2,3,1,3,2,2,1,3,2,3,1,3,1,2,3,2,1/
!
c      irotnode=0
      irotnode=0
      irefnode=0
      if((icascade.eq.1).and.(newstep.ne.1).and.(ncont.eq.0)) icascade=0
c      isochoric=.false.
!
      ii=0
      loop: do
         ii=ii+1
         if(ii.gt.nmpc) exit
         if(labmpc(ii)(1:5).eq.'RIGID') then
!
            index=ipompc(ii)
            inode=nodempc(1,index)
            idir=nodempc(2,index)
            coefmpc(index)=1.d0
!
            index=nodempc(3,index)
            irefnode=nodempc(1,index)
            coefmpc(index)=-1.d0
!
            index=nodempc(3,index)
            node=nodempc(1,index)
!
!           check whether the rotational node is the same as in
!           the last rigid body MPC
!
            if(node.ne.irotnode) then
               irotnode=node
               w(1)=vold(1,node)
               w(2)=vold(2,node)
               w(3)=vold(3,node)
c               write(*,*) 'w   ',w(1),w(2),w(3)
               ww=dsqrt(w(1)*w(1)+w(2)*w(2)+w(3)*w(3))
!
               c1=dcos(ww)
               if(ww.gt.1.d-10) then
                  c2=dsin(ww)/ww
               else
                  c2=1.d0
               endif
               if(ww.gt.1.d-5) then
                  c3=(1.d0-c1)/ww**2
               else
                  c3=0.5d0
               endif
!
!              rotation matrix c
!
               do i=1,3
                  do j=1,3
                     c(i,j)=c1*d(i,j)+
     &                  c2*(e(i,1,j)*w(1)+e(i,2,j)*w(2)+e(i,3,j)*w(3))+
     &                  c3*w(i)*w(j)
                  enddo
               enddo
!
               c4=-c2
               if(ww.gt.0.00464159) then
                  c5=(ww*dcos(ww)-dsin(ww))/ww**3
               else
                  c5=-1.d0/3.d0
               endif
               if(ww.gt.0.0031623) then
                  c6=(ww*dsin(ww)-2.d0+2.d0*dcos(ww))/ww**4
               else
                  c6=-1.d0/12.d0
               endif
!
!              derivative of the rotation matrix c with respect to
!              the rotation vector w
!
               do i=1,3
                  do j=1,3
                     do k=1,3
                        dc(i,j,k)=c4*w(k)*d(i,j)+
     &                            c5*w(k)*(e(i,1,j)*w(1)+
     &                                     e(i,2,j)*w(2)+e(i,3,j)*w(3))+
     &                            c2*e(i,k,j)+
     &                            c6*w(k)*w(i)*w(j)+
     &                            c3*(d(i,k)*w(j)+d(j,k)*w(i))
                     enddo
                  enddo
               enddo
!
!              dummy variable
!
               do i=1,3
                  do j=1,3
c                     f(i,j)=c(i,j)-d(i,j)-dc(i,j,1)*w(1)-dc(i,j,2)*w(2)-
c     &                                    dc(i,j,3)*w(3)
                     f(i,j)=c(i,j)-d(i,j)
                  enddo
               enddo
            endif
!
!           determining the coefficients of the rotational degrees
!           of freedom
!
            coefmpc(index)=dc(idir,1,1)*(co(1,irefnode)-co(1,inode))+
     &           dc(idir,2,1)*(co(2,irefnode)-co(2,inode))+
     &           dc(idir,3,1)*(co(3,irefnode)-co(3,inode))
!
            index=nodempc(3,index)
            coefmpc(index)=dc(idir,1,2)*(co(1,irefnode)-co(1,inode))+
     &           dc(idir,2,2)*(co(2,irefnode)-co(2,inode))+
     &           dc(idir,3,2)*(co(3,irefnode)-co(3,inode))
!
            index=nodempc(3,index)
            coefmpc(index)=dc(idir,1,3)*(co(1,irefnode)-co(1,inode))+
     &           dc(idir,2,3)*(co(2,irefnode)-co(2,inode))+
     &           dc(idir,3,3)*(co(3,irefnode)-co(3,inode))
!
!           determining the nonhomogeneous part
!
            index=nodempc(3,index)
            coefmpc(index)=1.d0
!
!           old value of the nonhomogeneous term must be zero
!
            vold(nodempc(2,index),nodempc(1,index))=0.d0
            idof=8*(nodempc(1,index)-1)+nodempc(2,index)
            call nident(ikboun,idof,nboun,id)
            xbounact(ilboun(id))=f(idir,1)*(co(1,irefnode)-co(1,inode))+
     &           f(idir,2)*(co(2,irefnode)-co(2,inode))+
     &           f(idir,3)*(co(3,irefnode)-co(3,inode))-
     &           vold(idir,irefnode)+vold(idir,inode)
!
         elseif(labmpc(ii)(1:4).eq.'KNOT') then
!
!           dependent node
!
            index=ipompc(ii)
            inode=nodempc(1,index)
            idir=nodempc(2,index)
            coefmpc(index)=1.d0
!
!           translation node
!
            index=nodempc(3,index)
            node=nodempc(1,index)
            coefmpc(index)=-1.d0
!
!           check whether knot is the same as in the previous MPC
!
            if(node.ne.irefnode) then
               newknot=1
               irefnode=node
            else
               newknot=0
            endif
!
            read(labmpc(ii)(5:5),'(i1)') idim
!
!           expansion node
!
            index=nodempc(3,index)
            iexpnode=nodempc(1,index)
!
            if((idim.eq.1).or.(idim.eq.3)) then
!
!              nodes of knot lie on a straight line (1 term in MPC)
!
               indexexp=index
            elseif(idim.eq.2) then
!
!              node of knot lie in a plane (3 terms in MPC)
!
               indexexp1=index
               index=nodempc(3,index)
!
               indexexp2=index
               index=nodempc(3,index)
!
               indexexp3=index
!
            endif
!
            if(newknot.eq.1) then
               if((idim.eq.1).or.(idim.eq.3)) then
                  expan=1.d0+vold(1,iexpnode)
               elseif(idim.eq.2) then
                  xi1=vold(1,iexpnode)
                  xi2=1.d0+vold(2,iexpnode)
                  xi3=1.d0+vold(3,iexpnode)
                  dco=dcos(xi1)
                  dsi=dsin(xi1)
                  dco2=dcos(2.d0*xi1)
                  dsi2=dsin(2.d0*xi1)
                  dd=xi2**2-xi3**2
               endif
            endif
!
!           rotation node
!
            index=nodempc(3,index)
            irotnode=nodempc(1,index)
!
            if(newknot.eq.1) then
               w(1)=vold(1,irotnode)
               w(2)=vold(2,irotnode)
               w(3)=vold(3,irotnode)
               ww=dsqrt(w(1)*w(1)+w(2)*w(2)+w(3)*w(3))
!
               c1=dcos(ww)
               if(ww.gt.1.d-10) then
                  c2=dsin(ww)/ww
               else
                  c2=1.d0
               endif
               if(ww.gt.1.d-5) then
                  c3=(1.d0-c1)/ww**2
               else
                  c3=0.5d0
               endif
!
!              rotation matrix c
!
               do i=1,3
                  do j=1,3
                     c(i,j)=c1*d(i,j)+
     &                  c2*(e(i,1,j)*w(1)+e(i,2,j)*w(2)+e(i,3,j)*w(3))+
     &                  c3*w(i)*w(j)
                  enddo
               enddo
!
               c4=-c2
               if(ww.gt.0.00464159) then
                  c5=(ww*dcos(ww)-dsin(ww))/ww**3
               else
                  c5=-1.d0/3.d0
               endif
               if(ww.gt.0.0031623) then
                  c6=(ww*dsin(ww)-2.d0+2.d0*dcos(ww))/ww**4
               else
                  c6=-1.d0/12.d0
               endif
!
!              derivative of the rotation matrix c with respect to
!              the rotation vector w
!
               do i=1,3
                  do j=1,3
                     do k=1,3
                        dc(i,j,k)=c4*w(k)*d(i,j)+
     &                            c5*w(k)*(e(i,1,j)*w(1)+
     &                                     e(i,2,j)*w(2)+e(i,3,j)*w(3))+
     &                            c2*e(i,k,j)+
     &                            c6*w(k)*w(i)*w(j)+
     &                            c3*(d(i,k)*w(j)+d(j,k)*w(i))
                     enddo
                  enddo
               enddo
!
               if((idim.eq.1).or.(idim.eq.3)) then
!
!                 dummy variable for constant term
!
                  do i=1,3
                     do j=1,3
                        f(i,j)=expan*c(i,j)-d(i,j)
                     enddo
                  enddo
!
!                 derivative of the rotation matrix w.r.t the rotation
!                 vector multiplied by the expansion coefficient
!
                  do i=1,3
                     do j=1,3
                        do k=1,3
                           dc(i,j,k)=dc(i,j,k)*expan
                        enddo
                     enddo
                  enddo
               elseif(idim.eq.2) then
!
!                 local unit vectors
!
                  do i=1,3
                     t1(i)=co(i,irotnode)
                     e1(i)=co(i,iexpnode)
                  enddo
                  e2(1)=t1(2)*e1(3)-t1(3)*e1(2)
                  e2(2)=t1(3)*e1(1)-t1(1)*e1(3)
                  e2(3)=t1(1)*e1(2)-t1(2)*e1(1)
!
                  do i=1,3
                     do j=1,3
                        e11(i,j)=e1(i)*e1(j)
                        e22(i,j)=e2(i)*e2(j)
                        e12(i,j)=e1(i)*e2(j)+e2(i)*e1(j)
!
                        u(i,j)=t1(i)*t1(j)
     &                        +((xi2*dco)**2+(xi3*dsi)**2)*e11(i,j)
     &                        +((xi2*dsi)**2+(xi3*dco)**2)*e22(i,j)
     &                        +dd*dco*dsi*e12(i,j)
                        u1(i,j)=dd*(dco2*e12(i,j)
     &                              -dsi2*(e11(i,j)-e22(i,j)))
                        u2(i,j)=2.d0*xi2*(dco*dco*e11(i,j)
     &                         +dsi*dsi*e22(i,j))+xi2*dsi2*e12(i,j)
                        u3(i,j)=2.d0*xi3*(dsi*dsi*e11(i,j)
     &                         +dco*dco*e22(i,j))-xi3*dsi2*e12(i,j)
                     enddo
                  enddo
!
!                 calculating r.u, r.u1, r.u2 and r.u3
!
                  do i=1,3
                     do j=1,3
                        ru(i,j)=0.d0
                        ru1(i,j)=0.d0
                        ru2(i,j)=0.d0
                        ru3(i,j)=0.d0
!
                        do k=1,3
                           ru(i,j)=ru(i,j)+c(i,k)*u(k,j)
                           ru1(i,j)=ru1(i,j)+c(i,k)*u1(k,j)
                           ru2(i,j)=ru2(i,j)+c(i,k)*u2(k,j)
                           ru3(i,j)=ru3(i,j)+c(i,k)*u3(k,j)
                        enddo
                        f(i,j)=ru(i,j)-d(i,j)
                     enddo
                  enddo
!
!                 calculating dc.u
!
                  do i=1,3
                     do j=1,3
                        do k=1,3
                           dcu(i,j,k)=0.d0
                           do l=1,3
                              dcu(i,j,k)=dcu(i,j,k)+dc(i,l,k)*u(l,j)
                           enddo
                           dc(i,j,k)=dcu(i,j,k)
                        enddo
                     enddo
                  enddo
               endif
!
            endif
!
!           determining the coefficients of the expansion degrees
!           of freedom
!
            if((idim.eq.1).or.(idim.eq.3)) then
               coefmpc(indexexp)=c(idir,1)*(co(1,irefnode)-co(1,inode))+
     &           c(idir,2)*(co(2,irefnode)-co(2,inode))+
     &           c(idir,3)*(co(3,irefnode)-co(3,inode))
            elseif(idim.eq.2) then
!
!              if xi2=xi3 xi1 cannot be determined, since its coefficient
!              is always zero (cf. definition of u1)
!
               if(dabs(xi2-xi3).lt.1.d-10) then
                  coefmpc(indexexp1)=0.d0
c                  if(nodempc(2,indexexp1).ne.2) then
                     if(icascade.lt.1) icascade=1
                     nodempc(2,indexexp1)=2
c                  endif
               else
                  coefmpc(indexexp1)=ru1(idir,1)*
     &                      (co(1,irefnode)-co(1,inode))+
     &                 ru1(idir,2)*(co(2,irefnode)-co(2,inode))+
     &                 ru1(idir,3)*(co(3,irefnode)-co(3,inode))
c                  if(nodempc(2,indexexp1).ne.1) then
                     if(icascade.lt.1) icascade=1
                     nodempc(2,indexexp1)=1
c                  endif
               endif
               coefmpc(indexexp2)=ru2(idir,1)*
     &                      (co(1,irefnode)-co(1,inode))+
     &           ru2(idir,2)*(co(2,irefnode)-co(2,inode))+
     &           ru2(idir,3)*(co(3,irefnode)-co(3,inode))
               coefmpc(indexexp3)=ru3(idir,1)*
     &                      (co(1,irefnode)-co(1,inode))+
     &           ru3(idir,2)*(co(2,irefnode)-co(2,inode))+
     &           ru3(idir,3)*(co(3,irefnode)-co(3,inode))
            endif
!
!           determining the coefficients of the rotational degrees
!           of freedom
!
            coefmpc(index)=(dc(idir,1,1)*(co(1,irefnode)-co(1,inode))+
     &           dc(idir,2,1)*(co(2,irefnode)-co(2,inode))+
     &           dc(idir,3,1)*(co(3,irefnode)-co(3,inode)))
!
            index=nodempc(3,index)
            coefmpc(index)=(dc(idir,1,2)*(co(1,irefnode)-co(1,inode))+
     &           dc(idir,2,2)*(co(2,irefnode)-co(2,inode))+
     &           dc(idir,3,2)*(co(3,irefnode)-co(3,inode)))
!
            index=nodempc(3,index)
            coefmpc(index)=(dc(idir,1,3)*(co(1,irefnode)-co(1,inode))+
     &           dc(idir,2,3)*(co(2,irefnode)-co(2,inode))+
     &           dc(idir,3,3)*(co(3,irefnode)-co(3,inode)))
!
!           determining the nonhomogeneous part
!
            index=nodempc(3,index)
            coefmpc(index)=1.d0
!
!           old value of the nonhomogeneous term must be zero
!
            vold(nodempc(2,index),nodempc(1,index))=0.d0
            idof=8*(nodempc(1,index)-1)+nodempc(2,index)
            call nident(ikboun,idof,nboun,id)
            xbounact(ilboun(id))=f(idir,1)*(co(1,irefnode)-co(1,inode))+
     &           f(idir,2)*(co(2,irefnode)-co(2,inode))+
     &           f(idir,3)*(co(3,irefnode)-co(3,inode))-
     &           vold(idir,irefnode)+vold(idir,inode)
!
         elseif(labmpc(ii)(1:8).eq.'STRAIGHT') then
!
!           determining nodes and directions involved in MPC
!
            index=ipompc(ii)
            np=nodempc(1,index)
            j=nodempc(2,index)
            index=nodempc(3,index)
            i=nodempc(2,index)
            index=nodempc(3,index)
            na=nodempc(1,index)
            index=nodempc(3,nodempc(3,index))
            nb=nodempc(1,index)
!
!           determining the coefficients
!
            index=ipompc(ii)
            c2=co(i,na)+vold(i,na)-co(i,nb)-vold(i,nb)
            if(dabs(c2).lt.1.d-5) then
               write(*,*) '*WARNING in nonlinmpc: coefficient of'
               write(*,*)
     &     '          dependent node in STRAIGHT MPC is zero'
               idofrem=8*(np-1)+j
!
!              determining a new dependent term
!
               ddmax=abs(c2)
               l=i
               m=j
               do k=1,2
                  l=l+1
                  m=m+1
                  if(l.gt.3) l=l-3
                  if(m.gt.3) m=m-3
                  dd=dabs(co(l,na)+vold(l,na)-co(l,nb)-vold(l,nb))
                  if(dd.gt.ddmax) then
                     ddmax=dd
                     lmax=l
                     mmax=m
                  endif
               enddo
               i=lmax
               j=mmax
               idofins=8*(np-1)+j
!
               call changedepterm(ikmpc,ilmpc,nmpc,ii,idofrem,idofins)
!
               index=ipompc(ii)
               nodempc(2,index)=j
               index=nodempc(3,index)
               nodempc(2,index)=i
               index=nodempc(3,index)
               nodempc(2,index)=j
               index=nodempc(3,index)
               nodempc(2,index)=i
               index=nodempc(3,index)
               nodempc(2,index)=j
               index=nodempc(3,index)
               nodempc(2,index)=i
               index=nodempc(3,index)
               nodempc(2,index)=j
               if(icascade.eq.0) icascade=1
               c2=co(i,na)+vold(i,na)-co(i,nb)-vold(i,nb)
            endif
            coefmpc(index)=c2
            index=nodempc(3,index)
            c3=co(j,nb)+vold(j,nb)-co(j,na)-vold(j,na)
            coefmpc(index)=c3
            index=nodempc(3,index)
            c5=co(i,nb)+vold(i,nb)-co(i,np)-vold(i,np)
            coefmpc(index)=c5
            index=nodempc(3,index)
            c6=co(j,np)+vold(j,np)-co(j,nb)-vold(j,nb)
            coefmpc(index)=c6
            index=nodempc(3,index)
            c4=co(i,np)+vold(i,np)-co(i,na)-vold(i,na)
            coefmpc(index)=c4
            index=nodempc(3,index)
            c1=co(j,na)+vold(j,na)-co(j,np)-vold(j,np)
            coefmpc(index)=c1
            index=nodempc(3,index)
!
!           nonhomogeneous term
!
            coefmpc(index)=1.d0
!
!           old value of the nonhomogeneous term must be zero
!
            idof=8*(nodempc(1,index)-1)+nodempc(2,index)
            call nident(ikboun,idof,nboun,id)
            xbounact(ilboun(id))=-c1*c2+c3*c4
            if(newstep.eq.1) xboun(ilboun(id))=xbounact(ilboun(id))
            vold(nodempc(2,index),nodempc(1,index))=
     &          (1.d0-reltime)*xboun(ilboun(id))
         elseif(labmpc(ii)(1:5).eq.'PLANE') then
!
!           determining nodes and directions involved in MPC
!
            index=ipompc(ii)
            np=nodempc(1,index)
            i1=nodempc(2,index)
            index=nodempc(3,index)
            i2=nodempc(2,index)
            index=nodempc(3,index)
            i3=nodempc(2,index)
            index=nodempc(3,index)
            na=nodempc(1,index)
            index=nodempc(3,nodempc(3,nodempc(3,index)))
            nb=nodempc(1,index)
            index=nodempc(3,nodempc(3,nodempc(3,index)))
            nc=nodempc(1,index)
!
!           determining the coefficients
!
            a11=co(i1,np)+vold(i1,np)-co(i1,nc)-vold(i1,nc)
            a12=co(i2,np)+vold(i2,np)-co(i2,nc)-vold(i2,nc)
            a13=co(i3,np)+vold(i3,np)-co(i3,nc)-vold(i3,nc)
            a21=co(i1,na)+vold(i1,na)-co(i1,nc)-vold(i1,nc)
            a22=co(i2,na)+vold(i2,na)-co(i2,nc)-vold(i2,nc)
            a23=co(i3,na)+vold(i3,na)-co(i3,nc)-vold(i3,nc)
            a31=co(i1,nb)+vold(i1,nb)-co(i1,nc)-vold(i1,nc)
            a32=co(i2,nb)+vold(i2,nb)-co(i2,nc)-vold(i2,nc)
            a33=co(i3,nb)+vold(i3,nb)-co(i3,nc)-vold(i3,nc)
!
            b11=a22*a33-a23*a32
            b12=a31*a23-a21*a33
            b13=a21*a32-a31*a22
            b21=a32*a13-a12*a33
            b22=a11*a33-a31*a13
            b23=a31*a12-a11*a32
            b31=a12*a23-a22*a13
            b32=a21*a13-a11*a23
            b33=a11*a22-a12*a21
!
            index=ipompc(ii)
            if(dabs(b11).lt.1.d-5) then
               write(*,*) '*WARNING in nonlinmpc: coefficient of'
               write(*,*) '         dependent node in PLANE MPC is zero'
!
               idofrem=8*(nodempc(1,index)-1)+i1
!
               if(dabs(b12).gt.dabs(b13)) then
                  idofins=8*(nodempc(1,index)-1)+i2
                  call changedepterm
     &                   (ikmpc,ilmpc,nmpc,ii,idofrem,idofins)
                  coefmpc(index)=b12
                  nodempc(2,index)=i2
                  index=nodempc(3,index)
                  coefmpc(index)=b11
                  nodempc(2,index)=i1
                  index=nodempc(3,index)
                  coefmpc(index)=b13
                  index=nodempc(3,index)
                  coefmpc(index)=b22
                  nodempc(2,index)=i2
                  index=nodempc(3,index)
                  coefmpc(index)=b21
                  nodempc(2,index)=i1
                  index=nodempc(3,index)
                  coefmpc(index)=b23
                  index=nodempc(3,index)
                  coefmpc(index)=b32
                  nodempc(2,index)=i2
                  index=nodempc(3,index)
                  coefmpc(index)=b31
                  nodempc(2,index)=i1
                  index=nodempc(3,index)
                  coefmpc(index)=b33
                  index=nodempc(3,index)
                  coefmpc(index)=-b12-b22-b32
                  nodempc(2,index)=i2
                  index=nodempc(3,index)
                  coefmpc(index)=-b11-b21-b31
                  nodempc(2,index)=i1
                  index=nodempc(3,index)
                  coefmpc(index)=-b13-b23-b33
                  if(icascade.eq.0) icascade=1
               else
                  idofins=8*(nodempc(1,index)-1)+i3
                  call changedepterm
     &                  (ikmpc,ilmpc,nmpc,ii,idofrem,idofins)
                  coefmpc(index)=b13
                  nodempc(2,index)=i3
                  index=nodempc(3,index)
                  coefmpc(index)=b12
                  index=nodempc(3,index)
                  coefmpc(index)=b11
                  nodempc(2,index)=i1
                  index=nodempc(3,index)
                  coefmpc(index)=b23
                  nodempc(2,index)=i3
                  index=nodempc(3,index)
                  coefmpc(index)=b22
                  index=nodempc(3,index)
                  coefmpc(index)=b21
                  nodempc(2,index)=i1
                  index=nodempc(3,index)
                  coefmpc(index)=b33
                  nodempc(2,index)=i3
                  index=nodempc(3,index)
                  coefmpc(index)=b32
                  index=nodempc(3,index)
                  coefmpc(index)=b31
                  nodempc(2,index)=i1
                  index=nodempc(3,index)
                  coefmpc(index)=-b13-b23-b33
                  nodempc(2,index)=i3
                  index=nodempc(3,index)
                  coefmpc(index)=-b12-b22-b32
                  index=nodempc(3,index)
                  coefmpc(index)=-b11-b21-b31
                  nodempc(2,index)=i1
                  if(icascade.eq.0) icascade=1
               endif
            else
               coefmpc(index)=b11
               index=nodempc(3,index)
               coefmpc(index)=b12
               index=nodempc(3,index)
               coefmpc(index)=b13
               index=nodempc(3,index)
               coefmpc(index)=b21
               index=nodempc(3,index)
               coefmpc(index)=b22
               index=nodempc(3,index)
               coefmpc(index)=b23
               index=nodempc(3,index)
               coefmpc(index)=b31
               index=nodempc(3,index)
               coefmpc(index)=b32
               index=nodempc(3,index)
               coefmpc(index)=b33
               index=nodempc(3,index)
               coefmpc(index)=-b11-b21-b31
               index=nodempc(3,index)
               coefmpc(index)=-b12-b22-b32
               index=nodempc(3,index)
               coefmpc(index)=-b13-b23-b33
            endif
            index=nodempc(3,index)
            coefmpc(index)=1.d0
            idof=8*(nodempc(1,index)-1)+nodempc(2,index)
!
!           old value of the nonhomogeneous term must be zero
!
            call nident(ikboun,idof,nboun,id)
            xbounact(ilboun(id))=a11*b11+a12*b12+a13*b13
            if(newstep.eq.1) xboun(ilboun(id))=xbounact(ilboun(id))
            vold(nodempc(2,index),nodempc(1,index))=0.d0
c         elseif(labmpc(ii)(1:9).eq.'ISOCHORIC') then
c            isochoric=.true.
c!
c!           next segment is deactivated (CYCLID instead of CYCLIC):
c!           cylic MPC's are considered to be linear
c!
c         elseif((labmpc(ii)(1:6).eq.'CYCLID').and.(ithermal.ne.2)) then
c            index=ipompc(ii)
c            noded=nodempc(1,index)
c!
c!           check for thermal MPC
c!
c            if(nodempc(2,index).eq.0) cycle loop
c!
c!           check whether the next two MPC's are cyclic MPC's
c!           applied to the same dependent node
c!
c            if((nodempc(1,ipompc(ii+1)).ne.noded).or.
c     &         (labmpc(ii+1)(1:6).ne.'CYCLIC').or.
c     &         (nodempc(1,ipompc(ii+2)).ne.noded).or.
c     &         (labmpc(ii+2)(1:6).ne.'CYCLIC')) then
c               write(*,*) '*WARNING in nonlinmpc: no three'
c               write(*,*) '         cyclic MPCs pertaining'
c               write(*,*) '         to the same dependent node;'
c               write(*,*) '         no update'
c               cycle loop
c            endif
c!
c!           finding the cyclic symmetry axis
c!
c            do i=1,ntrans
c               if(trab(7,i).eq.2) exit
c            enddo
c            if(i.gt.ntrans) then
c               write(*,*) '*ERROR in nonlinmpc: cyclic symmetry'
c               write(*,*) '       axis not found'
c               call exit(201)
c            endif
c            do j=1,6
c               csab(j)=trab(j,i)
c            enddo
c            csab(7)=-1
c!
c!           determining the independent node
c!
c            nodei=0
c            do
c               if(nodempc(1,index).ne.noded) then
c                  if(nodei.eq.0) then
c                     nodei=nodempc(1,index)
c                  elseif(nodei.ne.nodempc(1,index)) then
c                     write(*,*) '*WARNING in nonlinmpc:'
c                     write(*,*) '          cyclic symmetry conditions'
c                     write(*,*) '          between unequal meshes'
c                     write(*,*) '          no update'
c                     cycle loop
c                  endif
c               endif
c               index=nodempc(3,index)
c               if(index.eq.0) then
c                  if(nodei.eq.0) then
c                     write(*,*) '*ERROR in nonlinmpc:'
c                     write(*,*) '       no independent node found'
c                     call exit(201)
c                  else
c                     exit
c                  endif
c               endif
c            enddo
c!
c!           actual location of dependent and independent node
c!
c            do i=1,3
c               pd(i)=co(i,noded)+vold(i,noded)
c               pi(i)=co(i,nodei)+vold(i,nodei)
c            enddo
c!
c!           update transformation matrix
c!
c            call transformatrix(csab,pd,ad)
c            call transformatrix(csab,pi,ai)
c!
c!     checking for latin hypercube positions in matrix al none of
c!     which are zero
c!
c            do inum=1,6
c               if((dabs(ad(lathyp(1,inum),1)).gt.1.d-3).and.
c     &            (dabs(ad(lathyp(2,inum),2)).gt.1.d-3).and.
c     &            (dabs(ad(lathyp(3,inum),3)).gt.1.d-3)) exit
c            enddo
c!
c!           remove old DOFs
c!
c            do j=1,3
c               idof=8*(noded-1)+j
c               call nident(ikmpc,idof,nmpc,id)
c               if(id.lt.0) then
c                  write(*,*) '*ERROR in nonlinmpc: error in'
c                  write(*,*) '       MPC database'
c                  call exit(201)
c               elseif(ikmpc(id).ne.idof) then
c                  write(*,*) '*ERROR in nonlinmpc: error in'
c                  write(*,*) '       MPC database'
c                  call exit(201)
c               endif
c!
c               do k=id,nmpc-1
c                  ikmpc(k)=ikmpc(k+1)
c                  ilmpc(k)=ilmpc(k+1)
c               enddo
c            enddo
c!
c!           add new MPCs
c!
c            ii=ii-1
c            do ndir=1,3
c               ii=ii+1
c               number=lathyp(ndir,inum)
c               idof=8*(noded-1)+number
c               call nident(ikmpc,idof,nmpc-1,id)
c               if(id.gt.0) then
c                  if(ikmpc(id).eq.idof) then
c                     write(*,*) '*WARNING in nonlinmpc: cyclic MPC
c     & in node'
c                     write(*,*) '         ',noded,' and direction ',ndir
c                     write(*,*) '         cannot be created: the'
c                     write(*,*) '         DOF in this node is already us
c     &ed'
c                     cycle
c                  endif
c               endif
c               number=number-1
c!
c!     updating ikmpc and ilmpc
c!
c               do j=nmpc,id+2,-1
c                  ikmpc(j)=ikmpc(j-1)
c                  ilmpc(j)=ilmpc(j-1)
c               enddo
c               ikmpc(id+1)=idof
c               ilmpc(id+1)=nmpc
c!
c!              update the MPC coefficients
c!
c               index=ipompc(ii)
c               do j=1,3
c                  number=number+1
c                  if(number.gt.3) number=1
c                  if(dabs(ad(number,ndir)).lt.1.d-5) cycle
c                  if(index.eq.0) then
c                     write(*,*)'*ERROR in nonlinmpc: index=0'
c                     call exit(201)
c                  endif
c                  nodempc(1,index)=noded
c                  nodempc(2,index)=number
c                  coefmpc(index)=ad(number,ndir)
c                  index=nodempc(3,index)
c               enddo
c               do j=1,3
c                  number=number+1
c                  if(number.gt.3) number=1
c                  if(dabs(ai(number,ndir)).lt.1.d-5) cycle
c                  if(index.eq.0) then
c                     write(*,*)'*ERROR in nonlinmpc: index=0'
c                     call exit(201)
c                  endif
c                  nodempc(1,index)=nodei
c                  nodempc(2,index)=number
c                  coefmpc(index)=-ai(number,ndir)
c                  index=nodempc(3,index)
c               enddo
c            enddo
         elseif((labmpc(ii)(1:20).ne.'                    ').and.
     &          (labmpc(ii)(1:10).ne.'PRETENSION').and.
     &          (labmpc(ii)(1:7).ne.'CONTACT').and.
     &          (labmpc(ii)(1:7).ne.'NETWORK').and.
     &          (labmpc(ii)(1:5).ne.'FLUID').and.
     &          (labmpc(ii)(1:6).ne.'CYCLIC').and.
     &          (labmpc(ii)(1:9).ne.'SUBCYCLIC')) then
            index=ipompc(ii)
            i=0
            do
               if(index.eq.0) exit
               node=nodempc(1,index)
               i=i+1
               iaux(i)=nodempc(2,index)
               iaux(maxlenmpc+i)=node
               aux(6*maxlenmpc+i)=coefmpc(index)
               do j=1,3
                  aux(3*(i-1)+j)=co(j,node)
                  aux(3*(maxlenmpc+i-1)+j)=vold(j,node)
               enddo
               index=nodempc(3,index)
            enddo
            n=i-1
            if((labmpc(ii)(1:7).eq.'MEANROT').or.
     &         (labmpc(ii)(1:1).eq.'1')) then
               call umpc_mean_rot(aux,aux(3*maxlenmpc+1),const,
     &            aux(6*maxlenmpc+1),iaux,n,fmpc(ii),iit,idiscon,
     &            iaux(maxlenmpc+1),ikmpc,nmpc,ikboun,nboun,
     &            labmpc(ii))
            elseif(labmpc(ii)(1:4).eq.'DIST') then
               call umpc_dist(aux,aux(3*maxlenmpc+1),const,
     &            aux(6*maxlenmpc+1),iaux,n,fmpc(ii),iit,idiscon)
            elseif(labmpc(ii)(1:3).eq.'GAP') then
               call umpc_gap(aux,aux(3*maxlenmpc+1),const,
     &            aux(6*maxlenmpc+1),iaux,n,fmpc(ii),iit,idiscon)
            elseif(labmpc(ii)(1:4).eq.'USER') then
               call umpc_user(aux,aux(3*maxlenmpc+1),const,
     &            aux(6*maxlenmpc+1),iaux,n,fmpc(ii),iit,idiscon)
            else
               write(*,*) '*ERROR in nonlinmpc: mpc of type ',labmpc(ii)
               write(*,*) '       is unknown'
               call exit(201)
            endif
            index=ipompc(ii)
!
            if((iaux(1).ne.nodempc(2,index)).or.
     &         (iaux(maxlenmpc+1).ne.nodempc(1,index))) then
!
!              dependent MPC has changed
!
               idofrem=8*(nodempc(1,index)-1)+nodempc(2,index)
               idofins=8*(nodempc(1,index)-1)+iaux(1)
               call changedepterm(ikmpc,ilmpc,nmpc,ii,idofrem,idofins)
               if(icascade.eq.0) icascade=1
            endif
!
            i=0
            do
               if(index.eq.0) exit
               i=i+1
               if(i.le.n) then
!
!                 check whether any directions have changed:
!                 necessitates calling of remastruct
!
                  if(iaux(i).ne.nodempc(2,index)) then
                     if(icascade.eq.0) icascade=1
                  endif
                  nodempc(2,index)=iaux(i)
                  coefmpc(index)=aux(6*maxlenmpc+i)
               else
                  coefmpc(index)=1.d0
!
!                 old value of the nonhomogeneous term must be zero
!
                  vold(nodempc(2,index),nodempc(1,index))=0.d0
                  idof=8*(nodempc(1,index)-1)+nodempc(2,index)
                  call nident(ikboun,idof,nboun,id)
                  xbounact(ilboun(id))=const
               endif
               index=nodempc(3,index)
            enddo
         endif
      enddo loop
c!
c!     incompressible material
c!
c      if(.not.isochoric) return
c!
c!     initialization of the mpc's
c!
c      nmpc01=0
c      nmpcdif=0
c      do i=1,nmpc
c         if(labmpc(i)(1:9).eq.'ISOCHORIC') then
c            if(nmpc01.eq.0) nmpc01=i
c            nmpcdif=i
c            index=ipompc(i)
c            do
c               if(nodempc(3,index).eq.0) then
c                  idof=8*(nodempc(1,index)-1)+nodempc(2,index)
c                  call nident(ikboun,idof,nboun,id)
c                  xbounact(ilboun(id))=0.d0
c                  exit
c               endif
c               coefmpc(index)=0.d0
c               index=nodempc(3,index)
c            enddo
c         endif
c      enddo
c      nmpc0=nmpc01-1
c      nmpcdif=nmpcdif-nmpc0
c!
c      do i=1,ne
c         if(ipkon(i).lt.0) cycle
c         if(lakon(i)(1:7).eq.'C3D20RI') then
c            indexe=ipkon(i)
c!
c            do j=1,20
c               node=kon(indexe+j)
c               do k=1,3
c                  xlag(k,j)=co(k,node)
c                  xeul(k,j)=xlag(k,j)+vold(k,node)
c               enddo
c            enddo
c!
c            do j=1,8
c               mpc=0
c               node=kon(indexe+j)
c               label(1:9)='ISOCHORIC'
c               write(label(10:20),'(i11)') node
cc               write(*,*) 'nonlinmpclab ',label
c               call cident20(labmpc(nmpc01),label,nmpcdif,id)
c               id=id+nmpc0
cc               write(*,*) 'nonlinmpclab ',id,label,labmpc(id)
c               if(id.gt.0) then
c                  if(labmpc(id).eq.label) then
c                     mpc=id
c                  endif
c               endif
c               if(mpc.eq.0) cycle
c!
c               xi=coloc(1,j)
c               et=coloc(2,j)
c               ze=coloc(3,j)
c!
c               call deuldlag(xi,et,ze,xlag,xeul,xj,a)
c!
c               b(1,1)=a(2,2)*a(3,3)-a(2,3)*a(3,2)
c               b(1,2)=a(3,1)*a(2,3)-a(2,1)*a(3,3)
c               b(1,3)=a(2,1)*a(3,2)-a(3,1)*a(2,2)
c               b(2,1)=a(3,2)*a(1,3)-a(1,2)*a(3,3)
c               b(2,2)=a(1,1)*a(3,3)-a(3,1)*a(1,3)
c               b(2,3)=a(3,1)*a(1,2)-a(1,1)*a(3,2)
c               b(3,1)=a(1,2)*a(2,3)-a(2,2)*a(1,3)
c               b(3,2)=a(2,1)*a(1,3)-a(1,1)*a(2,3)
c               b(3,3)=a(1,1)*a(2,2)-a(1,2)*a(2,1)
c!
c               index=ipompc(mpc)
c               do
c                  if(nodempc(3,index).eq.0) then
c                     coefmpc(index)=1.d0
c                     idof=8*(nodempc(1,index)-1)+nodempc(2,index)
c                     call nident(ikboun,idof,nboun,id)
c                     xbounact(ilboun(id))=xbounact(ilboun(id))+
c     &                    a(1,1)*b(1,1)+a(1,2)*b(1,2)+a(1,3)*b(1,3)
c     &                    -1.d0/xj
cc                     write(*,*) 'nonlinmpcboun ',nodempc(1,index),
cc     &                        nodempc(2,index),ilboun(id),
cc     &                       xbounact(ilboun(id))
c                     exit
c                  else
c                     node=nodempc(1,index)
c                     idir=nodempc(2,index)
c                     do k=1,7
c                        if(kon(indexe+neigh(k,j)).eq.node) then
c                           if(k.eq.1) then
c                              if(idir.eq.1) then
c                                 coefmpc(index)=coefmpc(index)+
c     &                             1.5d0*(xi*b(1,1)+et*b(1,2)+ze*b(1,3))
c                              elseif(idir.eq.2) then
c                                 coefmpc(index)=coefmpc(index)+
c     &                             1.5d0*(xi*b(2,1)+et*b(2,2)+ze*b(2,3))
c                              elseif(idir.eq.3) then
c                                 coefmpc(index)=coefmpc(index)+
c     &                             1.5d0*(xi*b(3,1)+et*b(3,2)+ze*b(3,3))
c                              endif
c                           elseif(k.eq.2) then
c                              if(idir.eq.1) then
c                               coefmpc(index)=coefmpc(index)-
c     &                                        2.d0*xi*b(1,1)
c                              elseif(idir.eq.2) then
c                               coefmpc(index)=coefmpc(index)-
c     &                                        2.d0*xi*b(2,1)
c                              elseif(idir.eq.3) then
c                               coefmpc(index)=coefmpc(index)-
c     &                                        2.d0*xi*b(3,1)
c                              endif
c                           elseif(k.eq.3) then
c                              if(idir.eq.1) then
c                              coefmpc(index)=coefmpc(index)+
c     &                                       0.5d0*xi*b(1,1)
c                              elseif(idir.eq.2) then
c                              coefmpc(index)=coefmpc(index)+
c     &                                       0.5d0*xi*b(2,1)
c                              elseif(idir.eq.3) then
c                              coefmpc(index)=coefmpc(index)+
c     &                                       0.5d0*xi*b(3,1)
c                              endif
c                           elseif(k.eq.4) then
c                              if(idir.eq.1) then
c                               coefmpc(index)=coefmpc(index)-
c     &                                        2.d0*et*b(1,2)
c                              elseif(idir.eq.2) then
c                               coefmpc(index)=coefmpc(index)-
c     &                                        2.d0*et*b(2,2)
c                              elseif(idir.eq.3) then
c                               coefmpc(index)=coefmpc(index)-
c     &                                        2.d0*et*b(3,2)
c                              endif
c                           elseif(k.eq.5) then
c                              if(idir.eq.1) then
c                              coefmpc(index)=coefmpc(index)+
c     &                                       0.5d0*et*b(1,2)
c                              elseif(idir.eq.2) then
c                              coefmpc(index)=coefmpc(index)+
c     &                                       0.5d0*et*b(2,2)
c                              elseif(idir.eq.3) then
c                              coefmpc(index)=coefmpc(index)+
c     &                                       0.5d0*et*b(3,2)
c                              endif
c                           elseif(k.eq.6) then
c                              if(idir.eq.1) then
c                              coefmpc(index)=coefmpc(index)-
c     &                                       2.d0*ze*b(1,3)
c                              elseif(idir.eq.2) then
c                              coefmpc(index)=coefmpc(index)-
c     &                                       2.d0*ze*b(2,3)
c                              elseif(idir.eq.3) then
c                              coefmpc(index)=coefmpc(index)-
c     &                                       2.d0*ze*b(3,3)
c                              endif
c                           elseif(k.eq.7) then
c                              if(idir.eq.1) then
c                              coefmpc(index)=coefmpc(index)+
c     &                                       0.5d0*ze*b(1,3)
c                              elseif(idir.eq.2) then
c                              coefmpc(index)=coefmpc(index)+
c     &                                       0.5d0*ze*b(2,3)
c                              elseif(idir.eq.3) then
c                              coefmpc(index)=coefmpc(index)+
c     &                                       0.5d0*ze*b(3,3)
c                              endif
c                           endif
c                           exit
c                        endif
c                     enddo
c                  endif
c                  index=nodempc(3,index)
c               enddo
c!
c            enddo
c         endif
c      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


submodule(forlab_linalg) forlab_linalg_outer

    implicit none

contains

    module procedure outer_int8
    integer :: m, n

    m = size(x)
    n = size(y)
    allocate (outer_int8, &
              source=spread(x, 2, n)*spread(y, 1, m))

    end procedure outer_int8
    module procedure outer_int16
    integer :: m, n

    m = size(x)
    n = size(y)
    allocate (outer_int16, &
              source=spread(x, 2, n)*spread(y, 1, m))

    end procedure outer_int16
    module procedure outer_int32
    integer :: m, n

    m = size(x)
    n = size(y)
    allocate (outer_int32, &
              source=spread(x, 2, n)*spread(y, 1, m))

    end procedure outer_int32
    module procedure outer_int64
    integer :: m, n

    m = size(x)
    n = size(y)
    allocate (outer_int64, &
              source=spread(x, 2, n)*spread(y, 1, m))

    end procedure outer_int64

    module procedure outer_sp
    integer :: m, n

    m = size(x)
    n = size(y)
    allocate (outer_sp, &
              source=spread(x, 2, n)*spread(y, 1, m))

    end procedure outer_sp
    module procedure outer_dp
    integer :: m, n

    m = size(x)
    n = size(y)
    allocate (outer_dp, &
              source=spread(x, 2, n)*spread(y, 1, m))

    end procedure outer_dp

end submodule forlab_linalg_outer


text:
*     -------------------------------------------------------------
*
*     A PROGRAM FOR ROTATING A SET OF RADIAL ORBITALS
*
*                By P. Jonsson and M. Godefroid
*
*                1993
*
*     ------------------------------------------------------------
*
      PROGRAM ROTATE
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (NOD=220,NWD=60)
      CHARACTER*1 END,PP,ORBSYM
      CHARACTER*3 EL(75),ELC(3)
      CHARACTER*6 ATOM,TERM,ALABEL(2)
      COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,MASS,NCFG,IB,IC,ID
     :   ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,FINE,NSCF,NCLOSD,RMASS
      COMMON /RADIAL/R(NOD),RR(NOD),R2(NOD),P(NOD,NWD),YK(NOD),
     :   YR(NOD),X(NOD),AZ(NWD),L(NWD),MAX(NWD),N(NWD)
      DIMENSION IPT(3),ET(75),EK(75),PROT(NOD)
      LOGICAL PRINT
      PARAMETER (IREAD=5,IWRITE=6)
*
*
C
C     Open wfn files
C
      OPEN(UNIT=1,FILE='wfn.inp',STATUS='OLD',
     :     FORM='UNFORMATTED')
      OPEN(UNIT=2,FILE='wfn.rot',STATUS='UNKNOWN',
     :     FORM='UNFORMATTED')

      CALL INITR
C      write(0,*) ' Enter tolerance for keeping an orbital'
C      read(5,*) tol
C
      WRITE(*,*)
     : ' Enter orbitals symmetry (s,p,d ..) to be rotated'
      READ(*,'(A1)') ORBSYM
      LROT= LVAL(ORBSYM)
      write(*,*)
      WRITE(*,*) ' Original orbitals : '
      write(*,*) ' ------------------- '
      write(*,*)
C
C     Read the radial functions.
C
      NUMBERW = 0
      I = 1
111   READ(1,END=112) ATOM,TERM,EL(I),MR,Z,ET(I),EK(I),AZ(I),
     :     (P(J,I),J=1,MR)
      IF (EL(I)(1:1) .NE. ' ') THEN
         N(I) = ICHAR(EL(I)(1:1)) - ICHAR('1') + 1
         L(I) = LVAL(EL(I)(2:2))
      ELSE
         N(I) = ICHAR(EL(I)(2:2)) - ICHAR('1') + 1
         L(I) = LVAL(EL(I)(3:3))
      ENDIF
      NUMBERW = NUMBERW + 1
      MM = MR+1
      DO 113 J = MM,NO
         P(J,I) = D0
113   CONTINUE
      MAX(I) = MR
C
C     Initialize radial grid
C
      IF (I.EQ.1) THEN
         DO 2 J=1,NO
            R(J)=EXP(RHO)/Z
            RR(J)=R(J)*R(J)
            R2(J)=SQRT(R(J))
            RHO=RHO+H
2        CONTINUE
      ENDIF
      rad = qquadr(i,i,1)
      write(*,1001) el(i),rad,qquadr(i,i,0)
      I = I + 1
      GOTO 111
112   CONTINUE
C
C     Find first orbital of right symmetry
C
      write(*,*) ' How many do you want to skip ? '
      read(*,*) iskip
      IC = 0
      DO 114 J = 1,NUMBERW
         IF (L(J).EQ.LROT) THEN
            IFIRST=J
            IC = IC + 1
            IF (IC.EQ.iskip+1) GOTO 126
         ENDIF
114   CONTINUE
126   write(*,1002) ifirst
C
C     Rotate the orbital
C
      nfo = 1
      DO 115 J = IFIRST+1,NUMBERW
         IF (L(J).EQ.LROT) THEN
            nfo = nfo + 1
            DO 116 K = 1,NO
               P(K,IFIRST) = P(K,IFIRST) + P(K,J)
116         CONTINUE
            AZ(IFIRST) = AZ(IFIRST) +  AZ(J)
            maxrot = max0(max(ifirst),max(j))
         ENDIF
115   CONTINUE
      print*,' coucou '
      if (nfo.eq.1) then
        write(*,*) ' Only one orbital of this l-symmetry... stop...'
        call exit
      end if
      print*,' nfo = ',nfo,' el = ',el,' ifirst = ',ifirst
      WRITE(*,1004) el(ifirst),nfo-1
      write(*,*)
      WRITE(*,*) ' Rotated  orbitals : '
      write(*,*) ' ------------------- '
      write(*,*)
C
C     Normalize the rotated orbital and save it
C
      QNORM = DSQRT(QQUADR(IFIRST,IFIRST,0))
      DO 117 J = 1,NO
         P(J,IFIRST) = P(J,IFIRST)/QNORM
         PROT(J) = P(J,IFIRST)
117   CONTINUE
      AZ(IFIRST) = AZ(IFIRST)/DSQRT(DFLOAT(nfo))
      AZROT = AZ(IFIRST)
      QNORM = DSQRT(QQUADR(IFIRST,IFIRST,0))
C
C     Read the radial functions again and replace the
C     ifirst orbital with the rotated one.
C
      REWIND(UNIT=1)
      I = 1
      M = 1
   12 READ(1,END=13) ATOM,TERM,EL(I),MR,Z,ET(I),EK(I),AZ(I),
     :     (P(J,I),J=1,MR)
      IF (EL(I)(1:1) .NE. ' ') THEN
         N(I) = ICHAR(EL(I)(1:1)) - ICHAR('1') + 1
         L(I) = LVAL(EL(I)(2:2))
      ELSE
         N(I) = ICHAR(EL(I)(2:2)) - ICHAR('1') + 1
         L(I) = LVAL(EL(I)(3:3))
      ENDIF
      MM = MR+1
      DO 24 J = MM,NO
24    P(J,I) = D0
      MAX(I)= MR

      IF (I.EQ.IFIRST) THEN
         DO 118 J = 1,NO
            P(J,I) = PROT(J)
118      CONTINUE
         AZ(I) = AZROT
         MAX(I) = maxrot
      ENDIF
C
C  *****  ORTHOGONALIZE TO INNER FUNCTIONS
C
      IM = I - 1
      DO 6 II =1,IM
      IF (L(i) .ne. L(ii)) GO TO 6
      If( (EL(i)(1:1) .NE. ' ' .and. EL(ii)(1:1) .NE. ' ')
     :  .and. (EL(i)(3:3) .NE. EL(ii)(3:3))) Go to 6
      PN = QQUADR(I,II,0)
      IF (ABS(ABS(PN) - 1.d0) .lt. 1.D-08) GO TO 12
      IF ( DABS(PN) .GT. 1.D-08 ) THEN
         PNN = DSQRT(1.d0-PN*PN)
         IF (P(50,I) - PN*P(50,II) .LT. D0) PNN = -PNN
         write(*,*) i,ii,pn
	 mr = max0(max(i),max(ii))
         DO 25 J = 1,MR
            P(J,I) =(P(J,I) - PN*P(J,II))/PNN
25       CONTINUE
         AZ(I) = (AZ(I) - PN*AZ(II))/PNN
      END IF
6     CONTINUE
       pnn = QQUADR(I,I,0)
       pnn = dsqrt(pnn)
C       if (pnn .lt. tol) go to 9
       do 30 j = 1, MR
	  P(J,I) = P(J,I)/pnn
30     continue
       pnn = QQUADR(I,I,0)
       rad = qquadr(i,i,1)
       write(*,*) ' Orbital ',el(i), ' with radius',rad, ' <|> = ',pnn
       WRITE(2) ATOM,TERM,EL(I),MR,Z,ET(I),EK(I),AZ(I),
     :     (P(J,I),J=1,MR)
       I = I + 1
9     IF ( I .LT. NWD) THEN
	GO TO 12
      else
	write(*,*) ' TOO many orbitals'
      end if
13    continue

      REWIND(UNIT=2)
      I = 1
      M = 1
  912 READ(2,END=913) ATOM,TERM,EL(I),MR,Z,ET(I),EK(I),AZ(I),
     :     (P(J,I),J=1,MR)
      IF (EL(I)(1:1) .NE. ' ') THEN
         N(I) = ICHAR(EL(I)(1:1)) - ICHAR('1') + 1
         L(I) = LVAL(EL(I)(2:2))
      ELSE
         N(I) = ICHAR(EL(I)(2:2)) - ICHAR('1') + 1
         L(I) = LVAL(EL(I)(3:3))
      ENDIF
      MM = MR+1
      DO 924 J = MM,NO
924    P(J,I) = D0
      MAX(I)=MR
      i = i + 1
      GOTO 912
913   CONTINUE
      write(*,*)
      WRITE(*,*) ' Overlap matrix after rotation'
      write(*,*)
      DO 812 K = 1,NUMBERW
         DO 813 I = 1,K
           if (l(k).ne.l(i)) go to 813
            write(*,1000) el(k),el(i),qquadr(k,i,0)
813      CONTINUE
812   CONTINUE
      write(*,*)
      WRITE(*,*) ' The rotated orbitals are written on wfn.rot !'
1000  FORMAT(1H ,'<',A3,'|',A3,'> = ',F18.10)
1001  FORMAT(1H ,' <r> of original ',A3,' electron = ',F18.10
     :,' <|> = ',F18.10)
1002  FORMAT(//,1H ,'The first orbital of this symmetry has # ',I3)
1004  FORMAT(//1H ,'The orbital ',A3,' are  mixed (identical weights)
     : with the ',I4,//,' orbitals of the same symmetry. The latter
     : will be orthogonalized to the inner functions.',//)
      end
*
*     ------------------------------------------------------------------
*              Q U A D R
*     ------------------------------------------------------------------
*
*                                   kk
*       Evaluates the integral of  r   P (r) P (r) with respect to r
*                                       i     j
*
      DOUBLE PRECISION  FUNCTION QQUADR(I,J,KK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NOD=220,NWD=60)
      COMMON /PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,MASS,NCFG,IB,IC,ID
     :   ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,FINE,NSCF,NCLOSD,RMASS
      COMMON /RADIAL/R(NOD),RR(NOD),R2(NOD),P(NOD,NWD),YK(NOD),
     :   YR(NOD),X(NOD),AZ(NWD),L(NWD),MAX(NWD),N(NWD)
      COMMON/NEL/IQP,IQN,IQL,IQAZ,IQMAX,IQ(7)
      DIMENSION G(NOD)
*
      K = KK + 2
      LI = L(I)
      LJ = L(J)
      DEN = LI + LJ + 1 + K
      ZR = Z*R(4)
      BI = (P(4,I)/(AZ(I)*R2(4)*R(4)**LI) - D1+ZR/(LI+1) )/ZR**2
      BJ = (P(4,J)/(AZ(J)*R2(4)*R(4)**LJ) - D1+ZR/(LJ+1) )/ZR**2
      ALPHA= (D1/(LI + 1) + D1/(LJ + 1))/(DEN + D1)
      ZR = Z*R(1)
      BETA = (DEN+D1)*ALPHA**2 - D2*(BI+BJ+D1/((LI+1)*(LJ+1)))/(DEN+D2)
      D = P(1,I)*P(1,J)*R(1)**K*(((BETA*ZR+ALPHA)*ZR+D1)/(DEN*H1)+D5)
      DD = D0
      M = MIN0(MAX(I),MAX(J))
      DO 1 JJ = 1,M
        G(JJ) = P(JJ,I)*P(JJ,J)*R(JJ)**K
    1 CONTINUE
      DO 2 JJ= 3,M,2
        D = D + G(JJ)
    2 CONTINUE
      DO 3 JJ = 2,M,2
        DD = DD + G(JJ)
    3 CONTINUE
      QQUADR = H1*(D + D2*DD)
      RETURN
      END


text:
      SUBROUTINE DPSE4
C
C     PRESSURE STIFFNESS CALCULATIONS FOR A QUADRILATERAL MEMBRANE
C     ELEMENT, WHICH HAS 4 GRID POINTS.
C     THREE 6X6 STIFFNESS MATRICES FOR THE PIVOT POINT ARE INSERTED.
C
C     DOUBLE PRECISION VERSION
C
C     WRITTEN BY E. R. CHRISTENSEN/SVERDRUP,  9/91,  VERSION 1.1
C     INSTALLED IN NASTRAN AS ELEMENT DPSE4 BY G.CHAN/UNISYS, 2/92
C
C     REFERENCE - E. CHRISTENEN: 'ADVACED SOLID ROCKET MOTOR (ASRM)
C                 MATH MODELS - PRESSURE STIFFNESS EFFECTS ANALYSIS',
C                 NASA TD 612-001-02, AUGUST 1991
C
C     LIMITATION -
C     (1) ALL GRID POINTS USED BY ANY OF THE CPSE2/3/4 ELEMENTS MUST BE
C         IN BASIC COORDINATE SYSTEM!!!
C     (2) CONSTANT PRESSURE APPLIED OVER AN ENCLOSED VOLUMN ENCOMPASSED
C         BY THE CPSE2/3/4 ELEMENTRS
C     (3) PRESSURE ACTS NORMALLY TO THE CPSE2/3/4 SURFACES
C
C     SEE NASTRAN DEMONSTRATION PROBLEM -  T13022A
C
      DOUBLE PRECISION GAMMA,KIJ,DP,C,SIGN
      DIMENSION        NECPT(6)
C     COMMON /SYSTEM/  IBUF,NOUT
      COMMON /DS1AAA/  NPVT,ICSTM,NCSTM
      COMMON /DS1AET/  ECPT(26),DUM2(2),DUM12(12)
      COMMON /DS1ADP/  GAMMA,KIJ(36),DP(26),C(12),SIGN(3),NK(3),IK(3)
      EQUIVALENCE      (NECPT(1),ECPT(1))
C
C     ECPT FOR THE PRESSURE STIFFNESS CPES4 ELEMENT
C
C     ECPT( 1) = ELEMENT ID
C     ECPT( 2) = SIL FOR GRID POINT A OR 1
C     ECPT( 3) = SIL FOR GRID POINT B OR 2
C     ECPT( 4) = SIL FOR GRID POINT C OR 3
C     ECPT( 5) = SIL FOR GRID POINT C OR 4
C     ECPT( 6) = PRESSURE
C     ECPT( 7) = NOT USED
C     ECPT( 8) = NOT USED
C     ECPT( 9) = NOT USED
C     ECPT(10) = COORD. SYSTEM ID 1
C     ECPT(11) = X1
C     ECPT(12) = Y1
C     ECPT(13) = Z1
C     ECPT(14) = COORD. SYSTEM ID 2
C     ECPT(15) = X2
C     ECPT(16) = Y2
C     ECPT(17) = Z2
C     ECPT(18) = COORD. SYSTEM ID 3
C     ECPT(19) = X3
C     ECPT(20) = Y3
C     ECPT(21) = Z3
C     ECPT(22) = COORD. SYSTEM ID 4
C     ECPT(23) = X4
C     ECPT(24) = Y4
C     ECPT(25) = Z4
C     ECPT(26) = ELEMENT TEMPERATURE
C     ECPT(27) THRU ECPT(40) = DUM2 AND DUM12, NOT USED IN THIS ROUTINE
C
C     STORE ECPT IN DOUBLE PRECISION
C
      DP(6) = ECPT(6)
      K = 10
      DO 20 I = 1,4
      DO 10 J = 1,3
      K = K + 1
   10 DP(K) = ECPT(K)
   20 K = K + 1
C
C     CALCULATE THE FOUR VECTORS GAB, GAC, GAD, AND GBD USED IN
C     COMPUTING THE PRESSURE STIFFNESS MATRIC
C
C     GAB = RA + RB - RC - RD
C     GAC = RB - RD
C     GAD =-RA + RB + RC - RD
C     GBD =-RA + RC
C
C     GAB STORED IN C( 1), C( 2), C( 3)
C     GAC STORED IN C( 4), C( 5), C( 6)
C     GAD STORED IN C( 7), C( 8), C( 9)
C     GBD STORED IN C(10), C(11), C(12)
C
      C(1) = DP(11) + DP(15) - DP(19) - DP(23)
      C(2) = DP(12) + DP(16) - DP(20) - DP(24)
      C(3) = DP(13) + DP(17) - DP(21) - DP(25)
C
      C(4) = DP(15) - DP(23)
      C(5) = DP(16) - DP(24)
      C(6) = DP(17) - DP(25)
C
      C(7) =-DP(11) + DP(15) + DP(19) - DP(23)
      C(8) =-DP(12) + DP(16) + DP(20) - DP(24)
      C(9) =-DP(13) + DP(17) + DP(21) - DP(25)
C
      C(10)=-DP(11) + DP(19)
      C(11)=-DP(12) + DP(20)
      C(12)=-DP(13) + DP(21)
C
      DO 30 I = 1,4
      IF (NECPT(I+1) .NE. NPVT) GO TO 30
      NPIVOT = I
      GO TO 40
   30 CONTINUE
      RETURN
C
C     GENERATE THE THREE BY THREE PARTITIONS IN GLOBAL COORDINATES HERE
C
C     SET COUNTERS ACCORDING TO WHICH GRID POINT IS THE PIVOT
C
   40 IF (NPIVOT .EQ. 4) GO TO 80
      IF (NPIVOT-2) 50,60,70
C
C     SET COUNTERS AND POINTERS FOR CALCULATING KAB, KAC, KAD
C
   50 NK(1) = 2
      NK(2) = 3
      NK(3) = 4
      IK(1) = 1
      IK(2) = 4
      IK(3) = 7
      SIGN(1) = 1.0D0
      SIGN(2) = 1.0D0
      SIGN(3) = 1.0D0
      GO TO 90
C
C     SET COUNTERS AND POINTERS FOR CALCULATING KBA, KBC, KBD
C     NOTE THAT KBA = -KAB
C
   60 NK(1) = 1
      NK(2) = 3
      NK(3) = 4
      IK(1) = 1
      IK(2) = 7
      IK(3) = 10
      SIGN(1) =-1.0D0
      SIGN(2) = 1.0D0
      SIGN(3) = 1.0D0
      GO TO 90
C
C     SET COUNTERS AND POINTERS FOR CALCULATING KCA, KCB, KCD
C     NOTE THAT KCA = -KAC, KCB = -KBC
C
   70 NK(1) = 1
      NK(2) = 2
      NK(3) = 4
      IK(1) = 4
      IK(2) = 7
      IK(3) = 1
      SIGN(1) =-1.0D0
      SIGN(2) =-1.0D0
      SIGN(3) =-1.0D0
      GO TO 90
C
   80 NK(1) = 1
      NK(2) = 2
      NK(3) = 3
      IK(1) = 7
      IK(2) = 10
      IK(3) = 1
      SIGN(1) =-1.0D0
      SIGN(2) =-1.0D0
      SIGN(3) = 1.0D0
C
   90 GAMMA =-DP(6)/12.0D0
      DO 110 I = 1,3
      DO 100 J = 1,36
  100 KIJ(J) = 0.0D0
      K1 = IK(I)
      K2 = K1 + 1
      K3 = K1 + 2
      SG = GAMMA*SIGN(I)
      KIJ( 2) =-C(K3)*SG
      KIJ( 3) = C(K2)*SG
      KIJ( 7) = C(K3)*SG
      KIJ( 9) =-C(K1)*SG
      KIJ(13) =-C(K2)*SG
      KIJ(14) = C(K1)*SG
C
C     ASSEMBLE INTO THE GLOBAL STIFFNESS MATRIX
C
      IAS = NK(I)
      CALL DS1B (KIJ(1),NECPT(IAS+1))
  110 CONTINUE
C
      RETURN
      END


program matrix_multiplication
    implicit none
    integer :: m, n, p
    integer, allocatable :: a(:, :), b(:, :), c(:, :)
    integer :: i, j, k, summa
    character(len=128) :: tmpstr

    ! Get extents of the matrices.
    read *, m, n, p

    allocate(a(m, n), b(n, p), c(m, p))

    do i = 1, m
        read *, a(i, :)
    end do
    do i = 1, n
        read *, b(i, :)
    end do

    ! 1st way: multiplication in explicit loops.
    first_1: do i = 1, m
        last_1: do j = 1, p
            summa = 0
            middle_1: do k = 1, n
                summa = summa + a(i, k) * b(k, j)
            end do middle_1
            c(i, j) = summa
        end do last_1
    end do first_1

    ! 2nd way: multiplication in two outer loops with a slice in between.
    first_2: do i = 1, m
        last_2: do j = 1, p
            c(i, j) = sum(a(i, :) * b(:, j))
        end do last_2
    end do first_2

    ! 3rd way: multiplication sing intrinsic function matmul().  It applies only
    ! for matrices up to rank 2.
    c = matmul(a, b)

    ! Print the resulting matrix.
    do i = 1, m
        print *, c(i, :)
    end do

    deallocate(a, b, c)
end program matrix_multiplication


text:
!****h* ROBODoc/H5R
!
! NAME
!  MODULE H5R
!
! FILE
!  fortran/src/H5Rff.f90
!
! PURPOSE
!  This file contains Fortran interfaces for H5R functions. It includes
!  all the functions that are independent on whether the Fortran 2003 functions
!  are enabled or disabled.
!
!
! COPYRIGHT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   Copyright by The HDF Group.                                               *
!   Copyright by the Board of Trustees of the University of Illinois.         *
!   All rights reserved.                                                      *
!                                                                             *
!   This file is part of HDF5.  The full HDF5 copyright notice, including     *
!   terms governing use, modification, and redistribution, is contained in    *
!   the COPYING file, which can be found at the root of the source code       *
!   distribution tree, or in https://support.hdfgroup.org/ftp/HDF5/releases.  *
!   If you do not have access to either file, you may request a copy from     *
!   help@hdfgroup.org.                                                        *
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
! NOTES
!                         *** IMPORTANT ***
!  If you add a new H5R function you must add the function name to the
!  Windows dll file 'hdf5_fortrandll.def.in' in the fortran/src directory.
!  This is needed for Windows based operating systems.
!
!*****

MODULE H5R
  USE H5GLOBAL

  ! If you change the value of these parameters, do not forget to change corresponding
  ! values in the H5f90.h file.
  !        INTEGER, PARAMETER :: REF_OBJ_BUF_LEN = 2
  !        INTEGER, PARAMETER :: REF_REG_BUF_LEN = 3
  !
  !        TYPE hobj_ref_t_f
  !             INTEGER ref(REF_OBJ_BUF_LEN)
  !        END TYPE
  !
  !        TYPE hdset_reg_ref_t_f
  !             INTEGER ref(REF_REG_BUF_LEN)
  !        END TYPE
  !

  INTERFACE h5rget_object_type_f

     MODULE PROCEDURE h5rget_object_type_obj_f

  END INTERFACE

CONTAINS

!****s* H5R/h5rget_object_type_obj_f
!
! NAME
!  h5rget_object_type_obj_f
!
! PURPOSE
!  Retrieves the type of object that an object reference points to.
!
! INPUTS
!  dset_id 	 - identifier of the dataset containing
!                  reference to the objects
!  ref 	         - reference to open
! OUTPUTS
!  obj_type 	 - object_type, possible values:
!                    H5G_UNKNOWN_F
!                    H5G_GROUP_F
!                    H5G_DATASET_F
!                    H5G_TYPE_F
!  hdferr 	 - Returns 0 if successful and -1 if fails
!
! AUTHOR
!  Elena Pourmal
!  August 12, 1999
!
! HISTORY
!  Explicit Fortran interfaces were added for
!  called C functions (it is needed for Windows
!  port).  February 28, 2001
!
! NOTES
!  This is a module procedure for the h5rget_object_type_f
!  subroutine.
! SOURCE
  SUBROUTINE h5rget_object_type_obj_f(dset_id, ref, obj_type, hdferr)
    IMPLICIT NONE
    INTEGER(HID_T), INTENT(IN) :: dset_id   ! Dataset identifier
    TYPE(hobj_ref_t_f), INTENT(IN) :: ref   ! Object reference
    INTEGER, INTENT(OUT) :: obj_type     ! Object type
                                         !  H5G_UNKNOWN_F
                                         !  H5G_GROUP_F
                                         !  H5G_DATASET_F
                                         !  H5G_TYPE_F
    INTEGER, INTENT(OUT) :: hdferr       ! Error code
!*****
    INTEGER(HADDR_T) :: ref_f          ! Local buffer to pass reference

    INTERFACE
       INTEGER FUNCTION h5rget_object_type_obj_c(dset_id, ref_f, obj_type)
         USE H5GLOBAL
         !DEC$IF DEFINED(HDF5F90_WINDOWS)
         !DEC$ATTRIBUTES C,reference,decorate,alias:'H5RGET_OBJECT_TYPE_OBJ_C':: h5rget_object_type_obj_c
         !DEC$ENDIF
         !              INTEGER, PARAMETER :: REF_OBJ_BUF_LEN = 2
         INTEGER(HID_T), INTENT(IN) :: dset_id
         INTEGER(HADDR_T) :: ref_f
         INTEGER, INTENT(OUT) :: obj_type
       END FUNCTION h5rget_object_type_obj_c
    END INTERFACE

    ref_f = ref%ref
    hdferr = h5rget_object_type_obj_c(dset_id, ref_f, obj_type )

  END SUBROUTINE h5rget_object_type_obj_f

END MODULE H5R

max_stars_repo_name: GWU-CFD/flash-distro

max_issues_repo_name: GWU-CFD/flash-distro

max_forks_repo_name: GWU-CFD/flash-distro


program functionsDemo

    real :: ans;
    call sub1();
    call sub2(10, 20);

    ans = getsum(10.5, 2.2);
    print *, "sum : ", ans;

end program functionsDemo

subroutine sub1
    print *, "Hello world";
end subroutine sub1

subroutine sub2(a, b)
    integer :: a, b;
    print *, "a : ", a;
    print *, "b : ", b;
end subroutine sub2

function getsum(a, b)
    real :: getsum, a, b;
    getsum = a + b
 end function getsum


text:
\ header.f      - headerless word creation
\ ------------------------------------------------------------------------

 .( loading header.f ) cr

\ ------------------------------------------------------------------------

\ this extension is used to create headerless words.  this allows you to
\ hide a definition from global scope and helps to keep forths name space
\ clean.
\
\ to start creating headerless words you would use the word <headers which
\ is pronounced "from headers".  to switch back into headerfull mode you
\ would use the word headers> (pronounced headers to).  Each of these
\ words can be though of as an arrow pointing towards those words that
\ have headers.
\
\ e.g.
\
\    <headers    \ turn off headers (points back to previous code)
\
\    : foo .... ;
\    : bar .... ;
\
\    headers>    \ turn headers on. points towards new definitions below
\
\    : bam  10 foo 3 bar ;
\
\ once you have reached the end of your module you would execute the word
\ behead.  this word removes all headers for all words defined as being
\ headerless.  the code for headerless words remains but their headers
\ are gone.

\ ------------------------------------------------------------------------

  vocabulary h-voc compiler definitions

\ ------------------------------------------------------------------------
\ make current headerless state headerless!

  0 var h-current         \ real current vocabulary
  0 var h-hp              \ hp true address
  0 var h-last            \ most recent headerfull word
  0 var beheading?        \ true if beheading disabled
  0 var h-state           \ current headerless state

\ 0 = no headerless words defined
\ 1 = headerless words created, headers are on
\ 2 = headerless words created, headers are off

\ ------------------------------------------------------------------------

\ when creating headerless words we store the temporary headers some
\ distance above the current head space pointer address.  we can now
\ switch between headerless and headerfull mode.  all headerfull word
\ headers are written into head space proper.  all headerless headers are
\ written into the headerless pad area.

\ when we have headerless word headers we can only create a certain amount
\ of headerfull words. if we create too many then we can overflow the
\ gap between headers proper and the temporary headerless header pad area.

  8192 const h-pad        \ hp offset to temp store for headerless headers

\ ------------------------------------------------------------------------

  ' h-voc     >body const 'h-voc
  ' h-hp      >body const 'h-hp
  ' h-last    >body const 'h-last
  ' h-current >body const 'h-current
  ' hp        >body const 'hp
  ' last      >body const 'last

\ ------------------------------------------------------------------------
\ enable or disable going headerless

: +heads      ( h-state ?exit)  off> beheading? ;
: -heads      ( h-state ?exit)  on> beheading? ;

\ ------------------------------------------------------------------------
\ swap pointers to real head space and headerless head space etc

: swap-hp       ( --- )  'h-hp   'hp   juggle ;
: swap-last     ( --- )  'h-last 'last juggle ;

\ ------------------------------------------------------------------------
\ state is not headerless but has been before.  go headerless again

: h1        ( --- )
  swap-last                 \ remember most recent headerfull word
  current !> h-current      \ remember true current
  swap-hp                   \ set hp to headerless space.  save real hp
  2 !> h-state              \ all words are created headerless
  h-voc definitions ;       \ adds h-voc to context and current

\ ------------------------------------------------------------------------
\ going headerless for first time (since last behead)

: h0
  off> h-last               \ no previous headerless words yet
  hhere h-pad + !> h-hp     \ point hp 8k beyond where it realy is
  h1 ;                      \ switch to headerless mode

\ ------------------------------------------------------------------------
\ turn headers off (silently ignores switch from headerless to headerless)

: <headers
  beheading? ?exit          \ dont go headerless if beheading disabled
  h-state
  exec: h0 h1 noop ;

\ ------------------------------------------------------------------------
\ turn headers on (silently ignores switch from headerfull to headerfull)

: headers>      ( --- )
  beheading? ?exit          \ dont go headerfull if beheading disabled
  h-state 2 <> ?exit        \ if were headerless go headerfull

  decr> h-state             \ headers are on again now
  swap-hp swap-last
  h-current !> current ;    \ h-voc is still in context though

\ ------------------------------------------------------------------------
\ beheading...
\ ------------------------------------------------------------------------

\ ------------------------------------------------------------------------
\ zero pointers to nfa at cfa -4 for all words in a headerless thread

: (nonames)     ( thread --- )
  begin                     \ for each header in thread do
    @ ?dup
  while
    dup name>               \ go from nfa to cfa
    cell- off               \ erase pointer to nfa at cfa -4
    n>link                  \ point to previous word in thread chain
  repeat ;

\ ------------------------------------------------------------------------
\ zero pointers to nfa at cfa -4 for all words in the h-voc vocabulary

: nonames
  'h-voc #threads
  for
    dup (nonames)           \ blank out nfa pointer at cfa -4 for thread
    cell+                   \ blank this thread out within h-voc
  nxt
  'h-voc 256 erase
  drop ;

\ ------------------------------------------------------------------------
\ erase all headers - gone forever

: behead
  headers>                  \ turn headers on again
  h-voc previous            \ remove h-voc from context
  off> h-state              \ no longer headerless
  off> h-current
  off> h-hp
  off> h-last
  nonames ;                 \ make all beheaded words noname

\ ------------------------------------------------------------------------
\ allows switching current while headerless

  root definitions

: h-defs
  context                   \ point to last item in context stack
  16 #context -
  []@ !> h-current ;        \ copy that to the saved current

\ ------------------------------------------------------------------------
\ wrap the 'definitions' word to account for headerless modes

: definitions               \ create new definition for this word
  h-state                   \ if we are headerless or have been
  if
    h-defs
  else
    definitions             \ this redefinition references the original
  then ;

\ ------------------------------------------------------------------------

  forth definitions

\ ========================================================================


text:
module mod_string
#include <messenger.h>

    type, public :: string_t
        character(len=:), allocatable   :: str
    contains
        procedure   :: lower
        procedure   :: upper
        procedure   :: set
        procedure   :: get
        procedure   :: strip
        procedure   :: replace
!        procedure, private :: write_formatted
!        generic :: write(formatted) => write_formatted
    end type string_t

!    interface string
!        module procedure new_string
!    end interface

    public str
    interface str
        module procedure to_string
    end interface

contains

    !>  Convert data to string representation and return string.
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   11/20/2019
    !!
    !-------------------------------------------------------------------------------------------------
    function to_string(data_in) result(string_out)
        class(*),    intent(in)  :: data_in

        character(100)              :: write_internal
        character(:),   allocatable :: string_out

        string_out = ''

        select type(data_in)

            type is(character(len=*))
                string_out = data_in

            type is(integer)

                if ( data_in == 0 ) then
                    write_internal = '0'

                ! Positive integers
                else if ( data_in < 1e1 .and. data_in > 0) then
                    write(write_internal, '(I1.0)') data_in
                else if ( data_in < 1e2 .and. data_in > 0) then
                    write(write_internal, '(I2.0)') data_in
                else if ( data_in < 1e3 .and. data_in > 0) then
                    write(write_internal, '(I3.0)') data_in
                else if ( data_in < 1e4 .and. data_in > 0) then
                    write(write_internal, '(I4.0)') data_in
                else if ( data_in < 1e5 .and. data_in > 0) then
                    write(write_internal, '(I5.0)') data_in
                else if ( data_in < 1e6 .and. data_in > 0) then
                    write(write_internal, '(I6.0)') data_in
                else if ( data_in < 1e7 .and. data_in > 0) then
                    write(write_internal, '(I7.0)') data_in
                else if ( data_in < 1e8 .and. data_in > 0) then
                    write(write_internal, '(I8.0)') data_in
                else if ( data_in < 1e9 .and. data_in > 0) then
                    write(write_internal, '(I9.0)') data_in
                else if ( data_in < 1e10 .and. data_in > 0) then
                    write(write_internal, '(I10.0)') data_in
                else if ( data_in < 1e11 .and. data_in > 0) then
                    write(write_internal, '(I11.0)') data_in

                ! Negative integers: extra char width for neg sign
                else if ( abs(data_in) < 1e1 .and. data_in < 0) then
                    write(write_internal, '(I2.0)') data_in
                else if ( abs(data_in) < 1e2 .and. data_in < 0) then
                    write(write_internal, '(I3.0)') data_in
                else if ( abs(data_in) < 1e3 .and. data_in < 0) then
                    write(write_internal, '(I4.0)') data_in
                else if ( abs(data_in) < 1e4 .and. data_in < 0) then
                    write(write_internal, '(I5.0)') data_in
                else if ( abs(data_in) < 1e5 .and. data_in < 0) then
                    write(write_internal, '(I6.0)') data_in
                else if ( abs(data_in) < 1e6 .and. data_in < 0) then
                    write(write_internal, '(I7.0)') data_in
                else if ( abs(data_in) < 1e7 .and. data_in < 0) then
                    write(write_internal, '(I8.0)') data_in
                else if ( abs(data_in) < 1e8 .and. data_in < 0) then
                    write(write_internal, '(I9.0)') data_in
                else if ( abs(data_in) < 1e9 .and. data_in < 0) then
                    write(write_internal, '(I10.0)') data_in
                else if ( abs(data_in) < 1e10 .and. data_in < 0) then
                    write(write_internal, '(I11.0)') data_in
                else if ( abs(data_in) < 1e11 .and. data_in < 0) then
                    write(write_internal, '(I12.0)') data_in
                else
                    write(write_internal, '(I10.0)') data_in
                end if
                string_out = write_internal

            type is(integer(8))

                if ( data_in == 0 ) then
                    write_internal = '0'

                ! Positive integers
                else if ( data_in < 1e1 .and. data_in > 0) then
                    write(write_internal, '(I1.0)') data_in
                else if ( data_in < 1e2 .and. data_in > 0) then
                    write(write_internal, '(I2.0)') data_in
                else if ( data_in < 1e3 .and. data_in > 0) then
                    write(write_internal, '(I3.0)') data_in
                else if ( data_in < 1e4 .and. data_in > 0) then
                    write(write_internal, '(I4.0)') data_in
                else if ( data_in < 1e5 .and. data_in > 0) then
                    write(write_internal, '(I5.0)') data_in
                else if ( data_in < 1e6 .and. data_in > 0) then
                    write(write_internal, '(I6.0)') data_in
                else if ( data_in < 1e7 .and. data_in > 0) then
                    write(write_internal, '(I7.0)') data_in
                else if ( data_in < 1e8 .and. data_in > 0) then
                    write(write_internal, '(I8.0)') data_in
                else if ( data_in < 1e9 .and. data_in > 0) then
                    write(write_internal, '(I9.0)') data_in
                else if ( data_in < 1e10 .and. data_in > 0) then
                    write(write_internal, '(I10.0)') data_in
                else if ( data_in < 1e11 .and. data_in > 0) then
                    write(write_internal, '(I11.0)') data_in

                ! Negative integers: extra char width for neg sign
                else if ( abs(data_in) < 1e1 .and. data_in < 0) then
                    write(write_internal, '(I2.0)') data_in
                else if ( abs(data_in) < 1e2 .and. data_in < 0) then
                    write(write_internal, '(I3.0)') data_in
                else if ( abs(data_in) < 1e3 .and. data_in < 0) then
                    write(write_internal, '(I4.0)') data_in
                else if ( abs(data_in) < 1e4 .and. data_in < 0) then
                    write(write_internal, '(I5.0)') data_in
                else if ( abs(data_in) < 1e5 .and. data_in < 0) then
                    write(write_internal, '(I6.0)') data_in
                else if ( abs(data_in) < 1e6 .and. data_in < 0) then
                    write(write_internal, '(I7.0)') data_in
                else if ( abs(data_in) < 1e7 .and. data_in < 0) then
                    write(write_internal, '(I8.0)') data_in
                else if ( abs(data_in) < 1e8 .and. data_in < 0) then
                    write(write_internal, '(I9.0)') data_in
                else if ( abs(data_in) < 1e9 .and. data_in < 0) then
                    write(write_internal, '(I10.0)') data_in
                else if ( abs(data_in) < 1e10 .and. data_in < 0) then
                    write(write_internal, '(I11.0)') data_in
                else if ( abs(data_in) < 1e11 .and. data_in < 0) then
                    write(write_internal, '(I12.0)') data_in
                else
                    write(write_internal, '(I10.0)') data_in
                end if
                string_out = write_internal

            type is(real)
                if (abs(data_in) > 0.1) then
                    write(write_internal, '(F24.14)') data_in
                else
                    write(write_internal, '(E24.14)') data_in
                end if
                string_out = write_internal

            type is(real(8))
                if (abs(data_in) < 0.1) then
                    write(write_internal, '(E24.14)') data_in
                else if ( (abs(data_in) > 0.1) .and. (abs(data_in) < 1.e9) ) then
                    write(write_internal, '(F24.14)') data_in
                else
                    write(write_internal, '(E24.14)') data_in
                end if
                string_out = write_internal

            type is(logical(1))
                if (data_in) then
                    string_out = 'True'
                else
                    string_out = 'False'
                end if

            type is (logical(2))
                if (data_in) then
                    string_out = 'True'
                else
                    string_out = 'False'
                end if

            type is (logical(4))
                if (data_in) then
                    string_out = 'True'
                else
                    string_out = 'False'
                end if

            class default
                print*, 'Error: no IO rule for provided data in add_to_line'
                stop
        end select

    end function to_string
    !***********************************************************************************************

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/30/2016
    !!
    !---------------------------------------------------------------------------------------------
    subroutine set(self,str)
        class(string_t),    intent(inout)   :: self
        character(len=*),   intent(in)      :: str

        self%str = str

    end subroutine set
    !*********************************************************************************************

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/30/2016
    !!
    !---------------------------------------------------------------------------------------------
    function get(self) result(str)
        class(string_t),    intent(in)  :: self

        character(len=:),   allocatable :: str

        str = self%str

    end function get
    !*********************************************************************************************

    !>
    !!
    !!
    !---------------------------------------------------------------------------------------------
    function lower(self) result(str_lower)
        class(string_t),    intent(in)  :: self

        character(len=:), allocatable   :: str_lower

        str_lower = string_to_lower(self%str)

    end function lower
    !*********************************************************************************************

    !>
    !!
    !!
    !---------------------------------------------------------------------------------------------
    function upper(self) result(str_upper)
        class(string_t),    intent(in)  :: self

        character(len=:), allocatable   :: str_upper

        str_upper = string_to_upper(self%str)

    end function upper
    !*********************************************************************************************

    !> It reads in a file name and returns only the prefix, based on the expected extension
    !!
    !!  @author Matteo Ugolotti
    !!  @date   2/27/2017
    !!
    !!  @param[in]  file        Character string containing a file name
    !!  @param[in]  extensions  Character string of the expected extension
    !!
    !---------------------------------------------------------------------------------------------
    function get_file_prefix(file,extension) result(prefix)
        character(*),       intent(in)  :: file
        character(*),       intent(in)  :: extension

        character(len=:), allocatable   :: prefix
        integer                         :: iext, extloc

        ! Check if the file name has already the extension
        extloc = index(file, trim(extension))

        if ( extloc == 0 ) then
            prefix = file
        else
            prefix = trim(file(1:extloc-1))
        end if

    end function get_file_prefix
    !**********************************************************************************************

    !> Given a file name and a list of possible extensions, return an extension if
    !! it is contained in the file string.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/8/2016
    !!
    !!  @param[in]  file        Character string containing a file name
    !!  @param[in]  extensions  Array of character strings. Each string is an accepted extension
    !!
    !---------------------------------------------------------------------------------------------
    function get_file_extension(file,extensions) result(extension)
        character(*),       intent(in)  :: file
        character(len=5),   intent(in)  :: extensions(:)

        character(len=:), allocatable   :: extension
        integer                         :: iext, extloc

        do iext = 1,size(extensions)

            ! Check for extension in grid file
            extloc = index(file, trim(extensions(iext)))

            ! If extloc is nonzero, then the extension was found in the filename
            if ( extloc /= 0 ) then
                extension = trim(extensions(iext))
                exit
            end if

        end do

    end function get_file_extension
    !**********************************************************************************************

    !> Given a file name and an extension, return true if the extention is found.
    !!
    !!  @author Matteo Ugolotti
    !!  @date   10/16/2019
    !!
    !!  @param[in]  file        Character string containing a file name
    !!  @param[in]  extensions  Character string of the extension
    !!
    !---------------------------------------------------------------------------------------------
    function check_file_extension(file,extension) result(extension_found)
        character(*),       intent(in)  :: file
        character(*),       intent(in)  :: extension

        integer     :: iext, extloc
        logical     :: extension_found

        ! Check for extension in grid file
        extloc = index(file, trim(extension))

        ! If extloc is nonzero, then the extension was found in the filename
        if ( extloc /= 0 ) then
            extension_found = .true.
        else
            extension_found = .false.
        end if

    end function check_file_extension
    !**********************************************************************************************

    !>  Converts a character string to all lower-case by comparing the index of each character in the
    !!  ASCII set and shifting the index to the appropriate case of the character.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/8/2016
    !!
    !!
    !-----------------------------------------------------------------------------------------------
    function string_to_lower(string) result(string_lower)
        character(*),   intent(in)  :: string

        character(len=len(string))  :: string_lower
        integer                         :: i, k

        do i = 1,len(trim(string))

            k = iachar(string(i:i))

            if ( k >= iachar('A') .and. k <= iachar('Z') ) then
                k = k + (iachar('a') - iachar('A'))
                string_lower(i:i) = achar(k)
            else
                string_lower(i:i) = string(i:i)
            end if

        end do ! i
    end function string_to_lower
    !************************************************************************************************

    !>  Converts a character string to all upper-case by comparing the index of each character in the
    !!  ASCII set and shifting the index to the appropriate case of the character.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/8/2016
    !!
    !!
    !-----------------------------------------------------------------------------------------------
    function string_to_upper(string) result(string_upper)
        character(*),   intent(in)  :: string

        character(len=len(string))  :: string_upper
        integer                         :: i, k

        do i = 1,len(trim(string))

            k = iachar(string(i:i))

            if ( k >= iachar('a') .and. k <= iachar('z') ) then
                k = k - (iachar('a') - iachar('A'))
                string_upper(i:i) = achar(k)
            else
                string_upper(i:i) = string(i:i)
            end if

        end do ! i
    end function string_to_upper
    !************************************************************************************************

    !>  Split the string at the given separator and return an array of strings
    !!
    !!  @author Matteo Ugolotti
    !!  @date   2/12/2017
    !!
    !-----------------------------------------------------------------------------------------------
    function strip(self,separator) result(string_list)
        class(string_t),    intent(in)  :: self
        character(*),       intent(in)  :: separator

        type(string_t), allocatable     :: string_list(:)
        character(len=100)              :: temp_list(100)
        integer                         :: ierr, i, index_l, pos
        logical                         :: searching
        character(len=:), allocatable   :: usr_msg_1, usr_msg_2, word, string

        searching = .true.
        index_l   = 0
        string = self%get()

        ! Split the string in words and add each word to temp_list
        do while (searching)

            index_l = index_l + 1
            pos = index(string,separator)

            if (pos == 0 ) then
                word = string
                searching = .false.
            else
                word = string(1:pos-1)
            end if

            ! Check if the number of words is greater than the size of temp_list
            usr_msg_1 = "string%strip: The number of words found in the stripped string "
            usr_msg_2 = "are too many."
            if (index_l > size(temp_list)) call chidg_signal_two(FATAL,usr_msg_1,self%str,usr_msg_2)

            ! Add word to temp_list
            temp_list(index_l) = word
            string = string(pos+1:)

        end do

        ! Allocate the vector of strings based on the word found
        allocate(string_list(index_l), stat=ierr)
        if (ierr /= 0) call chidg_signal(FATAL, "strip: error in allocating the string_list")

        ! Add words to the output vector
        do i = 1,index_l
            ! Adjust and trim the word and add it to string_list
            string_list(i)%str = trim(adjustl(temp_list(i)))
        end do

    end function strip
    !************************************************************************************************

    !>  Replace a charcter (before) with a new character (after)
    !!
    !!  @author Matteo Ugolotti
    !!  @date   2/12/2017
    !!
    !-----------------------------------------------------------------------------------------------
    function replace(self,before,after) result(string_out)
        class(string_t),    intent(in)  :: self
        character(*),       intent(in)  :: before
        character(*),       intent(in)  :: after

        type(string_t)                  :: string_out
        integer                         :: ierr, i, index_l, pos, after_len, before_len
        logical                         :: searching
        character(len=:), allocatable   :: word, string

        searching = .true.
        string = self%get()

        after_len  = len(after)
        before_len = len(before)
        i = 1

        ! Split the string in words and add each word to temp_list
        do while (searching)

            pos = index(string,before)

            if (pos == 0 ) then
                word = string
                searching = .false.
            else
                word = string(:pos-1)// after // string(pos+before_len:)
            end if

            string = word

            if (i == 100) then
                searching = .false.
            end if

            i = i + 1

        end do

        ! Adjust and trim the word and add it to string_list
        call string_out%set(trim(adjustl(word)))

    end function replace
    !************************************************************************************************

end module mod_string

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text: Music Scene Local rock band featuring Users/JoeFinkel on guitar, Brett Hart on guitar, and Users/TeddyBriggs on bosstone.