       program main
!      *************************************************************************
!      ** determines the permutation vectors defining the topologically       **
!      ** inequivalent, linked feynman diagrams of the specified order        **
!      **                                                                     **
!      ** provide the order of the diagram from standard input (integer)      **
!      **                                                                     **
!      ** interaction lines connect odd numbered vertices with the next higher**
!      ** even numbered index. the j-th component p_j of the permutation      **
!      ** specifies a green's function from vertex j to vertex p_j.           **
!      **                                                                     **
!      ** the value of a diagram is given in                                  **
!      **   p.Bloechl, PhiSX Advanced Solid-State Theory                      **
!      **   http://www2.pt.tu-clausthal.de/atp/phisx.html                     **
!      **                                                                     **
!      *************************************************************************
!      ** COPYRIGHT (C) 2017 PETER BLOECHL, (MODIFIED M.I.T. LICENSE)         **
!      **                                                                     **
!      ** PERMISSION IS HEREBY GRANTED, FREE OF CHARGE, TO ANY PERSON         **
!      ** OBTAINING A COPY OF THIS SOFTWARE (THE "SOFTWARE"), TO DEAL IN THE  **
!      ** SOFTWARE WITHOUT RESTRICTION, INCLUDING WITHOUT LIMITATION THE      **
!      ** RIGHTS TO USE, COPY, MODIFY, MERGE, PUBLISH, DISTRIBUTE, SUBLICENSE,**
!      ** AND/OR SELL COPIES OF THE SOFTWARE, AND TO PERMIT PERSONS TO WHOM   **
!      ** THE SOFTWARE IS FURNISHED TO DO SO, SUBJECT TO THE FOLLOWING        **
!      ** CONDITIONS:                                                         **
!      **                                                                     **
!      ** THE ABOVE COPYRIGHT NOTICE AND THIS PERMISSION NOTICE SHALL BE      **
!      ** INCLUDED IN ALL COPIES OR SUBSTANTIAL PORTIONS OF THE SOFTWARE.     **
!      ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,     **
!      ** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF  **
!      ** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND               **
!      ** NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS **
!      ** BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN  **
!      ** ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN   **
!      ** CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE    **
!      ** SOFTWARE.                                                           **
!      **                                                                     **
!      ** THIS LICENSE APPLIES ONLY TO THE SECTION OF FORTRAN CODE THAT       **
!      ** CONTAINS THIS LICENSE. IT DOES NOT APPLY TO THE PHISX LECTURE NOTES **
!      ** THAT MAY CONTAIN THE SOFTWARE.                                      **
!      *************************************************************************
       implicit none
       integer(4)             :: order     !order in the interaction
       integer(4)             :: len       !length of permutation vector
       integer(4)             :: nvecx     !X#(permutation vectors)
       integer(4)             :: nmapx
       integer(4),allocatable :: pvec(:,:) !(len,nvecx) permutation vector
       integer(4),allocatable :: map(:,:)  !(len,nmapx)
       integer(4),allocatable :: sofd(:)   !(nvecx) symmetry factor*sign
       integer(4),allocatable :: sgn(:)    !(nvecx)
       integer(4)             :: nvec      !#(permutation vectors)
       integer(4)             :: nmap
       integer(4)             :: i,j
       real(8)                :: svar
!      *************************************************************************
       write(*,fmt='("provide the order of the diagrams (integer):")')
       read(*,fmt=*)order
       write(*,fmt='(40("."),t1,"order of the diagrams",t40,i15)')order
       len=2*order
!
!      =========================================================================
!      == allocate array for permutation vector                               ==
!      == a permutation vector has the length len=2*order                     ==
!      == the number of permutation vectors in this order is factorial(len)   ==
!      =========================================================================
       nvecx=1
       do i=1,len
         nvecx=nvecx*i
       enddo
       allocate(pvec(len,nvecx))
       allocate(sgn(nvecx))
       allocate(sofd(nvecx))

!      =========================================================================
!      == determine all permutation vectors                                   ==
!      =========================================================================
       call permutationvectors(len,nvecx,nvec,pvec)
       write(*,fmt='(40("."),t1,"nr. of all diagrams (=n!)",t40,i15)')nvec
!
!      =========================================================================
!      == count the number of loops (sign of the diagram)                     ==
!      == identify unlinked diagrams (sgn=0)                                  ==
!      =========================================================================
       do i=1,nvec
         call loops(len,pvec(:,i),sgn(i))
       enddo
!
!      =========================================================================
!      == remove unlinked diagrams                                            ==
!      =========================================================================
       j=0
       do i=1,nvec
         if(sgn(i).eq.0) cycle
         j=j+1
         pvec(:,j)=pvec(:,i)
         sgn(j)=sgn(i)
!         write(*,fmt='("(",20i4,")",i10)')pvec(:,j),sgn(j)
       enddo
       nvec=j
       write(*,fmt='(40("."),t1,"nr. of linked diagrams",t40,i15)')nvec
!
!      =========================================================================
!      == determine symmetry factor (sofd=0 for duplicates)                   ==
!      == sign is integrated into the symmetry factor                         ==
!      =========================================================================
       nmapx=1
       do i=1,order
         nmapx=nmapx*2*i
       enddo
       allocate(map(len,nmapx))
       call vertexmaps(len,nmapx,nmap,map)
       do i=1,nvec
         call symmetryfactor(len,pvec(:,i),nmap,map,sofd(i))
         sofd(i)=sofd(i)*sgn(i)
       enddo
       deallocate(map)
!
!      =========================================================================
!      == remove topologically equivalent diagrams                            ==
!      =========================================================================
       j=0
       do i=1,nvec
         if(sofd(i).eq.0) cycle
         j=j+1
         pvec(:,j)=pvec(:,i)
         sofd(j)=sofd(i)
         sgn(j)=sgn(i)
       enddo
       nvec=j
       write(*,fmt='(40("."),t1,"nr. of top. distinct diagrams",t40,i15)')nvec
!
!      =========================================================================
!      == test sum rule                                                       ==
!      =========================================================================
       svar=sum(1.d0/abs(sofd(:nvec)))
       do i=1,len/2
         svar=svar*2.d0*real(i,kind=8)
       enddo
       write(*,fmt='(40("."),t1,"sum-rule test 2^n * n! * S(D) ",t40,f15.1)')svar
       write(*,fmt='(40("."),t1,"...must reproduce number of linked diagrams")')
!
!      =========================================================================
!      == print resulting permutation vectors and symmetry factors            ==
!      =========================================================================
       write(*,*)
       write(*,fmt='("sign, symmetry factor and permutation vectors:")')
       do i=1,nvec
         write(*,fmt='(i8,"  sign*sofd=",i8," p=",20i4)') &
      &         i,sofd(i),pvec(:,i)
       enddo
       stop
       end
!
!      ..1.........2.........3.........4.........5.........6.........7.........8
       subroutine vertexmaps(len,nmapx,nmap,map)
!      *************************************************************************
!      ** determine the mappings of vertices that produce topologically       **
!      ** equivalent diagrams.                                                **
!      ** the map array contains all permutations of the vertices among each  **
!      ** other.
!      *************************************************************************
       implicit none
       integer(4),intent(in) :: len    ! #(green's functions)
       integer(4),intent(in) :: nmapx
       integer(4),intent(out):: nmap
       integer(4),intent(out):: map(len,nmapx)
       integer(4),allocatable:: pvec(:,:)
       integer(4)            :: nvecx
       integer(4)            :: nvec
       integer(4)            :: nmap2,imap,jmap,ivec,i,j
!      *************************************************************************
       if(mod(len,2).eq.1) then
         stop 'argument len must be even'
       end if
!
!      == allocate dimension of map array ======================================
       nvecx=1
       do i=1,len/2
         nvecx=nvecx*i
       enddo
       allocate(pvec(len/2,nvecx))
       call permutationvectors(len/2,nvecx,nvec,pvec)
!
!      == construct all permutation of vertices
       nmap=0
       do ivec=1,nvec
         nmap=nmap+1
         do i=1,len/2
           map(2*i-1,nmap)=2*pvec(i,ivec)-1
           map(2*i  ,nmap)=2*pvec(i,ivec)
         enddo
         nmap2=nmap
         jmap=nmap
         do j=1,len/2
           do imap=nmap,nmap2
             jmap=jmap+1
             map(:,jmap)=map(:,imap)
             map(2*j-1,jmap)=map(2*j,imap)
             map(2*j,jmap)=map(2*j-1,imap)
           enddo
           nmap2=jmap
         enddo
         nmap=nmap2
       enddo
       return
       end
!
!      ..1.........2.........3.........4.........5.........6.........7.........8
       subroutine symmetryfactor(len,pvec,nmap,map,sofd)
!      *************************************************************************
!      ** determine symmetry factor                                           **
!      ** a zero value for the symmetry factor is set if a topologically      **
!      ** equivalent diagram has already been considered.                     **
!      **                                                                     **
!      ** The symmetry factor Sofd is the ratio of mere deformations among    **
!      ** all n!2^n diagrams obtained by remapping the vertices.              **
!      **                                                                     **
!      ** for each diagram there are n!2^n diagrams that give the same value. **
!      ** These diagrams are obtained by remapping the vertices. These        **
!      ** these mappings (1) permutate the interaction lines and              **
!      ** (2) they interchange the two vertices of one interaction.           **
!      **                                                                     **
!      ** A deformation is characterized by the same permutation vector and   **
!      ** is thus the identical diagram. Hence it must not be counted.        **
!      **                                                                     **
!      ** each diagram is converted into a value obtained from the permutation**
!      ** vector. Only those diagrams are retained, for which the             **
!      ** topologically equivalent partners have a higher number than the     **
!      ** diagram itself. This singles out one diagram from each class of     **
!      ** topologically distinct diagrams.                                    **
!      **                                                                     **
!      *************************************************************************
       implicit none
       integer(4),intent(in) :: len         ! #(green's function in the diagram)
       integer(4),intent(in) :: pvec(len)   ! permutation vector
       integer(4),intent(in) :: nmap
       integer(4),intent(in) :: map(len,nmap)
       integer(4),intent(out):: sofd        ! symmetry factor
       integer(4)            :: imap,i
       integer(4)            :: ndef        ! #(deformations)
       integer(4)            :: ovec(len)
       logical(4)            :: tident
       integer(4)            :: order
       integer(4)            :: ndistinct
!      *************************************************************************
       order=len/2
!
!      =========================================================================
!      == initial checks                                                      ==
!      =========================================================================
!      == check that len is even. It is the number of green's functions and   ==
!      == that number is twice the number of interaction lines                ==
       if(mod(len,2).eq.1) then
         stop 'argument len must be even'
       end if
!
!      =========================================================================
!      == loop over all mappings 
!      =========================================================================
       do imap=1,nmap
!        == remap vertices =====================================================
         do i=1,len
           ovec(map(i,imap))=map(pvec(i),imap)
         enddo
!
!        =======================================================================
!        == exclude all diagrams that map onto another one with a smaller count=
!        =======================================================================
         tident=.true.
         do i=1,len
           if(ovec(i).ne.pvec(i)) then
             if(tident.and.ovec(i).lt.pvec(i)) then
               sofd=0 ! is already counted by another pvec
               return
             end if 
             exit
           end if
         enddo
       enddo

!      =========================================================================
!      == count mere deformations                                             ==
!      =========================================================================
       ndef=0
       do imap=1,nmap
         do i=1,len
           ovec(map(i,imap))=map(pvec(i),imap)
         enddo
         tident=.true.
         do i=1,len
           if(ovec(i).ne.pvec(i)) then
             tident=.false.
             exit
           end if
         enddo
         if(tident)ndef=ndef+1
       enddo
!
!      =========================================================================
!      == calculate symmetry factor                                           ==
!      =========================================================================
!      == the diagram has ndef deformations. the number of topologically      ==
!      == equivalent diagrams which are not mere deformations is              ==
!      == ndistinct=nmap/ndef.                                                ==
!      == the inverse of symmetry factor is the number of                     ==
!      == topologically equivalent, which are not mere deformations,          ==
!      == divided by the factor n!2^n.                                        ==
!      == Hence, 1/sofd=ndistinct/nmap=1/ndef                                 ==
!      =========================================================================
       ndistinct=nmap/ndef
       sofd=ndef !=nmap/ndistinct
       return
       end
!
!      ..1.........2.........3.........4.........5.........6.........7.........8
       subroutine loops(len,pvec,sgn)
!      *************************************************************************
!      **  sgn=1 for an even number of green's function loops                 **
!      **  sgn=-1 for an odd number of green's function loops                 **
!      **  sgn=0 for a diagram which is not a linked cluster                  **
!      **                                                                     **
!      ** the sign theorem links the sign of a diagram to the number of       **
!      ** loops formed by Green's functions.                                  **
!      **                                                                     **
!      ** the code follows green's function lines and marks each vertex       **
!      ** visited by a loop number. If an un-marked vertex (iloop=0) is       **
!      ** encountered, the code searches for interaction lines, of which one  **
!      ** vertex has been visited and the other not. At such an interaction   **
!      ** the loop index in increased, and loop search  at the unmarked vertex**
!      ** is initiated.                                                       **
!      **                                                                     **
!      ** if interaction lines with two unmarked vertices are left over       **
!      ** the diagram consists of disconnected parts and sgn=0 is set.        **
!      ** otherwise the sign is determined from the the number of loops.      **
!      **                                                                     **
!      *************************************************************************
       implicit none
       integer(4),intent(in) :: len        ! #(greens functions)
       integer(4),intent(in) :: pvec(len)  ! permutation vector of diagram
       integer(4),intent(out):: sgn        ! sign or (zero for linked cluster)
       integer(4)            :: work(len)  ! holds loop indices for each vertex
       logical(4)            :: t1,t2
       integer(4)            :: iloop      ! loop index
       integer(4)            :: ip         ! vertex pointer
       integer(4)            :: i
!      *************************************************************************
       if(mod(len,2).eq.1) then
         write(*,*)'argument len must be even'
         stop 'in loops'
       end if
!
!      =========================================================================
!      == explore diagram for fermi loops                                     ==
!      =========================================================================
       work(:)=0
       iloop=1
       ip=1
       work(ip)=iloop
!
1000   continue
!
!      == walk around the first loop of greensfunctions and mark with loop nr.
       do 
         ip=pvec(ip)
         if(work(ip).ne.0) exit ! end of Fermi loop encountered
         work(ip)=iloop
       enddo
!
!      == find interaction with one vertex that has been visited and another 
!      == one that has not been visited.
       ip=0
       do i=1,len/2
         t1=work(2*i-1).eq.0
         t2=work(2*i).eq.0
         if(t1.eqv.t2) cycle  ! both are true or both are false.
!        -- set pointer to unmarked vertex -----------------------------------
         if(t1) then 
           ip=2*i-1
         else
           ip=2*i
         end if
       enddo
!      == jump back to next fermi-loop search or exit
       if(ip.ne.0) then  ! vertex for fermi-loop search found
         iloop=iloop+1   ! increase loop index
         work(ip)=iloop  ! mark actual vertex
         goto 1000
       end if
!
!      =========================================================================
!      == calculate sgn                                                      ==
!      =========================================================================
       sgn=1
       if(mod(iloop,2).eq.1) sgn=-1
!
!      =========================================================================
!      == check if this is a linked diagram                                   ==
!      =========================================================================
       do i=1,len
         if(work(i).eq.0) then
           sgn=0
           exit
         end if
       enddo
!
       return
       end
!
!      ..1.........2.........3.........4.........5.........6.........7.........8
       subroutine permutationvectors(len,nvecx,nvec,pvec)
!      *************************************************************************
!      ** construct all len-factorial permutation vectors of length len       **
!      ** as described in section 11.1.1 "Construct all permutation vectors"  **
!      ** of P.Bloechl PhiSX Advanced Solid-State Theory                      **
!      *************************************************************************
       implicit none
       integer(4),intent(in)  :: len
       integer(4),intent(in)  :: nvecx           ! X#(permutation vectors)
       integer(4),intent(out) :: nvec            ! #(permutation vectors)
       integer(4),intent(out) :: pvec(len,nvecx) ! permutation vectors
       integer(4)             :: i,j,k,ip1,ip2,isvar
       integer(4)             :: work(len)
!      *************************************************************************
       isvar=1
       do i=1,len
         isvar=isvar*i
       enddo
       if(nvecx.lt.isvar) then
         stop 'nvecx too small (a)'
       end if
!
!      =========================================================================
!      ==  start the recursive construction with trivial vector (1,2,.., len) ==
!      =========================================================================
       nvec=1
       do i=1,len
         pvec(i,nvec)=i
       enddo       
!
!      =========================================================================
!      == construct all other permutation vectors recusively
!      =========================================================================
       do
!
!        == find longest trailing descending sequence ==========================
!        == that is pvec(ip1) > pvec(ip1+1) > pvec(ip1+2) > ... ================
         j=len
         ip1=0
         do j=len,2,-1
           if(pvec(j-1,nvec).gt.pvec(j,nvec)) cycle
           ip1=j-1
           exit
         enddo
         if(ip1.eq.0) exit ! sequence is descending, exit
!
!        == prepare an entry for a new perturbation vector =====================
         nvec=nvec+1
         if(nvec.gt.nvecx) then
           stop 'nvec exceeds nvecx (b)'
         end if
!
!        == copy previous permutation vector as non-descending part ============
         pvec(:,nvec)=pvec(:,nvec-1)

!        -- ip1 is the start of the descending series. -------------------------        
!        -- identify the smallest component of pvec in the descending sequence -
!        -- that is still larger than the first element of the descending ------
!        -- sequence. ----------------------------------------------------------
         ip2=ip1+1
         do k=ip1+1,len
           if(pvec(k,nvec).lt.pvec(ip1,nvec))exit
           ip2=k
         end do
!
!        -- replace first element of descending seqence ------------------------
         isvar=pvec(ip1,nvec)
         pvec(ip1,nvec)=pvec(ip2,nvec)
         pvec(ip2,nvec)=isvar
!
!        -- change descending sequence into an ascending one -------------------
         work(ip1+1:)=pvec(ip1+1:,nvec)
         do j=1,len-ip1
           pvec(ip1+j,nvec)=work(len+1-j)
         enddo
       enddo
       return
       end
