!! DENSITY MATRIX SUBROUTINES !!

#include "Definitions.INC"

!! ALL MODULES

module denutilmod
contains

!! ASSUMES denmat,denvects is ready.

subroutine output_denmat( incmf, time)
  use parameters
  use mpimod
  use xxxmod
  implicit none
  integer :: i,j, incmf, numcalledhere=0, mytimingout
  real*8 ::  rsum, rsum2, time, rsum3

  if (myrank.ne.1) then
     return
  endif
  numcalledhere=numcalledhere+1

  if (incmf==1) then
     mytimingout=1
  else
     mytimingout=timingout
  endif
  if (numcalledhere==1) then
     open(853,file=dendatfile,status="unknown")
     write(853,*) ;           close(853)
     open(853,file=denrotfile,status="unknown")
     write(853,*) ;           close(853)
  endif

  if (mod(numcalledhere,mytimingout).eq.0)  then
     rsum=0.d0;  rsum2=0.d0;  rsum3=0.d0
     do i=1,nspf
        do j=1,nspf
           if (i /= j) then
              rsum=rsum+abs(yyy%denmat(i,j,0)**2)
           else
              rsum2=rsum2+abs(yyy%denmat(i,j,0)**2)
           endif
           if (shells(i) /= shells(j)) then
              rsum3=rsum3+abs(yyy%denmat(i,j,0)**2)
           endif
        enddo
     enddo

     open(853,file=dendatfile,status="unknown", position="append")
     write(853,'(F17.8, 100E17.7)') time, yyy%denvals(1:nspf); close(853)

     open(853,file=denrotfile,status="unknown", position="append")
     write(853,'(100F23.14)') time, rsum/rsum2, rsum3/rsum2;     close(853)

!!$     if (cdenflag.ne.0.and.rdenflag.ne.0) then
!!$        !!     call openfile(); write(mpifileptr,*) "Calling schmidt check"; call closefile()
!!$        call schmidtcheck(1)
!!$     endif
  endif
end subroutine output_denmat

subroutine denprint(iprop)
  use parameters
  use xxxmod
  implicit none
  integer :: i,iprop
  CNORMTYPE ::  sum2

  sum2=0d0
  do i=1,nspf
     sum2=sum2+yyy%denvals(i)
  enddo
  if (sum2.eq.0.d0) then
     print *, "SUM2 ERR!!"
  endif
  if (iprop.eq.1) then
     OFL; WRFL "|-------Natural orbital occupation numbers for each state-------------------|"
  endif
  write(mpifileptr,'(100F9.5)') (abs(yyy%denvals(i)/sum2*numelec),i=1,nspf); CFL

end subroutine denprint

end module denutilmod

module natrepbiomod
  use biorthotypemod
  implicit none
  type(biorthotype),target :: natrepbiovar
end module

module repnatmod
contains

subroutine replace_withnat(printflag)
  use natrepbiomod
  use biorthomod
  use parameters
  use configmod    !! bioww
  use xxxmod
  use spfsubmod
  implicit none

  DATATYPE,allocatable :: outspfs(:,:)
  integer :: i,j,printflag,imc
  real*8 :: errorval
  DATATYPE,target :: smo(nspf,nspf)

  allocate(outspfs(spfsize,nspf))
  outspfs=0d0

  do j=1,nspf  ! which natorb
     do i=1,nspf  ! which original
        outspfs(:,j)=outspfs(:,j)+ &
             yyy%cmfspfs((i-1)*spfsize+1:i*spfsize,0)*yyy%denvects(i,j)
     enddo
  enddo

  call spf_orthogit(outspfs, errorval)
  if (errorval.gt.1d-7) then
     OFLWR "ERROR IN REPLACENAT ", errorval; CFLST
  endif

  if (printflag==1) then
     OFLWR "REPLACING SPFS - DENMAT EIGS"
     do i=1,nspf
        write(mpifileptr,'(2E25.10)') yyy%denvals(i)/mcscfnum
     enddo
     WRFL; CFL
  endif

  call bioset(natrepbiovar,smo,numr,bioww)

  do imc=1,mcscfnum

!!$     tempavector(:)=yyy%cmfavec(:,imc,0)

     call biotransform(yyy%cmfspfs(:,0),outspfs, yyy%cmfavec(:,imc,0),natrepbiovar)

!!$! doesn't work.  permoverlaps broken presumably.
!!$#ifdef NOWAYDUDE
!!$     call autocorrelate_one(yyy%cmfavec(:,imc,0), outspfs, yyy%cmfspfs(:,0), tempavector, sum,1,needbiovar)
!!$!! CHECK
!!$     tempavector2(:)=yyy%cmfavec(:,imc,0)
!!$     call biotransform(outspfs,yyy%cmfspfs(:,0),tempavector2(:),natrepbiovar)
!!$     print *, DOT_PRODUCT(tempavector,tempavector),"AAAAAA"
!!$     tempavector2(:)=     tempavector2(:) - tempavector(:)
!!$     print *, DOT_PRODUCT(tempavector2,tempavector2),"BBBBB";
!!$     nullspfs(:,:)=0d0
!!$nullspfs(:,:)=nullspfs(:,:)-outspfs(:,:)
!!$print *, "BIOCHECK ", DOT_PRODUCT(RESHAPE(nullspfs, (/spfsize*nspf/)), RESHAPE(nullspfs, (/spfsize*nspf/)))
!!$!     call checkbio(yyy%cmfspfs(:,0),outspfs,tempavector,yyy%cmfavec(:,imc,0))
!!$!     call checkbio(yyy%cmfspfs(:,0),nullspfs,tempavector,yyy%cmfavec(:,imc,0))
!!$#endif

  enddo

  yyy%cmfspfs(:,0)=RESHAPE(outspfs,(/totspfdim/))

  deallocate(outspfs)

end subroutine replace_withnat

end module

!! denmat is the true denmat, not transposed.

module densubmod
contains

subroutine getdenmat00(inholeflag,wwin,avector1,in_avector2,rvector, denmat, numpoints,howmany)
  use walkmod
  use dotmod
  use mpisubmod
  implicit none
  type(walktype),intent(in) :: wwin
  integer,intent(in) ::  numpoints,howmany,inholeflag
  DATATYPE, intent(in) :: in_avector2(numpoints,wwin%firstconfig:wwin%lastconfig,howmany),&
       avector1(numpoints,wwin%firstconfig:wwin%lastconfig,howmany)
  DATAECS,intent(in) :: rvector(numpoints)
  DATATYPE,intent(out) :: denmat(wwin%nspf,wwin%nspf)
  DATATYPE :: a1(numpoints,howmany), a2(numpoints,howmany), &
       mydenmat(wwin%nspf,wwin%nspf), csum
  DATATYPE, allocatable :: avector2(:,:,:)
  integer :: config1,config2,  ispf,jspf,  dirphase, iwalk, ii, ihop

  allocate(avector2(numpoints,wwin%numconfig,howmany))
  avector2(:,:,:)=0d0
  if (wwin%lastconfig.ge.wwin%firstconfig) then
     avector2(:,wwin%firstconfig:wwin%lastconfig,:) = in_avector2(:,:,:)
  endif

!! DO SUMMA (parconsplit.ne.0 and sparsesummaflag.eq.2, "circ")

  if (wwin%parconsplit.ne.0) then
     do ii=1,howmany
        call mpiallgather(avector2(:,:,ii),wwin%numconfig*numpoints,&
             wwin%configsperproc(:)*numpoints,wwin%maxconfigsperproc*numpoints)
     enddo
  endif

  denmat(:,:)=0.d0

!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(config1,ii,iwalk,config2,dirphase,ispf,jspf,mydenmat,a1,a2,csum,ihop)

!! removed  REDUCTION(+:denmat) after problems in get_twoereduced, OMP CRITICAL instead

  mydenmat(:,:)=0.d0

!$OMP DO SCHEDULE(DYNAMIC)
  do config1=wwin%botconfig,wwin%topconfig

     do ii=1,howmany
        a1(:,ii)=avector1(:,config1,ii) * rvector(:)
     enddo

     do ihop=1,wwin%numsinglehops(config1)
        config2=wwin%singlehop(ihop,config1)

        a2(:,:)=avector2(:,config2,:)
        csum=dot(a2(:,:),a1(:,:),numpoints*howmany)

        do iwalk=wwin%singlehopwalkstart(ihop,config1),wwin%singlehopwalkend(ihop,config1)

           dirphase=wwin%singlewalkdirphase(iwalk+wwin%scol(config1))
           ispf=wwin%singlewalkopspf(1,iwalk+wwin%scol(config1))  !! goes with config1
           jspf=wwin%singlewalkopspf(2,iwalk+wwin%scol(config1))  !! goes with config2

           mydenmat(ispf,jspf)=mydenmat(ispf,jspf)+ &
                csum*dirphase

!!$!!  ONCE AND FOR ALL 2014           :                     RIGHT MULTIPLYING!!!
!!$           denmat(ispf,jspf)=denmat(ispf,jspf)+ &
!!$                dirphase*a1*CONJUGATE(a2)
!!$!! NOT
!!$!!           denmat(jspf,ispf)=denmat(jspf,ispf)+ &
!!$!!                dirphase*a1*CONJUGATE(a2)

        enddo

     enddo
  enddo
!$OMP END DO
!$OMP CRITICAL
  denmat(:,:)=denmat(:,:)+mydenmat(:,:)
!$OMP END CRITICAL
!$OMP END PARALLEL

  call mympireduce(denmat,wwin%nspf**2)

  if (wwin%holeflag.ne.0.and.inholeflag.ne.0) then
     csum=0d0
     if (wwin%lastconfig.ge.wwin%firstconfig) then
        csum=dot(in_avector2,avector1,numpoints*wwin%localnconfig*howmany)
     endif
     if (wwin%parconsplit.ne.0) then
        call mympireduceone(csum)
     endif
     denmat(:,:) = (-1) * denmat(:,:)
     do ispf=1,wwin%nspf
        denmat(ispf,ispf)=denmat(ispf,ispf) + csum * 2d0
     enddo
  endif

  deallocate(avector2)

end subroutine getdenmat00

end module densubmod

module denmatxmod
contains

subroutine getdenmatx()
  use parameters
  use configmod
  use xxxmod
  implicit none

  call getdenmatstuff(www,yyy%cmfavec(:,:,0), yyy%denmat(:,:,0) , &
       yyy%invdenmat(:,:,0) , yyy%denvals(:) , yyy%denvects(:,:), numr, mcscfnum)

!!  OFLWR "DENCHECK"
!!#ifdef REALGO
!!  write(mpifileptr,'(F20.10)') yyy%denmat(:,1,0)
!!#else
!!  write(mpifileptr,'(2F20.10)') yyy%denmat(:,1,0)
!!#endif
!!  WRFL "DENSTOP!"; CFLST

!!$  if (rdenflag==1) then
!!$     call getrdenmat()
!!$  endif
!!$  if (cdenflag==1) then
!!$     call getnatconfig()
!!$  endif
!!$  if (cdenflag==1.and.rdenflag==1) then
!!$     call schmidtcheck(0)
!!$  endif

contains
subroutine getdenmatstuff(wwin,avector, denmat, invdenmat, denvals, &
     denvects, numpoints,howmany)
  use class_parameters
  use denreg_parameters
  use walkmod
  use invsubmod
  use densubmod
  use eigenmod
  implicit none
  type(walktype),intent(in) :: wwin
  integer,intent(in) ::  numpoints,howmany
  CNORMTYPE,intent(out) :: denvals(wwin%nspf)
  DATATYPE, intent(in) :: avector(numpoints,wwin%firstconfig:wwin%lastconfig,howmany)
  DATATYPE,intent(out) :: denmat(wwin%nspf,wwin%nspf),&
       invdenmat(wwin%nspf,wwin%nspf),denvects(wwin%nspf,wwin%nspf)
  CNORMTYPE :: tempdenvals(wwin%nspf,numclasses)
  DATATYPE :: tempinvden(wwin%nspf,wwin%nspf,numclasses),&
       tempdenvects(wwin%nspf,wwin%nspf,numclasses), csum
  integer :: ispf,jspf,iclass
  DATAECS :: rvector(numpoints)

  rvector(:)=1d0
  call getdenmat00(1,wwin,avector,avector,rvector,denmat,numpoints,howmany)

  denvects(:,:)=0d0; invdenmat(:,:)=0d0

  tempinvden(:,:,:)=0d0; tempdenvals(:,:)=0d0; tempdenvects(:,:,:)=0d0;

  do iclass=1,numclasses
     do ispf=1,nperclass(iclass)
        do jspf=1,nperclass(iclass)
           tempinvden(ispf,jspf,iclass) = &
                denmat(classorb(ispf,iclass),classorb(jspf,iclass))
        enddo
     enddo
     tempinvden(:,:,iclass)= (-1) * tempinvden(:,:,iclass)
     call EIGEN(tempinvden(:,:,iclass),nperclass(iclass), wwin%nspf, &
          tempdenvects(:,:,iclass),tempdenvals(:,iclass))
     tempinvden(:,:,iclass)= (-1) * tempinvden(:,:,iclass)
     tempdenvals(:,iclass)= (-1) * tempdenvals(:,iclass)

     call invmatsmooth(tempinvden(:,:,iclass),nperclass(iclass),wwin%nspf,denreg,.false.)
     do ispf=1,nperclass(iclass)
        denvals(classorb(ispf,iclass))=tempdenvals(ispf,iclass)
        do jspf=1,nperclass(iclass)
           denvects(classorb(ispf,iclass),classorb(jspf,iclass)) = &
                tempdenvects(ispf,jspf,iclass)
           invdenmat(classorb(ispf,iclass),classorb(jspf,iclass)) = &
                tempinvden(ispf,jspf,iclass)
        enddo
     enddo
  enddo

  do ispf=1,wwin%nspf
     csum=denvects(ispf,ispf)
#ifdef CNORMFLAG
     if (real(csum,8).lt.0d0) then
        denvects(:,ispf) = denvects(:,ispf) * (-1)
     endif
#else
     if (csum.ne.0d0) then
        denvects(:,ispf) = denvects(:,ispf) * abs(csum) / csum
     endif
#endif
  enddo

end subroutine getdenmatstuff

end subroutine getdenmatx

end module denmatxmod

module getoccmod
contains

subroutine getoccupations(wwin,in_avector, numpoints, occupations)
  use walkmod
  use dotmod
  use mpisubmod
  implicit none
  type(walktype),intent(in) :: wwin
  integer,intent(in) ::  numpoints
  CNORMTYPE, intent(out) :: occupations(wwin%nspf)
  integer :: config1,  ispf,jspf,  iwalk,idiag
  DATATYPE, intent(in) :: in_avector(numpoints,wwin%firstconfig:wwin%lastconfig)
  DATATYPE, allocatable :: avector(:,:)

  allocate(avector(numpoints,wwin%numconfig))
  avector(:,:)=0d0

  if (wwin%lastconfig.ge.wwin%firstconfig) then
     avector(:,wwin%firstconfig:wwin%lastconfig) = in_avector(:,:)
  endif

!! DO SUMMA (parconsplit.ne.0 and sparsesummaflag.eq.2, "circ")

  if (wwin%parconsplit.ne.0) then
     call mpiallgather(avector,wwin%numconfig*numpoints,&
          wwin%configsperproc(:)*numpoints,wwin%maxconfigsperproc*numpoints)
  endif

  occupations(:)=0d0

  do config1=wwin%botconfig,wwin%topconfig
     do idiag=1,wwin%numsinglediagwalks(config1)
        iwalk=wwin%singlediag(idiag,config1)

        ispf=wwin%singlewalkopspf(1,iwalk+wwin%scol(config1))  !! goes with config1
        jspf=wwin%singlewalkopspf(2,iwalk+wwin%scol(config1))  !! goes with config2

        occupations(ispf)=occupations(ispf) + &  !! ok conversion
             dot(avector(:,config1),avector(:,config1),numpoints)

     enddo
  enddo

  deallocate(avector)

#ifndef REALGO
#ifndef CNORMFLAG
  call mympirealreduce(occupations,wwin%nspf)
#else
  call mympireduce(occupations,wwin%nspf)
#endif
#else
  call mympireduce(occupations,wwin%nspf)
#endif

  if (wwin%holeflag.ne.0) then
     occupations(:)=2d0-occupations(:)
  endif

end subroutine getoccupations

end module getoccmod

!!$subroutine getrdenmat()
!!$  use parameters
!!$  use xxxmod
!!$  implicit none
!!$
!!$  DATATYPE :: avector(numr,numconfig,mcscfnum), a1(mcscfnum), a2(mcscfnum), &
!!$       csum, csum2, sum,dot
!!$  integer ::     ii,jj, config1
!!$
!!$  avector(:,:,:)=RESHAPE(yyy%cmfavec(:,:,0),(/numr,numconfig,mcscfnum/))
!!$  do ii=1,numr
!!$     do jj=1,numr
!!$        sum=0.d0
!!$        do config1=1,numconfig  !! ok walks
!!$           a1(:)=avector(ii,config1,:);
!!$           a2(:)=avector(jj,config1,:)
!!$           sum=sum+dot(a2(:),a1(:),mcscfnum)
!!$        enddo
!!$        yyy%rdenmat(ii,jj)=sum
!!$     enddo
!!$  enddo
!!$
!!$  yyy%rdenmat=(-1)*yyy%rdenmat
!!$  call EIGEN(yyy%rdenmat,numr, numr, yyy%rdenvects(:,:), yyy%rdenvals(:))
!!$  yyy%rdenvals(:)=(-1)*yyy%rdenvals(:) ;  yyy%rdenmat=(-1)*yyy%rdenmat
!!$
!!$  do ii=1,numr
!!$     csum=0.d0;     csum2=0.d0
!!$     do jj=1,numr
!!$
!!$!! 111510 this is just for output, no big deal, but rdenexpect should be
!!$        CNORMTYPE, really DATATYPE if want complex expect for chmctdh
!!$
!!$        csum=csum+              ALLCON(yyy%rdenvects(jj,ii))*     yyy%rdenvects(jj,ii)
!!$        csum2=csum2+ ALLCON(yyy%rdenvects(jj,ii))*yyy%rdenvects(jj,ii)*bondpoints(ii)
!!$     enddo
!!$     yyy%rdenexpect(ii)=csum2/csum   !! ok for implicit conv (chmctdh, pmctdh, mctdh)
!!$  enddo
!!$
!!$end subroutine getrdenmat

!!$  KEEPME KEEPME KEEMPE
!!$
!!$!! checks whether projections of natural configurations are
!!$!! proportional to natural orbitals in r
!!$!! only valid for mcscfnum=1....   projections are for each state,
!!$!! r natorbs are averaged
!!$
!!$subroutine schmidtcheck(schmidtprint)
!!$  use parameters
!!$  use xxxmod
!!$  use natprojmod
!!$  implicit none
!!$
!!$  DATATYPE :: dot
!!$  real*8 :: rsum
!!$  integer ::  i,j,ssize, schmidtprint
!!$  DATATYPE, allocatable :: schmidt(:,:)
!!$
!!$  if (numr.eq.1.or.mcscfnum.ne.1) then
!!$     return
!!$  endif
!!$  do i=1,numr
!!$     if (abs(yyy%rdenvals(i)).lt.1.d-8) then
!!$        exit
!!$     endif
!!$  enddo
!!$  ssize=i-1;  allocate(schmidt(ssize,ssize))
!!$  do i=1,ssize
!!$     do j=1,ssize
!!$        schmidt(i,j)=dot(natproj(:,i,1),yyy%rdenvects(:,j),numr)
!!$     enddo
!!$  enddo
!!$  if (schmidtprint.ne.0) then
!!$     call openfile()
!!$     write(mpifileptr,*) "------- Schmidt Overlaps ----------"
!!$     do i=1,ssize
!!$        write(mpifileptr,'(100E12.4)') schmidt(:,i)
!!$     enddo
!!$     write(mpifileptr,*) "-----------------------------------"
!!$  endif
!!$  do i=1,ssize
!!$     do j=1,ssize
!!$        rsum=abs(schmidt(i,j))**2
!!$        if (i==j) then
!!$#ifndef CNORMFLAG
!!$           rsum=rsum-yyy%rdenvals(i)
!!$        endif
!!$#else
!!$        else
!!$#endif
!!$           if (abs(rsum).gt.1.d-12) then
!!$              OFLWR "WARNING!  Schmidt dot fail", rsum, i,j,schmidt(i,j),yyy%rdenvals(i);CFL
!!$           endif
!!$#ifdef CNORMFLAG
!!$        endif
!!$#endif
!!$     enddo
!!$  enddo
!!$
!!$  deallocate(schmidt)
!!$end subroutine schmidtcheck
!!$
!!$
!!$subroutine getnatconfig()
!!$  use parameters
!!$  use natprojmod
!!$  use xxxmod
!!$  implicit none
!!$
!!$  DATATYPE :: avector(numr,numconfig,mcscfnum), a1(mcscfnum), a2(mcscfnum),dot
!!$  integer ::  config1,config2,     ii,imc,i,ir,numvects
!!$  DATATYPE :: natmat(numconfig,numconfig), natvects(numconfig,numconfig),tempconfig(numconfig)
!!$
!!$  avector(:,:,:)=RESHAPE(yyy%cmfavec(:,:,0),(/numr,numconfig,numr,mcscfnum/))
!!$
!!$  if (numconfig.gt.2000) then
!!$     OFLWR "Error: getting natconfigs but numconfig gt 2000: ", numconfig; CFLST
!!$  endif
!!$  natmat=0.d0
!!$  do ii=1,numr
!!$     do config1=1,numconfig    !! ok walks
!!$        a1(:)=avector(ii,config1,:)
!!$        do config2=1,numconfig !! ok walks
!!$           a2(:)=avector(ii,config2,:)
!!$           natmat(config1,config2) = natmat(config1,config2) + dot(a2,a1,mcscfnum)
!!$        enddo
!!$     enddo
!!$  enddo
!!$
!!$  natmat=(-1)*natmat
!!$  call EIGEN(natmat,numconfig, numconfig, natvects, natvals)
!!$
!!$  natmat=(-1)*natmat;  natvals=(-1)*natvals
!!$  numvects=min(numr,numconfig);  natconfigs=0.d0
!!$
!!$  do imc=1,mcscfnum
!!$     do i=1,numvects
!!$        do ir=1,numr
!!$        tempvector(:)=avector(ir,:,imc)
!!$           natproj(ir,i,imc)=dot(natvects(:,i),  & !! ok imp conv (p,ch)
!!$                tempvector(:),numconfig)
!!$        enddo
!!$        natdot(i,imc)=dot(natproj(:,i,imc),       & !! ok imp conv (p,ch)
!!$             natproj(:,i,imc),numr)
!!$     enddo
!!$  enddo
!!$
!!$  natconfigs(:,1:numvects)=natvects(:,1:numvects)
!!$
!!$end subroutine getnatconfig


text:
              Character filename*20,res*3,atom*4,aothers*20,a*80
              character b(30)*80,atot*80,temp*80,a1(30)*80,arg(10)*32
              dimension x(1000000),y(1000000),z(1000000),num(1000000),
     *        iatm(1000000),ires(1000000),res(1000000),atom(1000000),
     *        aothers(10000000)

              do i = 1,iargc()
                call getarg(i, arg(i))
                  write(*,*)i,arg(i)
              enddo
              write(*,*)"input file:",arg(1)
              write(*,*)"outputfile:", arg(2)
              open(unit=10,file=arg(1))
              open(unit=4,file=arg(2))
              do i=1,5
                read(10,'(a80)')a
                write(*,'(a80)')a
                write(4,'(a80)')a
              enddo
              irescopy=1
              read(10,'(6x,i5)')iresnum
              if (iresnum.eq.1) then
               ii=1
              else
               ii=0
              endif
              rewind(unit=10)
              do i=1,5
                read(10,'(a80)')a
              enddo
              iii=1
              do i=1,100000
               read(10,1,END=99)iatm(i),ires(i),res(i),atom(i),
     *         x(i),y(i),z(i),aothers(i)
c              if (ires(1).eq.1)  then
c                    ii=1
c              else
c                    ii=0
c              endif
               if (ires(i).ne.irescopy)  then
                ii=ii+1
               endif
               if (ii.eq.1) then
                 if (atom(i).eq."HT1 ") then
                 write(*,*)"check HT1"
                    iii=1
                    iatm(i)=iii
                     write(3,2)iatm(i),ii,res(i),atom(i),
     *               x(i),y(i),z(i),aothers(i)
                 elseif (atom(i).eq."N   ") then
                 write(*,*)"check N"
                    iii=2
                    iatm(i)=iii
                     write(3,2)iatm(i),ii,res(i),atom(i),
     *               x(i),y(i),z(i),aothers(i)
                 elseif (atom(i).eq."HT2 ") then
                 write(*,*)"check HT2"
                    iii=3
                    iatm(i)=iii
                     write(3,2)iatm(i),ii,res(i),atom(i),
     *               x(i),y(i),z(i),aothers(i)
                 elseif (atom(i).eq."HT3 ") then
                 write(*,*)"check HT3"
                    iii=4
                    iatm(i)=iii
                     write(3,2)iatm(i),ii,res(i),atom(i),
     *               x(i),y(i),z(i),aothers(i)
                 elseif (atom(i).eq."CA  ") then
                 write(*,*)"check CA"
                    iii=5
                    iatm(i)=iii
                     write(3,2)iatm(i),ii,res(i),atom(i),
     *               x(i),y(i),z(i),aothers(i)
                 else
                     iii=iii+1
                    iatm(i)=iii
                     write(3,2)iatm(i),ii,res(i),atom(i),
     *               x(i),y(i),z(i),aothers(i)
                    ilist=iii
                 endif
                 write(*,*)"iii",iii
               else
                     write(3,2)i,ii,res(i),atom(i),x(i),y(i),z(i),
     *               aothers(i)
               endif
               irescopy=ires(i)
              enddo
99            continue
              iserial=i
                 write(*,*)"ilist",ilist
              close(unit=3)
              open (unit=3,file='fort.3')
              do i5=1,ilist
               read(3,'(a80)')a1(i5)
               write(*,'(a80)')a1(i5)
                read(a1(i5)(:5),'(i5)')num(i5)
              enddo
               DO 15 I=1,ilist
                AMIN=num(i)
                temp=a1(i)
                DO 25 K=i+1,ilist
                 IF (AMIN.GE.num(K)) then
                  Zx=AMIN
                  AMIN=num(K)
                  num(K)=Zx
                  atot=temp
                  temp=a1(k)
                  a1(k)=atot
                 endif
25              CONTINUE
                num(I5)=AMIN
                a1(i5)=temp
                write(4,'(a80)')a1(i5)
15             continue
               write(*,*)"ilist, iserial-1",ilist,iserial-1
               read (3,2) (iatm(i5),ires(i5),res(i5),atom(i5),
     *         x(i5),y(i5),z(i5),aothers(i5),I5=ilist+1,iserial-1)
               WRITE (4,2) (iatm(i5),ires(i5),res(i5),atom(i5),
     *         x(i5),y(i5),z(i5),aothers(i5),I5=ilist+1,iserial-1)
1             format(i5,3x,i3,1x,a3,2x,a4,f8.3,3x,f8.3,2x,f8.3,a20)
2             format(i5,1x,i4,1x,a3,2x,a4,3f10.5,a20)
              stop
              end

ext: f

max_stars_repo_name: vam-sin/deepcys

max_issues_repo_name: vam-sin/deepcys

max_forks_repo_name: vam-sin/deepcys


text:
C--------------------------------------------------------------------------
        subroutine ftgbtp(ibuff,keynam,value,status)

C       Get Binary Table Parameter
C       test if the keyword is one of the table column definition keywords
C       of a binary table. If so, decode it and update the values in the common
C       block

C       ibuff   i sequence number of the data buffer
C       OUTPUT PARAMETERS:
C       keynam  c name of the keyword
C       value   c value of the keyword
C       status  i  returned error status (0=ok)
C
C       written by Wm Pence, HEASARC/GSFC, June 1991

        integer ibuff,status,width
        character keynam*8,value*70

C-------COMMON BLOCK DEFINITIONS:--------------------------------------------
C       nb = number of file buffers = max. number of FITS file opened at once
C       nf = maximum number of fields allowed in a table
        integer nf,nb
        parameter (nb = 20)
        parameter (nf = 3000)
        integer tfield,tstart,tbcol,rowlen,tdtype,trept,tnull,scount
        integer theap,nxheap
        double precision tscale,tzero
        common/ft0002/tfield(nb),tstart(nb),tbcol(nf),rowlen(nb),
     &  tdtype(nf),trept(nf),tscale(nf),tzero(nf),tnull(nf),scount(nb)
     &  ,theap(nb),nxheap(nb)
C-------END OF COMMON BLOCK DEFINITIONS-----------------------------------

        integer nfield,tstat
        character tform*16

        if (status .gt. 0)return
        tstat=status

        if (keynam(1:5) .eq. 'TFORM')then
C               get the field number
                call ftc2ii(keynam(6:8),nfield,status)
                if (status .gt. 0)then
C                   this must not have been a TFORMn keyword
                    status=tstat
                else
C                   get the TFORM character string, without quotes
                    call ftc2s(value,tform,status)
C                   get the datatype code and repeat count
                    call ftbnfm(tform,tdtype(nfield+tstart(ibuff)),
     &                 trept(nfield+tstart(ibuff)),width,status)
                    if (tdtype(nfield+tstart(ibuff)) .eq. 1)then
C                       treat Bit datatype as if it were a Byte datatype
                        tdtype(nfield+tstart(ibuff))=11
                        trept(nfield+tstart(ibuff))=(trept(nfield+
     &                  tstart(ibuff))+7)/8
                    else if (tdtype(nfield+tstart(ibuff)) .eq. 16)then
C                      store the width of the ASCII field in the TNULL parameter
                        tnull(nfield+tstart(ibuff))=width
                   end if
                end if
        else if (keynam(1:5) .eq. 'TSCAL')then
C               get the field number
                call ftc2ii(keynam(6:8),nfield,status)
                if (status .gt. 0)then
C                   this must not have been a TSCALn keyword
                    status=tstat
                else
C                   get the scale factor
                    call ftc2dd(value,tscale(nfield+tstart(ibuff)),
     &                          status)
                    if (status .gt. 0)then
                         call ftpmsg('Error reading value of'//keynam
     &                //' as a Double: '//value)
                    end if
                end if
        else if (keynam(1:5) .eq. 'TZERO')then
C               get the field number
                call ftc2ii(keynam(6:8),nfield,status)
                if (status .gt. 0)then
C                   this must not have been a TZEROn keyword
                    status=tstat
                else
C                   get the scaling zero point
                    call ftc2dd(value,tzero(nfield+tstart(ibuff)),
     &                          status)
                    if (status .gt. 0)then
                         call ftpmsg('Error reading value of'//keynam
     &                //' as a Double: '//value)
                    end if
                end if
        else if (keynam(1:5) .eq. 'TNULL')then
C               get the field number
                call ftc2ii(keynam(6:8),nfield,status)
                if (status .gt. 0)then
C                   this must not have been a TNULLn keyword
                    status=tstat
                else
C                   make sure this is not an ASCII column (the tnull
C                   variable is use to store the ASCII column width)
                    if (tdtype(nfield+tstart(ibuff)) .ne. 16)then
C                       get the Null value flag (Integer)
                        call ftc2ii(value,tnull(nfield+tstart(ibuff)),
     &                              status)
                        if (status .gt. 0)then
                            call ftpmsg('Error reading value of '//
     &                      keynam//' as an integer: '//value)
                        end if
                    end if
                end if
        else if (keynam(1:8) .eq. 'THEAP   ')then
C               get the heap offset value
                call ftc2ii(value,theap(ibuff),status)
                if (status .gt. 0)then
                        call ftpmsg('Error reading value of '//keynam
     &                  //' as an integer: '//value)
                end if
        end if
        end

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


module GwfGwtExchangeModule

  use KindModule,              only: DP, I4B
  use ListsModule,             only: basemodellist, baseexchangelist
  use SimVariablesModule,      only: errmsg
  use BaseExchangeModule,      only: BaseExchangeType, AddBaseExchangeToList
  use BaseModelModule,         only: BaseModelType, GetBaseModelFromList
  use GwfModule,               only: GwfModelType
  use GwtModule,               only: GwtModelType
  use BndModule,               only: BndType, GetBndFromList

  implicit none
  public :: GwfGwtExchangeType
  public :: gwfgwt_cr

  type, extends(BaseExchangeType) :: GwfGwtExchangeType

    integer(I4B), pointer :: m1id => null()
    integer(I4B), pointer :: m2id => null()

  contains

    procedure :: exg_df
    procedure :: exg_ar
    procedure :: exg_da
    procedure, private :: set_model_pointers
    procedure, private :: allocate_scalars
    procedure, private :: gwfbnd2gwtfmi

  end type GwfGwtExchangeType

  contains

  subroutine gwfgwt_cr(filename, id, m1id, m2id)
! ******************************************************************************
! gwfgwt_cr -- Create a new GWF to GWT exchange object
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    ! -- dummy
    character(len=*), intent(in) :: filename
    integer(I4B), intent(in) :: id
    integer(I4B), intent(in) :: m1id
    integer(I4B), intent(in) :: m2id
    ! -- local
    class(BaseExchangeType), pointer :: baseexchange => null()
    type(GwfGwtExchangeType), pointer :: exchange => null()
    character(len=20) :: cint
! ------------------------------------------------------------------------------
    !
    ! -- Create a new exchange and add it to the baseexchangelist container
    allocate(exchange)
    baseexchange => exchange
    call AddBaseExchangeToList(baseexchangelist, baseexchange)
    !
    ! -- Assign id and name
    exchange%id = id
    write(cint, '(i0)') id
    exchange%name = 'GWF-GWT_' // trim(adjustl(cint))
    exchange%memoryPath = exchange%name
    !
    ! -- allocate scalars
    call exchange%allocate_scalars()
    exchange%m1id = m1id
    exchange%m2id = m2id
    !
    ! -- set model pointers
    call exchange%set_model_pointers()
    !
    ! -- return
    return
  end subroutine gwfgwt_cr

  subroutine set_model_pointers(this)
! ******************************************************************************
! set_model_pointers -- allocate and read
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    ! -- dummy
    class(GwfGwtExchangeType) :: this
    ! -- local
    class(BaseModelType), pointer :: mb => null()
    type(GwfModelType), pointer :: gwfmodel => null()
    type(GwtModelType), pointer :: gwtmodel => null()
    integer(I4B) :: ngwfpack, ip
    class(BndType), pointer :: packobj => null()
! ------------------------------------------------------------------------------
    !
    ! -- set gwfmodel
    mb => GetBaseModelFromList(basemodellist, this%m1id)
    select type (mb)
    type is (GwfModelType)
      gwfmodel => mb
    end select
    !
    ! -- set gwtmodel
    mb => GetBaseModelFromList(basemodellist, this%m2id)
    select type (mb)
    type is (GwtModelType)
      gwtmodel => mb
    end select
    !
    ! -- Tell transport model fmi flows are not read from file
    gwtmodel%fmi%flows_from_file = .false.
    !
    ! -- setup pointer to gwf variables that were allocated in gwf_df
    gwtmodel%fmi%gwfbndlist => gwfmodel%bndlist
    ngwfpack = gwfmodel%bndlist%Count()
    !
    ! -- allocate arrays in fmi of size ngwfpack
    call gwtmodel%fmi%allocate_gwfpackages(ngwfpack)
    !
    ! -- Assign values in the fmi package
    do ip = 1, ngwfpack
      packobj => GetBndFromList(gwfmodel%bndlist, ip)
      call gwtmodel%fmi%gwfpackages(ip)%set_name(packobj%packName)
      gwtmodel%fmi%flowpacknamearray(ip) = packobj%packName
    end do
    !
    ! -- return
    return
  end subroutine set_model_pointers

  subroutine exg_df(this)
! ******************************************************************************
! exg_df -- define
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    ! -- dummy
    class(GwfGwtExchangeType) :: this
    ! -- local
    class(BaseModelType), pointer :: mb => null()
    type(GwfModelType), pointer :: gwfmodel => null()
    type(GwtModelType), pointer :: gwtmodel => null()
    integer(I4B) :: ngwfpack, ip
    class(BndType), pointer :: packobj => null()
! ------------------------------------------------------------------------------
    !
    !
    ! -- set gwfmodel
    mb => GetBaseModelFromList(basemodellist, this%m1id)
    select type (mb)
    type is (GwfModelType)
      gwfmodel => mb
    end select
    !
    ! -- set gwtmodel
    mb => GetBaseModelFromList(basemodellist, this%m2id)
    select type (mb)
    type is (GwtModelType)
      gwtmodel => mb
    end select
    !
    ! -- Set pointer to flowja
    gwtmodel%fmi%gwfflowja => gwfmodel%flowja
    !
    ! -- Set the npf flag so that specific discharge is available for
    !    transport calculations if dispersion is active
    if (gwtmodel%indsp > 0) then
      gwfmodel%npf%icalcspdis = 1
    end if
    !
    ! -- Set the auxiliary names for gwf flow packages in gwt%fmi
    ngwfpack = gwfmodel%bndlist%Count()
    do ip = 1, ngwfpack
      packobj => GetBndFromList(gwfmodel%bndlist, ip)
      call gwtmodel%fmi%gwfpackages(ip)%set_auxname(packobj%naux,              &
                                                    packobj%auxname)
    end do
    !
    ! -- return
    return
  end subroutine exg_df

  subroutine exg_ar(this)
! ******************************************************************************
! exg_ar --
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use SimModule, only: store_error, ustop
    ! -- dummy
    class(GwfGwtExchangeType) :: this
    ! -- local
    class(BaseModelType), pointer :: mb => null()
    type(GwfModelType), pointer :: gwfmodel => null()
    type(GwtModelType), pointer :: gwtmodel => null()
    ! -- formats
    character(len=*),parameter :: fmtdiserr = &
      "('GWF and GWT Models do not have the same discretization for exchange&
      & ',a,'.&
      &  GWF Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
      &  GWT Model has ', i0, ' user nodes and ', i0, ' reduced nodes.&
      &  Ensure discretization packages, including IDOMAIN, are identical.')"
! ------------------------------------------------------------------------------
    !
    ! -- set gwfmodel
    mb => GetBaseModelFromList(basemodellist, this%m1id)
    select type (mb)
    type is (GwfModelType)
      gwfmodel => mb
    end select
    !
    ! -- set gwtmodel
    mb => GetBaseModelFromList(basemodellist, this%m2id)
    select type (mb)
    type is (GwtModelType)
      gwtmodel => mb
    end select
    !
    ! -- Check to make sure sizes are identical
    if (gwtmodel%dis%nodes /= gwfmodel%dis%nodes .or.&
        gwtmodel%dis%nodesuser /= gwfmodel%dis%nodesuser) then
      write(errmsg, fmtdiserr) trim(this%name), &
                               gwfmodel%dis%nodesuser, &
                               gwfmodel%dis%nodes, &
                               gwtmodel%dis%nodesuser, &
                               gwtmodel%dis%nodes
      call store_error(errmsg)
      call ustop()
    end if
    !
    ! -- setup pointers to gwf variables allocated in gwf_ar
    gwtmodel%fmi%gwfhead   => gwfmodel%x
    gwtmodel%fmi%gwfsat    => gwfmodel%npf%sat
    gwtmodel%fmi%gwfspdis  => gwfmodel%npf%spdis
    !
    ! -- setup pointers to the flow storage rates. GWF strg arrays are
    !    available after the gwf_ar routine is called.
    if(gwtmodel%inmst > 0) then
      if (gwfmodel%insto > 0) then
        gwtmodel%fmi%gwfstrgss => gwfmodel%sto%strgss
        gwtmodel%fmi%igwfstrgss = 1
        if (gwfmodel%sto%iusesy == 1) then
          gwtmodel%fmi%gwfstrgsy => gwfmodel%sto%strgsy
          gwtmodel%fmi%igwfstrgsy = 1
        endif
      endif
    endif
    !
    ! -- Set a pointer to conc
    if (gwfmodel%inbuy > 0) then
      call gwfmodel%buy%set_concentration_pointer(gwtmodel%name, gwtmodel%x, &
                                                  gwtmodel%ibound)
    endif
    !
    ! -- transfer the boundary package information from gwf to gwt
    call this%gwfbnd2gwtfmi()
    !
    ! -- if mover package is active, then set a pointer to it's budget object
    if (gwfmodel%inmvr /= 0) then
      gwtmodel%fmi%mvrbudobj => gwfmodel%mvr%budobj
    end if
    !
    ! -- return
    return
  end subroutine exg_ar

  subroutine exg_da(this)
! ******************************************************************************
! allocate_scalars
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MemoryManagerModule, only: mem_deallocate
    ! -- dummy
    class(GwfGwtExchangeType) :: this
    ! -- local
! ------------------------------------------------------------------------------
    !
    call mem_deallocate(this%m1id)
    call mem_deallocate(this%m2id)
    !
    ! -- return
    return
  end subroutine exg_da

  subroutine allocate_scalars(this)
! ******************************************************************************
! allocate_scalars
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MemoryManagerModule, only: mem_allocate
    ! -- dummy
    class(GwfGwtExchangeType) :: this
    ! -- local
! ------------------------------------------------------------------------------
    !
    call mem_allocate(this%m1id, 'M1ID', this%memoryPath)
    call mem_allocate(this%m2id, 'M2ID', this%memoryPath)
    this%m1id = 0
    this%m2id = 0
    !
    ! -- return
    return
  end subroutine allocate_scalars

  subroutine gwfbnd2gwtfmi(this)
! ******************************************************************************
! gwfbnd2gwtfmi
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    ! -- dummy
    class(GwfGwtExchangeType) :: this
    ! -- local
    integer(I4B) :: ngwfpack, ip
    class(BaseModelType), pointer :: mb => null()
    type(GwfModelType), pointer :: gwfmodel => null()
    type(GwtModelType), pointer :: gwtmodel => null()
    class(BndType), pointer :: packobj => null()
! ------------------------------------------------------------------------------
    !
    ! -- set gwfmodel
    mb => GetBaseModelFromList(basemodellist, this%m1id)
    select type (mb)
    type is (GwfModelType)
      gwfmodel => mb
    end select
    !
    ! -- set gwtmodel
    mb => GetBaseModelFromList(basemodellist, this%m2id)
    select type (mb)
    type is (GwtModelType)
      gwtmodel => mb
    end select
    !
    ! -- Allocate the gwfpackages in fmi and transfer information
    ngwfpack = gwfmodel%bndlist%Count()
    do ip = 1, ngwfpack
      packobj => GetBndFromList(gwfmodel%bndlist, ip)
      call gwtmodel%fmi%gwfpackages(ip)%set_pointers( &
                           packobj%packName, &
                           packobj%auxname, &
                           packobj%nbound, &
                           packobj%naux, &
                           packobj%nodelist, &
                           packobj%hcof, &
                           packobj%rhs, &
                           packobj%auxvar, &
                           packobj%xnew)
    end do
    !
    ! -- return
    return
  end subroutine gwfbnd2gwtfmi

end module GwfGwtExchangeModule


text: #redirect Ghausi Hall

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
subroutine opt_Turnon(a)
   use typre
   use Mod_Optics
   implicit none

   class(OpticsProblem) :: a

   call a%Timer%Total%Tic
   call a%Timer%Turnon%Tic

   call a%SetExmod
   call a%OpenFiles
   if (a%MPIrank == a%MPIroot) then

      !Read Physical Problem
      call a%Reaphy

      !Read Numerical Treatment
      call a%Reanut

      !Read Output Strategy
      call a%Reaous

   endif
   call a%ReaMPI

   !Warnings and errors
   call a%Outerr(1)

   !Memory allocation
   call a%Memall

   !Warnings and errors
   call a%Outerr(2)

   call a%CloseDataFile

   call a%Timer%Turnon%Toc
   call a%Timer%Total%Toc

end subroutine opt_Turnon

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
module map_mod
    character, dimension(:,:), allocatable :: game_map
    integer :: map_w, map_h
contains
    ! Collision Check
    logical function check_coll(x, y) result(can_pass)
        implicit none
        integer, value :: x, y
        if (x .lt. 1 .or. x .gt. map_w .or. y .lt. 1 .or. y .gt. map_h) then
            can_pass = .false.
        else
            can_pass = .not. (game_map(x,y) .eq. '#')
        end if
    end function
end module map_mod

max_stars_repo_name: Rezznov/FortranRoguelike

max_issues_repo_name: Rezznov/FortranRoguelike

max_forks_repo_name: Rezznov/FortranRoguelike


PROGRAM Master
  ! This program demonstrates how to use the MPMD (Multiple Programs Multiple
  ! Data) model. In this case, a master/worker MPMD program is demonstrated.

  IMPLICIT NONE
  INCLUDE 'mpif.h'

  ! Configuration
  INTEGER, PARAMETER :: njobmax = 10

  ! MPI variables
  INTEGER :: nprocs, myrank, tag, dummy, err, dest
  INTEGER, DIMENSION(MPI_STATUS_SIZE) :: status

  ! Count variables
  INTEGER job, i

  CALL MPI_INIT(err)
  CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, err)
  CALL MPI_COMM_RANK(MPI_COMM_WORLD, myrank, err)

  ! Distribute the work among the workers
  tag = 1
  DO job = 1, njobmax
    CALL MPI_RECV(dummy, 1, MPI_INTEGER, MPI_ANY_SOURCE, tag, MPI_COMM_WORLD, &
                  status, err)
    dest = status(MPI_SOURCE)
    CALL MPI_SEND(job, 1, MPI_INTEGER, dest, tag, MPI_COMM_WORLD, err)
  ENDDO

  ! Tell the workers to stop
  DO i = 1, nprocs-1
    CALL MPI_RECV(dummy, 1, MPI_INTEGER, MPI_ANY_SOURCE, tag, MPI_COMM_WORLD, &
                  status, err)
    dest = status(MPI_SOURCE)
    CALL MPI_SEND(-1, 1, MPI_INTEGER, dest, tag, MPI_COMM_WORLD, err)
  ENDDO

  CALL MPI_FINALIZE(err)
END PROGRAM Master


text:
      program DRSNML
c>> 2001-05-25 DRSCDNML Krogh Added comma to format.
c>> 1996-05-28 DRSCDNML Krogh  Changes to use M77CON
c>> 1994-07-06 DRSCDNML WV Snyder JPL set up for CHGTYP
c>> 1994-04-12 DRSCDNML WV Snyder JPL repair format to display sign
c
c     Evaluate the Cumulative Normal Distribution using SCDNML.
c
c--S replaces "?": DR?NML, DR?CDNML, ?CDNML, ?PPNML
      real             SCDNML, SPPNML
      external SCDNML, SPPNML
      real             X, C, P, MU, SIGMA
c
10    format ('        X       C=SCDNML(X)     SPPNML(C)')
20    format (1x,1p,2g14.7,2x,g14.7)
      x = -4.0e0
      mu = 0.0e0
      sigma = 1.0e0
      print 10
30    if (x .le. 4.0e0) then
         c = scdnml(x,mu,sigma)
         p = sppnml(c,mu,sigma)
         print 20, x, c, p
         x = x + 0.5e0
         go to 30
      end if
      stop
      end


text:
      INCLUDE 'VICMAIN_FOR'

      SUBROUTINE MAIN44

C-----PROGRAM LABSWTCH                  LABSWTCH (A,B) C SIZE
C-----THIS PROGRAM WILL OUTPUT THE DATA OF THE 2ND INPUT WITH
C-----THE USER LABELS OF THE 1ST INPUT, AND THE SYSTEM LABEL
C-----DETERMINED BY THE 2ND INPUT OR THE SIZE FIELD.
C-----IF THE 1ST AND 2ND INPUTS ARE NOT THE SAME SIZE, THE SIZE
C-----FIELD MUST BE USED OR THE OUTPUT MUST BE PREALLOCATED.

C   REVISION HISTORY
C      3-87   SP  ADDED XVSIGNAL CALLS FOR ERROR CHECKING (FOR FR 14986).
C      9-94   AS  (CRI) MSTP S/W CONVERSION (VICAR PORTING)
C
      CHARACTER*32 FORMAT
      INTEGER*2 BUF(16000)
      INTEGER SL,SS,NLO,NSO,NLI,NSI,NL,NS,LASTLINE,STATUS
      INTEGER INUNIT(2),OUTUNIT
C
      CALL IFMESSAGE('LABSWTCH version 05-SEP-94')
      CALL XVEACTION('SA',' ')
C--Open input files
      DO I = 1,2
          CALL XVUNIT(INUNIT(I),'INP',I,STATUS,' ')
           CALL XVSIGNAL(INUNIT( I), STATUS, .TRUE.)
          CALL XVOPEN(INUNIT(I),STATUS,' ')
           CALL XVSIGNAL(INUNIT( I), STATUS, .TRUE.)
      END DO
C
C--Get SIZE values if specified and determine number of
C--lines to be copied
      CALL XVSIZE(SL,SS,NLO,NSO,NLI,NSI)
C
C--Get system information from 2nd input
C--Get FORMAT, NL, and NS from internal control block
      CALL XVGET(INUNIT(2),STATUS,'FORMAT',FORMAT,'NL',NL,'NS',NS,' ')
       CALL XVSIGNAL(INUNIT( 2), STATUS, .TRUE.)
      IF (NLO .NE. NLI) NL = NLO
      IF (NSO .NE. NSI) NS = NSO
C
C--Now produce output file with system label of 2nd input and history
C--labels of the 1st input.
C
C--Open output with specified FORMAT, NL, and NS
      CALL XVUNIT(OUTUNIT,'OUT',1,STATUS,' ')
       CALL XVSIGNAL(OUTUNIT, STATUS, .TRUE.)
      CALL XVOPEN(OUTUNIT,STATUS,'OP','WRITE','U_NL',NL,'U_NS',NS,
     +            'U_FORMAT',FORMAT,'O_FORMAT',FORMAT,' ')
       CALL XVSIGNAL(OUTUNIT, STATUS, .TRUE.)
C
C--Copy dataset
      LASTLINE = SL + NL - 1
      DO LINE = SL, LASTLINE
          CALL XVREAD(INUNIT(2),BUF,STATUS,'LINE',LINE,'SAMP',SS,
     +                'NSAMPS',NS,' ')
            CALL XVSIGNAL(INUNIT( 2), STATUS, .TRUE.)
          CALL XVWRIT(OUTUNIT,BUF,STATUS,'LINE',LINE,'SAMP',SS,
     +                'NSAMPS',NS,' ')
            CALL XVSIGNAL(OUTUNIT, STATUS, .TRUE.)
      END DO
C
C--Close files and return
      CALL XVCLOSE(INUNIT(1),STATUS,' ')
       CALL XVSIGNAL(INUNIT(1), STATUS, .TRUE.)
      CALL XVCLOSE(INUNIT(2),STATUS,' ')
       CALL XVSIGNAL(INUNIT(2), STATUS, .TRUE.)
      CALL XVCLOSE(OUTUNIT,STATUS,' ')
       CALL XVSIGNAL(OUTUNIT, STATUS, .TRUE.)
      RETURN
      END

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
      PROGRAM xrc
C     driver for routine rc
      INTEGER i,nval
      REAL rc,value,x,y
      CHARACTER text*31
      open(7,file='FNCVAL.DAT',status='OLD')
10    read(7,'(a)') text
      if (text.ne.'Elliptic Integral Degenerate RC') goto 10
      read(7,*) nval
      write(*,*) text
      write(*,'(1x,t7,a1,t15,a1,t28,a6,t47,a7)')
     *     'X','Y','Actual','RC(X,Y)'
      do 11 i=1,nval
        read(7,*) x,y,value
        write(*,'(2f8.2,2e20.6)') x,y,value,rc(x,y)
11    continue
      close(7)
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xrc.for


! -*- f90 -*-
!****************************************
! Routines for updating lam

SUBROUTINE TrustRegion(n,m, fvec, fjac, dtd, delta, lam)
  !! Calls dgqt supplied by minpack to calculate the step and Lagrange multiplier
  IMPLICIT NONE
  INTEGER n, m, i, itmax, info
  REAL (KIND=8) fvec(m), fjac(m,n), dtd(n,n), delta, lam, v(n)
  REAL (KIND=8) z(n), wa1(n), wa2(n) !! Work arrays for dgqt
  REAL (KIND=8) rtol, atol, f
  REAL (KIND=8) jtilde(m,n), gradCtilde(n), g(n,n)

  !! Parameters for dgqt
  rtol = 1.0d-03
  atol = 1.0d-03
  itmax = 10

  DO i = 1,n
     jtilde(:,i) = fjac(:,i)/SQRT(dtd(i,i)) !! This assumes that dtd is diagonal...
  END DO
  gradCtilde = MATMUL(fvec, jtilde)
  g = MATMUL(TRANSPOSE(jtilde), jtilde)
  CALL dgqt(n, g, n, gradCtilde, delta, rtol, atol, itmax, lam, f, v, info, i, z, wa1, wa2)
  !! Transform v back to non-dtd units
!!$  DO i = 1,n
!!$     v(i) = v(i)/SQRT(dtd(i,i))
!!$  END DO
  RETURN
END SUBROUTINE TrustRegion
!! traditional update methods

SUBROUTINE Updatelam_factor(lam, accepted, factoraccept, factorreject)
  !! Update lam based on accepted/rejected step
  IMPLICIT NONE
  INTEGER accepted
  REAL (KIND=8) lam, factoraccept, factorreject

  IF(accepted.GE.0) THEN
     lam = lam / factoraccept
  ELSE
     lam = lam * factorreject
  ENDIF
END SUBROUTINE Updatelam_factor

SUBROUTINE Updatelam_nelson(lam, accepted, factoraccept, factorreject, rho)
  !! Update method due to Nelson [ref]
  IMPLICIT NONE
  INTEGER accepted, i
  DOUBLE PRECISION lam, factoraccept, factorreject, nu, rho
  IF(accepted.GE.0) THEN
     lam = lam * MAX( 1.0d+0/factoraccept, 1.0d+0 - (2.0d+0*(rho - 0.5d+0))**3 )
  ELSE
     nu = factorreject
     DO i = 1,-1*accepted !! double nu for each rejection
        nu = nu*2.0d+0
     END DO
     lam = lam * nu
  ENDIF
END SUBROUTINE Updatelam_nelson

SUBROUTINE Updatelam_Umrigar(m,n,lam, accepted, v, vold, fvec, fjac, dtd, a_param,C,Cnew)
  !! Method due to Umrigar and Nightingale [unpublished]
  IMPLICIT NONE
  INTEGER n,m,accepted, info
  DOUBLE PRECISION lam,v(n),vold(n),fvec(m), fjac(m,n), g(n,n), dtd(n,n), C, Cnew
  DOUBLE PRECISION lamold, d1, d2, grad(n), factor, a_param
  DOUBLE PRECISION amemory, cos_on

  amemory = EXP(-1.0d+0/5.0d+0)
  cos_on = DOT_PRODUCT(v,MATMUL(dtd, vold))
  cos_on = cos_on/SQRT( DOT_PRODUCT(v, MATMUL(dtd,v))*DOT_PRODUCT(vold, MATMUL(dtd,vold)))
  IF( accepted.GE.0) THEN

     IF( Cnew .LE. C ) THEN
        IF( cos_on .GT. 0) THEN
           a_param = amemory*a_param + 1.0 - amemory
        ELSE
           a_param =  amemory*a_param + 0.5*(1.0 - amemory)
        END IF
     ELSE
        a_param =  amemory*a_param + 0.5*(1.0 - amemory)
     END IF

     factor = MIN( 100.0d+0, MAX(1.1d+0, 1.0d+0/(2.2D-16 + 1.0d+0-ABS(2.0d+0*a_param - 1.0d+0))**2))
     IF (Cnew .LE. C .AND. cos_on .GE. 0) THEN
        lam = lam/factor
     ELSEIF(Cnew .GT. C) THEN
        lam = lam*SQRT(factor)
     END IF

  ELSE
     a_param =  amemory*a_param
     factor = MIN( 100.0d+0, MAX(1.1d+0, 1.0d+0/(2.2D-16 + 1.0d+0-ABS(2.0d+0*a_param - 1.0d+0))**2))
     lamold = lam
     IF( cos_on .GT. 0 ) THEN
        lam = lam * SQRT(factor)
     ELSE
        lam = lam * factor
     END IF
     ! Check for a 10% change in drift
     ! Umrigar and Nightingal suggest a check that the the proposed change in lam actually produces a meaningful change in the step.
     ! But this code produces strange results in a few cases.  -MKT
!!$     IF(accepted .EQ. -1) THEN
!!$        d1 = SQRT(DOT_PRODUCT(v,v))
!!$        g = MATMUL(TRANSPOSE(fjac),fjac) + lam*dtd
!!$        CALL DPOTRF('U', n, g, n, info)
!!$        grad = MATMUL(fvec, fjac)
!!$        CALL DPOTRS('U', n, 1, g, n, grad, n, info)
!!$        d2 = SQRT(DOT_PRODUCT(grad, grad) )
!!$        IF( 10.0d+0*ABS(d2-d1) .LT. d2 ) lam = lam - 0.1*d2*(lamold - lam)/(d1 - d2)
!!$     END IF
  ENDIF
END SUBROUTINE Updatelam_Umrigar

!! Trust region update methods

SUBROUTINE Updatedelta_factor(delta, accepted, factoraccept, factorreject)
  !! Update lam based on accepted/rejected step
  IMPLICIT NONE
  INTEGER accepted
  REAL (KIND=8) delta, factoraccept, factorreject

  IF(accepted.GE.0) THEN
     delta = delta * factoraccept
  ELSE
     delta = delta / factorreject
  ENDIF
END SUBROUTINE Updatedelta_factor

SUBROUTINE Updatedelta_more(delta, lam, n, x, dtd, rho, C, Cnew, dirder, actred, av, avmax)
  IMPLICIT NONE
  INTEGER n
  DOUBLE PRECISION delta, lam, x(n), dtd(n,n), rho, C, Cnew, dirder, actred, av, avmax
  DOUBLE PRECISION pnorm, temp
  pnorm = SQRT(DOT_PRODUCT(x,MATMUL(dtd, x)))
  IF (rho .GT. 0.25d+0) THEN
     IF (lam .GT. 0.0d+0 .AND. rho .LT. 0.75d+0) THEN
        temp = 1.0d+0
     ELSE
        temp = 2.0d+0*pnorm/delta
     END IF
  ELSE
     IF ( actred .GE. 0.0d+0) THEN
        temp = 0.5d+0
     ELSE
        temp = 0.5d+0*dirder/(dirder + 0.5d+0*actred)
     END IF
     IF ( 0.01*Cnew .GE. C .OR. temp .LT. 0.1d+0) temp = 0.1d+0
  END IF
  !! We need to make sure that if acceleration is too big, we decrease teh step size
  IF (av .GT. avmax) THEN
     temp = MIN(temp,MAX(avmax/av,0.1d+0))
  END IF

  delta = temp*MIN(delta,10.0d+0*pnorm)
  lam = lam/temp
END SUBROUTINE Updatedelta_more


text:
c
c
c     =====================================================
      subroutine rpt2(ixy,maxm,meqn,mwaves,mbc,mx,
     &                  ql,qr,aux1,aux2,aux3,
     &                  imp,asdq,bmasdq,bpasdq)
c     =====================================================
      implicit double precision (a-h,o-z)
c
c     # Riemann solver in the transverse direction for 2D Burgers' equation
c     #  u_t + cos(theta)*(0.5*u^2)_x + sin(theta)*(0.5*u^2)_y = 0
c
c     # Split asdq into eigenvectors of Roe matrix B.
c     # For the scalar equation, this simply amounts to computing the
c     # transverse wave speed from the opposite Riemann problem.
c
      dimension    ql(1-mbc:maxm+mbc, meqn)
      dimension    qr(1-mbc:maxm+mbc, meqn)
      dimension   asdq(1-mbc:maxm+mbc, meqn)
      dimension bmasdq(1-mbc:maxm+mbc, meqn)
      dimension bpasdq(1-mbc:maxm+mbc, meqn)
      common /comrp/ theta
c
      if (ixy .eq. 1) then
          b = 0.5d0*dsin(theta)
        else
          b = 0.5d0*dcos(theta)
        endif
c
          do 10 i = 2-mbc, mx+mbc
             sb = b*(qr(i-1,1) + ql(i,1))
             bmasdq(i,1) = dmin1(sb, 0.d0) * asdq(i,1)
             bpasdq(i,1) = dmax1(sb, 0.d0) * asdq(i,1)
   10        continue
c
      return
      end

max_issues_repo_name: che-wenchao/D-Claw

max_forks_repo_name: che-wenchao/D-Claw


text:
! Program to test the shape intrinsic
program testbounds
   implicit none
   real, dimension(:, :), allocatable :: a
   integer, dimension(2) :: j
   integer i

   allocate (a(3:8, 6:7))

   j = shape (a);
   if (any (j .ne. (/ 6, 2 /))) call abort

   call test(a)
contains

subroutine test (a)
   real, dimension (1:, 1:) :: a

   if (any (shape (a) .ne. (/ 6, 2 /))) call abort
end subroutine
end program

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


! { dg-do run }
!
! Test the fix for PR83344.
!
! Contributed by Janne Blomqvist  <jb@gcc.gnu.org>
! and Steve Kargl  <kargl@gcc.gnu.org>
!
program foo
   implicit none
   character(len=1) a
   character(len=2) b
   character(len=3) c
   a = 'a'
   call bah(a, len (a))
   b = 'bb'
   call bah(b, len (b))
   c = 'ccc'
   call bah(c, len (c))
   contains
      subroutine bah(x, clen)
         implicit none
         integer :: clen
         character(len=*), intent(in) :: x
         associate(y => x)
            if (len(y) .ne. clen) stop 1
            if (y .ne. x) stop 2
         end associate
      end subroutine bah
end program foo


text:
module fluid_variables
    implicit none
    save

    integer,parameter :: ndfpad=5,nsdpad=3,nenpad=8,nquadpad=8
    integer :: iquad, nquad, nquad2d
    real(8) :: sq(0:nsdpad,nenpad,nquadpad)
    real(8) :: xq(nsdpad,nquadpad),wq(nquadpad)
    real(8) :: sq2d(0:nsdpad,nenpad,nquadpad*6)
    real(8) :: xq2d(nsdpad,nquadpad*6),wq2d(nquadpad*6)
    integer :: maxconn
    real(8) :: t_start,alpha,res_g,del_g,res_l,del_l,turb_kappa
    real(8) :: ref_lgt,ref_vel,ref_den,vis_liq,den_liq
    real(8) :: liq,vmin,vmax,hmin,hmax
    integer,parameter :: maxnsurf = 21
    integer :: bc(ndfpad,maxnsurf)
    real(8) :: bv(ndfpad,maxnsurf),ic(ndfpad)
    integer :: bcd(nsdpad,maxnsurf)
    real(8) :: bvd(nsdpad,maxnsurf),icd(nsdpad)
    real(8) :: landa_over_mu
    integer :: mapping(6,8,8)
    real(8) :: interface(1:nsdpad),gravity(1:nsdpad),delta(0:21)
    integer :: etype,inner,outer,iscaling,kinner,kouter
    logical :: hg_vol,static,taudt,stokes,steady,conserve
    integer :: restart
    logical :: twod
    logical :: calcforce
    integer :: nn,ne,nq,nen,ndf,nsd,nrng,neface,nnface
    integer :: iit,nit,idisk
    integer,parameter :: tri=1,qud=2,tet=3,hex=4,tris=5,quds=6,tets=7,hexs=8
    integer :: udf,vdf,wdf,pdf
    real(8),parameter :: epsilon = 1.0e-12
    integer :: ne_inflow, edge_inflow ! Neumann "pressure" boundary condition
    real(8) pin                       ! inflow pressure
    integer ptotflag ! flag for use total pressure BC or not 1---> yes, 0 ---> no
    integer outedge
    integer nbcTimeVar,ndfTimeVar     ! #BC and #dof of time-varying BC

end module fluid_variables

max_stars_repo_name: biofluids/IFEM-archive

max_issues_repo_name: biofluids/IFEM-archive

max_forks_repo_name: biofluids/IFEM-archive


text:
      PROGRAM xdecchk
C     driver for routine decchk
      INTEGER j,k,l,n,nbad,ntot
      LOGICAL decchk,iok,jok
      CHARACTER lin*128,ch*1,chh*1,tf*1
C     test all jump transpositions of the form 86jlk41
      ntot=0
      nbad=0
      do 13 j=48,57
        do 12 k=48,57
          do 11 l=48,57
            if (j.ne.k) then
              ntot=ntot+1
              lin(1:7)='86'//char(j)//char(l)//char(k)//'41'
              iok=decchk(lin(1:7),7,ch)
              iok=decchk(lin(1:7)//ch,8,chh)
              lin(1:7)='86'//char(k)//char(l)//char(j)//'41'
              jok=decchk(lin(1:7)//ch,8,chh)
              if ((.not.iok).or.jok) then
                nbad=nbad+1
              endif
            endif
11        continue
12      continue
13    continue
      write(*,'(1x,a,t30,i3)') 'Total tries:',ntot
      write(*,'(1x,a,t30,i3)') 'Bad tries:',nbad
      write(*,'(1x,a,t30,f3.2)') 'Fraction good:',
     *  float(ntot-nbad)/float(ntot)
C     construct check digits for some user-supplied strings
1     write(*,*) 'enter string terminated by x:'
      read(*,'(a20)',END=999) lin
      do 14 j=1,128
        if (lin(j:j).eq.'x') goto 2
14    continue
2     n=j-1
      if (n.eq.0) goto 999
      iok=decchk(lin(1:n),n,ch)
      jok=decchk(lin(1:n)//ch,n+1,chh)
      if (jok) then
        tf='T'
      else
        tf='F'
      endif
      write(*,*) lin(1:n)//ch,' checks as ',tf
      goto 1
999   write(*,*) 'NORMAL COMPLETION'
      STOP
      END

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xdecchk.for


text:
! { dg-options "-O2 -fgraphite-identity -fno-loop-block -fno-loop-interchange -fno-loop-strip-mine" }

MODULE beta_gamma_psi
  INTEGER, PARAMETER :: dp=KIND(0.0D0)
CONTAINS
  FUNCTION basym (a, b, lambda, eps) RESULT(fn_val)
    REAL(dp) :: a0(21), b0(21), bsum, c(21), d(21), dsum, &
         j0, j1, r, r0, r1, s, sum, t, t0, t1, &
         u, w, w0, z, z0, z2, zn, znm1
    DO n = 2, num, 2
       DO i = n, np1
          b0(1) = r*a0(1)
          DO m = 2, i
             bsum = 0.0e0_dp
             mm1 = m - 1
             DO j = 1, mm1
                mmj = m - j
                bsum = bsum + (j*r - mmj)*a0(j)*b0(mmj)
             END DO
             b0(m) = r*a0(m) + bsum/m
          END DO
          c(i) = b0(i)/(i + 1.0e0_dp)
          d(i) = -(dsum + c(i))
       END DO
       t0 = d(n)*w*j0
       sum = sum + (t0 + t1)
    END DO
    fn_val = e0*t*u*sum
  END FUNCTION basym
END MODULE beta_gamma_psi
! { dg-final { cleanup-modules "beta_gamma_psi" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
Web hosting is a service that allows a website to be accessible via the World Wide Web. Several different such businesses exist in Davis offer web hosting services. Most of these places keep their servers elsewhere, usually in Sacramento (Herakles data center being one example) or in the Bay Area. A few people have servers at Omsoft, who have a datacenter on 5th St.

##
##
##    STOP AND READ THIS:
##
##   This is a wiki for and about Davis California.  If you are not based in or physically connected to Davis,
##   then this is not an appropriate place to list your web hosting company.  It will be removed, as you are
##   spamming a community resource.
##
##   Vandalizing public works makes your company look like a jackass, and the discussion about how youre
##   not a local company and were trying to abuse local public resources will stick in peoples minds far more
##   than anything else.
##
##   Dont be a jackass.  Dont spam.
##
##

  Access Computer Solutions  Computer repair, networking, cloudbased services for business and consumer clients. Wed design and web hosting. Onsite or pickup and drop off services available.
  Angry Sam Productions  Offers web hosting services. Headquarters in Davis. Datacenters located in Seattle, Washington DC, Dallas, and Houston.
  Baytel  Manages websites for customers on their own server in Davis and at other hosts.
  Campus Data Center  This is where most ucdavis.edu web pages are served out of. Used by university departments and such.
  http://cloudbeast.net CloudBeast  Shared Webhosting for small businesses, nonprofits and individualsm, offers 1click install for various  web scripts.
  Davis Community Network  Offers webhosting for individuals and nonprofits through Omsoft.
  Omsoft  Offers webhosting and colocation on 5th st.
(http://www.cedant.com Cedant Web Hostings tech support office lives in Davis.)


text:
The City of Davis gives out http://publicworks.cityofdavis.org/generalnotices/19thannualenvironmentalrecognitionaward Environmental Recognition Awards. There are three categories for this award: Individual, Business and NonProfit Organizations NonProfit Organization. There is a form on the website to submit a nomination. More details can be found at the website.

 Past Winners
 Individuals
  1995  Robin Kulakow
  1996  Frank Maurer
  1997  Robert L. Bugg
  1998  Steve Chainey & Bob Cordrey
  1999  James Zanetto
  2000  John Kemper
  2001  Joe Krovoza & Oam Nieberg
  2002  Max Cadji
  2003  Ron & Petra Unger
  2004  Dorothy Peterson
  2005  Richard Marovich
  2006  Dave Feliz
  2007  Charles Rominger
  2008  Emily Griswold

 Businesses
  1995  Bridge Builders Group, Inc.
  1996  Davis Energy Group
  1997  Davis Food Coop
  1998  Tandem Properties, Inc.
  1999  Calgene LLC
  2001  Davis Food Coop
  2003  Screaming Squeegee Screaming Squeegee Screen Printing & Embroidery
  2004  Sunmart, Inc.
  2005  Harrington Place
  2006  Island InkJet Island Ink Jet

  2008  MAK Design+Build, Inc.

 NonProfit Organizations
  1996  Riparian Improvement Organization
  1997  Putah Creek Putah Creek Council
  1998  Yolo Land Trust
  1999  Duck Days Steering Committee
  2000  Environmental Policy and Planning Commission ASUCD Environmental Policy & Planning Commission, and Friends of the Davis Public Library
  2001  UCD  R4 Recycling R4 Recycling Program
  2002  Yolo Basin Foundation
  2003  ASUCD (for Unitrans bus system)
  2004  Association for Efficient Environmental Energy Systems
  2005  Yolo Transportation Management Association
  2006  TREE Davis
  2007  The Unitarian Universalist Church of Davis
  2008  Davis AYSO


!*==dlasq2.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated with the qd Array Z to high relative accuracy. Used by sbdsqr and sstegr.
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download DLASQ2 + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE DLASQ2( N, Z, INFO )
!
!       .. Scalar Arguments ..
!       INTEGER            INFO, N
!       ..
!       .. Array Arguments ..
!       DOUBLE PRECISION   Z( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> DLASQ2 computes all the eigenvalues of the symmetric positive
!> definite tridiagonal matrix associated with the qd array Z to high
!> relative accuracy are computed to high relative accuracy, in the
!> absence of denormalization, underflow and overflow.
!>
!> To see the relation of Z to the tridiagonal matrix, let L be a
!> unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
!> let U be an upper bidiagonal matrix with 1's above and diagonal
!> Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
!> symmetric tridiagonal to which it is similar.
!>
!> Note : DLASQ2 defines a logical variable, IEEE, which is true
!> on machines which follow ieee-754 floating-point standard in their
!> handling of infinities and NaNs, and false otherwise. This variable
!> is passed to DLASQ3.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>        The number of rows and columns in the matrix. N >= 0.
!> \endverbatim
!>
!> \param[in,out] Z
!> \verbatim
!>          Z is DOUBLE PRECISION array, dimension ( 4*N )
!>        On entry Z holds the qd array. On exit, entries 1 to N hold
!>        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
!>        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
!>        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
!>        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
!>        shifts that failed.
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>        = 0: successful exit
!>        < 0: if the i-th argument is a scalar and had an illegal
!>             value, then INFO = -i, if the i-th argument is an
!>             array and the j-entry had an illegal value, then
!>             INFO = -(i*100+j)
!>        > 0: the algorithm failed
!>              = 1, a split was marked by a positive value in E
!>              = 2, current block of Z not diagonalized after 100*N
!>                   iterations (in inner while loop).  On exit Z holds
!>                   a qd array with the same eigenvalues as the given Z.
!>              = 3, termination criterion of outer while loop not met
!>                   (program created more than N unreduced blocks)
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup auxOTHERcomputational
!
!> \par Further Details:
!  =====================
!>
!> \verbatim
!>
!>  Local Variables: I0:N0 defines a current unreduced segment of Z.
!>  The shifts are accumulated in SIGMA. Iteration count is in ITER.
!>  Ping-pong is controlled by PP (alternates between 0 and 1).
!> \endverbatim
!>
!  =====================================================================
      SUBROUTINE DLASQ2(N,Z,Info)
      IMPLICIT NONE
!*--DLASQ2116
!
!  -- LAPACK computational routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      INTEGER Info , N
!     ..
!     .. Array Arguments ..
      DOUBLE PRECISION Z(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      DOUBLE PRECISION CBIAS
      PARAMETER (CBIAS=1.50D0)
      DOUBLE PRECISION ZERO , HALF , ONE , TWO , FOUR , HUNDRD
      PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0,FOUR=4.0D0,  &
     &           HUNDRD=100.0D0)
!     ..
!     .. Local Scalars ..
      LOGICAL ieee
      INTEGER i0 , i1 , i4 , iinfo , ipn4 , iter , iwhila , iwhilb , k ,&
     &        kmin , n0 , n1 , nbig , ndiv , nfail , pp , splt , ttype
      DOUBLE PRECISION d , dee , deemin , desig , dmin , dmin1 , dmin2 ,&
     &                 dn , dn1 , dn2 , e , emax , emin , eps , g ,     &
     &                 oldemn , qmax , qmin , s , safmin , sigma , t ,  &
     &                 tau , temp , tol , tol2 , trace , zmax , tempe , &
     &                 tempq
!     ..
!     .. External Subroutines ..
      EXTERNAL DLASQ3 , DLASRT , XERBLA
!     ..
!     .. External Functions ..
      INTEGER ILAENV
      DOUBLE PRECISION DLAMCH
      EXTERNAL DLAMCH , ILAENV
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS , DBLE , MAX , MIN , SQRT
!     ..
!     .. Executable Statements ..
!
!     Test the input arguments.
!     (in case DLASQ2 is not called by DLASQ1)
!
      Info = 0
      eps = DLAMCH('Precision')
      safmin = DLAMCH('Safe minimum')
      tol = eps*HUNDRD
      tol2 = tol**2
!
      IF ( N<0 ) THEN
         Info = -1
         CALL XERBLA('DLASQ2',1)
         RETURN
      ELSEIF ( N==0 ) THEN
         RETURN
      ELSEIF ( N==1 ) THEN
!
!        1-by-1 case.
!
         IF ( Z(1)<ZERO ) THEN
            Info = -201
            CALL XERBLA('DLASQ2',2)
         ENDIF
         RETURN
      ELSEIF ( N==2 ) THEN
!
!        2-by-2 case.
!
         IF ( Z(1)<ZERO ) THEN
            Info = -201
            CALL XERBLA('DLASQ2',2)
            RETURN
         ELSEIF ( Z(2)<ZERO ) THEN
            Info = -202
            CALL XERBLA('DLASQ2',2)
            RETURN
         ELSEIF ( Z(3)<ZERO ) THEN
            Info = -203
            CALL XERBLA('DLASQ2',2)
            RETURN
         ELSEIF ( Z(3)>Z(1) ) THEN
            d = Z(3)
            Z(3) = Z(1)
            Z(1) = d
         ENDIF
         Z(5) = Z(1) + Z(2) + Z(3)
         IF ( Z(2)>Z(3)*tol2 ) THEN
            t = HALF*((Z(1)-Z(3))+Z(2))
            s = Z(3)*(Z(2)/t)
            IF ( s<=t ) THEN
               s = Z(3)*(Z(2)/(t*(ONE+SQRT(ONE+s/t))))
            ELSE
               s = Z(3)*(Z(2)/(t+SQRT(t)*SQRT(t+s)))
            ENDIF
            t = Z(1) + (s+Z(2))
            Z(3) = Z(3)*(Z(1)/t)
            Z(1) = t
         ENDIF
         Z(2) = Z(3)
         Z(6) = Z(2) + Z(1)
         RETURN
      ENDIF
!
!     Check for negative data and compute sums of q's and e's.
!
      Z(2*N) = ZERO
      emin = Z(2)
      qmax = ZERO
      zmax = ZERO
      d = ZERO
      e = ZERO
!
      DO k = 1 , 2*(N-1) , 2
         IF ( Z(k)<ZERO ) THEN
            Info = -(200+k)
            CALL XERBLA('DLASQ2',2)
            RETURN
         ELSEIF ( Z(k+1)<ZERO ) THEN
            Info = -(200+k+1)
            CALL XERBLA('DLASQ2',2)
            RETURN
         ENDIF
         d = d + Z(k)
         e = e + Z(k+1)
         qmax = MAX(qmax,Z(k))
         emin = MIN(emin,Z(k+1))
         zmax = MAX(qmax,zmax,Z(k+1))
      ENDDO
      IF ( Z(2*N-1)<ZERO ) THEN
         Info = -(200+2*N-1)
         CALL XERBLA('DLASQ2',2)
         RETURN
      ENDIF
      d = d + Z(2*N-1)
      qmax = MAX(qmax,Z(2*N-1))
      zmax = MAX(qmax,zmax)
!
!     Check for diagonality.
!
      IF ( e==ZERO ) THEN
         DO k = 2 , N
            Z(k) = Z(2*k-1)
         ENDDO
         CALL DLASRT('D',N,Z,iinfo)
         Z(2*N-1) = d
         RETURN
      ENDIF
!
      trace = d + e
!
!     Check for zero data.
!
      IF ( trace==ZERO ) THEN
         Z(2*N-1) = ZERO
         RETURN
      ENDIF
!
!     Check whether the machine is IEEE conformable.
!
      ieee = (ILAENV(10,'DLASQ2','N',1,2,3,4)==1)
!
!     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
!
      DO k = 2*N , 2 , -2
         Z(2*k) = ZERO
         Z(2*k-1) = Z(k)
         Z(2*k-2) = ZERO
         Z(2*k-3) = Z(k-1)
      ENDDO
!
      i0 = 1
      n0 = N
!
!     Reverse the qd-array, if warranted.
!
      IF ( CBIAS*Z(4*i0-3)<Z(4*n0-3) ) THEN
         ipn4 = 4*(i0+n0)
         DO i4 = 4*i0 , 2*(i0+n0-1) , 4
            temp = Z(i4-3)
            Z(i4-3) = Z(ipn4-i4-3)
            Z(ipn4-i4-3) = temp
            temp = Z(i4-1)
            Z(i4-1) = Z(ipn4-i4-5)
            Z(ipn4-i4-5) = temp
         ENDDO
      ENDIF
!
!     Initial split checking via dqd and Li's test.
!
      pp = 0
!
      DO k = 1 , 2
!
         d = Z(4*n0+pp-3)
         DO i4 = 4*(n0-1) + pp , 4*i0 + pp , -4
            IF ( Z(i4-1)<=tol2*d ) THEN
               Z(i4-1) = -ZERO
               d = Z(i4-3)
            ELSE
               d = Z(i4-3)*(d/(d+Z(i4-1)))
            ENDIF
         ENDDO
!
!        dqd maps Z to ZZ plus Li's test.
!
         emin = Z(4*i0+pp+1)
         d = Z(4*i0+pp-3)
         DO i4 = 4*i0 + pp , 4*(n0-1) + pp , 4
            Z(i4-2*pp-2) = d + Z(i4-1)
            IF ( Z(i4-1)<=tol2*d ) THEN
               Z(i4-1) = -ZERO
               Z(i4-2*pp-2) = d
               Z(i4-2*pp) = ZERO
               d = Z(i4+1)
            ELSEIF ( safmin*Z(i4+1)<Z(i4-2*pp-2) .AND.                  &
     &               safmin*Z(i4-2*pp-2)<Z(i4+1) ) THEN
               temp = Z(i4+1)/Z(i4-2*pp-2)
               Z(i4-2*pp) = Z(i4-1)*temp
               d = d*temp
            ELSE
               Z(i4-2*pp) = Z(i4+1)*(Z(i4-1)/Z(i4-2*pp-2))
               d = Z(i4+1)*(d/Z(i4-2*pp-2))
            ENDIF
            emin = MIN(emin,Z(i4-2*pp))
         ENDDO
         Z(4*n0-pp-2) = d
!
!        Now find qmax.
!
         qmax = Z(4*i0-pp-2)
         DO i4 = 4*i0 - pp + 2 , 4*n0 - pp - 2 , 4
            qmax = MAX(qmax,Z(i4))
         ENDDO
!
!        Prepare for the next iteration on K.
!
         pp = 1 - pp
      ENDDO
!
!     Initialise variables to pass to DLASQ3.
!
      ttype = 0
      dmin1 = ZERO
      dmin2 = ZERO
      dn = ZERO
      dn1 = ZERO
      dn2 = ZERO
      g = ZERO
      tau = ZERO
!
      iter = 2
      nfail = 0
      ndiv = 2*(n0-i0)
!
      DO iwhila = 1 , N + 1
         IF ( n0<1 ) GOTO 200
!
!        While array unfinished do
!
!        E(N0) holds the value of SIGMA when submatrix in I0:N0
!        splits from the rest of the array, but is negated.
!
         desig = ZERO
         IF ( n0==N ) THEN
            sigma = ZERO
         ELSE
            sigma = -Z(4*n0-1)
         ENDIF
         IF ( sigma<ZERO ) THEN
            Info = 1
            RETURN
         ENDIF
!
!        Find last unreduced submatrix's top index I0, find QMAX and
!        EMIN. Find Gershgorin-type bound if Q's much greater than E's.
!
         emax = ZERO
         IF ( n0>i0 ) THEN
            emin = ABS(Z(4*n0-5))
         ELSE
            emin = ZERO
         ENDIF
         qmin = Z(4*n0-3)
         qmax = qmin
         DO i4 = 4*n0 , 8 , -4
            IF ( Z(i4-5)<=ZERO ) GOTO 50
            IF ( qmin>=FOUR*emax ) THEN
               qmin = MIN(qmin,Z(i4-3))
               emax = MAX(emax,Z(i4-5))
            ENDIF
            qmax = MAX(qmax,Z(i4-7)+Z(i4-5))
            emin = MIN(emin,Z(i4-5))
         ENDDO
         i4 = 4
!
 50      i0 = i4/4
         pp = 0
!
         IF ( n0-i0>1 ) THEN
            dee = Z(4*i0-3)
            deemin = dee
            kmin = i0
            DO i4 = 4*i0 + 1 , 4*n0 - 3 , 4
               dee = Z(i4)*(dee/(dee+Z(i4-2)))
               IF ( dee<=deemin ) THEN
                  deemin = dee
                  kmin = (i4+3)/4
               ENDIF
            ENDDO
            IF ( (kmin-i0)*2<n0-kmin .AND. deemin<=HALF*Z(4*n0-3) ) THEN
               ipn4 = 4*(i0+n0)
               pp = 2
               DO i4 = 4*i0 , 2*(i0+n0-1) , 4
                  temp = Z(i4-3)
                  Z(i4-3) = Z(ipn4-i4-3)
                  Z(ipn4-i4-3) = temp
                  temp = Z(i4-2)
                  Z(i4-2) = Z(ipn4-i4-2)
                  Z(ipn4-i4-2) = temp
                  temp = Z(i4-1)
                  Z(i4-1) = Z(ipn4-i4-5)
                  Z(ipn4-i4-5) = temp
                  temp = Z(i4)
                  Z(i4) = Z(ipn4-i4-4)
                  Z(ipn4-i4-4) = temp
               ENDDO
            ENDIF
         ENDIF
!
!        Put -(initial shift) into DMIN.
!
         dmin = -MAX(ZERO,qmin-TWO*SQRT(qmin)*SQRT(emax))
!
!        Now I0:N0 is unreduced.
!        PP = 0 for ping, PP = 1 for pong.
!        PP = 2 indicates that flipping was applied to the Z array and
!               and that the tests for deflation upon entry in DLASQ3
!               should not be performed.
!
         nbig = 100*(n0-i0+1)
         DO iwhilb = 1 , nbig
            IF ( i0>n0 ) GOTO 100
!
!           While submatrix unfinished take a good dqds step.
!
            CALL DLASQ3(i0,n0,Z,pp,dmin,sigma,desig,qmax,nfail,iter,    &
     &                  ndiv,ieee,ttype,dmin1,dmin2,dn,dn1,dn2,g,tau)
!
            pp = 1 - pp
!
!           When EMIN is very small check for splits.
!
            IF ( pp==0 .AND. n0-i0>=3 ) THEN
               IF ( Z(4*n0)<=tol2*qmax .OR. Z(4*n0-1)<=tol2*sigma ) THEN
                  splt = i0 - 1
                  qmax = Z(4*i0-3)
                  emin = Z(4*i0-1)
                  oldemn = Z(4*i0)
                  DO i4 = 4*i0 , 4*(n0-3) , 4
                     IF ( Z(i4)<=tol2*Z(i4-3) .OR. Z(i4-1)<=tol2*sigma )&
     &                    THEN
                        Z(i4-1) = -sigma
                        splt = i4/4
                        qmax = ZERO
                        emin = Z(i4+3)
                        oldemn = Z(i4+4)
                     ELSE
                        qmax = MAX(qmax,Z(i4+1))
                        emin = MIN(emin,Z(i4-1))
                        oldemn = MIN(oldemn,Z(i4))
                     ENDIF
                  ENDDO
                  Z(4*n0-1) = emin
                  Z(4*n0) = oldemn
                  i0 = splt + 1
               ENDIF
            ENDIF
!
         ENDDO
!
         Info = 2
!
!        Maximum number of iterations exceeded, restore the shift
!        SIGMA and place the new d's and e's in a qd array.
!        This might need to be done for several blocks
!
         i1 = i0
         n1 = n0
         DO
            tempq = Z(4*i0-3)
            Z(4*i0-3) = Z(4*i0-3) + sigma
            DO k = i0 + 1 , n0
               tempe = Z(4*k-5)
               Z(4*k-5) = Z(4*k-5)*(tempq/Z(4*k-7))
               tempq = Z(4*k-3)
               Z(4*k-3) = Z(4*k-3) + sigma + tempe - Z(4*k-5)
            ENDDO
!
!        Prepare to do this on the previous block if there is one
!
            IF ( i1>1 ) THEN
               n1 = i1 - 1
               DO WHILE ( (i1>=2) .AND. (Z(4*i1-5)>=ZERO) )
                  i1 = i1 - 1
               ENDDO
               sigma = -Z(4*n1-1)
               CYCLE
            ENDIF

            DO k = 1 , N
               Z(2*k-1) = Z(4*k-3)
!
!        Only the block 1..N0 is unfinished.  The rest of the e's
!        must be essentially zero, although sometimes other data
!        has been stored in them.
!
               IF ( k<n0 ) THEN
                  Z(2*k) = Z(4*k-1)
               ELSE
                  Z(2*k) = 0
               ENDIF
            ENDDO
            RETURN
         ENDDO
!
!        end IWHILB
!
!
 100  ENDDO
!
      Info = 3
      RETURN
!
!     end IWHILA
!
!
!     Move q's to the front.
!
 200  DO k = 2 , N
         Z(k) = Z(4*k-3)
      ENDDO
!
!     Sort and compute sum of eigenvalues.
!
      CALL DLASRT('D',N,Z,iinfo)
!
      e = ZERO
      DO k = N , 1 , -1
         e = e + Z(k)
      ENDDO
!
!     Store trace, sum(eigenvalues) and information on performance.
!
      Z(2*N+1) = trace
      Z(2*N+2) = e
      Z(2*N+3) = DBLE(iter)
      Z(2*N+4) = DBLE(ndiv)/DBLE(N**2)
      Z(2*N+5) = HUNDRD*nfail/DBLE(iter)
!
!     End of DLASQ2
!
      END SUBROUTINE DLASQ2


! PR fortran/46874
! { dg-do run }

  interface
    subroutine sub (a, b, c, d, n)
      integer :: n
      integer, allocatable :: a(:), b(:), c(:), d(:)
    end subroutine
  end interface

  integer, allocatable :: a(:), b(:), c(:), d(:)
  integer :: i, j
  allocate (a(50), b(50), c(50), d(50))
  do i = 1, 50
    a(i) = 2 + modulo (i, 7)
    b(i) = 179 - modulo (i, 11)
  end do
  c = 0
  d = 2147483647
  call sub (a, b, c, d, 50)
  do i = 1, 50
    j = 0
    if (i .eq. 3) then
      j = 8
    else if (i .gt. 1 .and. i .lt. 9) then
      j = 7
    end if
    if (c(i) .ne. j) call abort
    j = 179 - modulo (i, 11)
    if (i .gt. 1 .and. i .lt. 9) j = i
    if (d(i) .ne. j) call abort
  end do
  deallocate (a, b, c, d)
end

subroutine sub (a, b, c, d, n)
  integer :: n
  integer, allocatable :: a(:), b(:), c(:), d(:)
!$omp parallel do shared(a, b) reduction(+:c) reduction(min:d)
  do i = 1, n
    c(a(i)) = c(a(i)) + 1
    d(i) = min(d(i), b(i))
    d(a(i)) = min(d(a(i)), a(i))
  end do
end


    Program dgesvj_example

!     DGESVJ Example Program Text

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     .. Use Statements ..
      Use lapack_example_aux, Only: nagf_file_print_matrix_real_gen
      Use lapack_interfaces, Only: ddisna, dgesvj
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: eps, serrbd
      Integer :: i, ifail, info, j, lda, ldv, lwork, m, n
!     .. Local Arrays ..
      Real (Kind=dp), Allocatable :: a(:, :), rcondu(:), rcondv(:), s(:), &
        v(:, :), work(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: abs, epsilon
!     .. Executable Statements ..
      Write (nout, *) 'DGESVJ Example Program Results'
      Write (nout, *)
      Flush (nout)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) m, n
      lda = m
      ldv = n
      lwork = n + m
      Allocate (a(lda,n), rcondu(m), rcondv(m), s(n), v(ldv,n), work(lwork))

!     Read the m by n matrix A from data file

      Read (nin, *)((a(i,j),j=1,n), i=1, m)

!     Compute the singular values and left and right singular vectors
!     of A (A = U*S*V, m.ge.n)

      Call dgesvj('G', 'U', 'V', m, n, a, lda, s, 0, v, ldv, work, lwork, &
        info)

      If (info==0) Then

!       Compute the approximate error bound for the computed singular values
!       using the 2-norm, s(1) = norm(A), and machine precision, eps.
        eps = epsilon(1.0E0_dp)
        serrbd = eps*s(1)

!       Print solution
        Write (nout, *) 'Singular values'
        Write (nout, 100)(s(j), j=1, n)

        If (abs(work(1)-1.0_dp)>eps) Then
          Write (nout, 130) 'Values need scaling by factor = ', work(1)
        End If
        Flush (nout)

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call nagf_file_print_matrix_real_gen('General', ' ', m, n, a, lda, &
          'Left singular vectors', ifail)

        Write (nout, *)
        Flush (nout)

        ifail = 0
        Call nagf_file_print_matrix_real_gen('General', ' ', n, n, v, ldv, &
          'Right singular vectors', ifail)

!       Call DDISNA to estimate reciprocal condition
!       numbers for the singular vectors
        Call ddisna('Left', m, n, s, rcondu, info)
        Call ddisna('Right', m, n, s, rcondv, info)

!       Print the approximate error bounds for the singular values
!       and vectors
        Write (nout, '(/1X,A)') 'Error estimate for the singular values'
        Write (nout, 110) serrbd
        Write (nout, '(/1X,A)') 'Error estimates for left singular vectors'
        Write (nout, 110)(serrbd/rcondu(i), i=1, n)
        Write (nout, '(/1X,A)') 'Error estimates for right singular vectors'
        Write (nout, 110)(serrbd/rcondv(i), i=1, n)
      Else
        Write (nout, 120) 'Failure in DGESVJ. INFO =', info
      End If

100   Format (3X, (8F8.4))
110   Format (4X, 1P, 6E11.1)
120   Format (1X, A, I4)
130   Format (/, 1X, A, 1P, E13.5)
    End Program


text:
program permutations

  implicit none
  integer, parameter :: value_min = 1
  integer, parameter :: value_max = 3
  integer, parameter :: position_min = value_min
  integer, parameter :: position_max = value_max
  integer, dimension (position_min : position_max) :: permutation

  call generate (position_min)

contains

  recursive subroutine generate (position)

    implicit none
    integer, intent (in) :: position
    integer :: value

    if (position > position_max) then
      write (*, *) permutation
    else
      do value = value_min, value_max
        if (.not. any (permutation (: position - 1) == value)) then
          permutation (position) = value
          call generate (position + 1)
        end if
      end do
    end if

  end subroutine generate

end program permutations

max_stars_repo_name: mullikine/RosettaCodeData

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData