text:
It happens.  You got laid off.  Or fired.  Or quit.  Whatever.  Regardless, youll most likely be facing an extreme shortage, if not a complete cessation, of incoming money.  Fortunately, there are different programs in place to help you until you can find your next job.

Luckily for Davis, the towns demographics unemployment rate has been significantly lower than the state average for many years.

 Temporary Financial Assistance
 Unemployment Benefits
http://www.edd.ca.gov/Unemployment/ Unemployment benefits are available to qualifying individuals who have recently been laid off due to no fault of their own.  Benefits come in cash (or debit card) payments, and amount to somewhere around half of your regular paycheck.  Benefits also come with the stipulation that you register with the http://www.ladpss.org/dpss/calworks/default.cfm CalWorks job search and keep track of your job search efforts.  They also have a few training opportunity programs with limited enrollment.

 General Assistance
The county provides emergency assistance to people in dire need.  The program provides financial assistance, substance abuse counseling, health care benefits and job training.  Visit the http://www.yolocounty.org/Index.aspx?page557 Yolo County General Assistance webpage for more info.

 Paying for Necessities
 General
  Davis Community Meals and STEAC offers not only food assistance, but also assistance regarding housing, family issues and substance abuse.
  Living Cheaply

 Food
  http://www.dss.cahwnet.gov/foodstamps/default.htm Supplemental Nutrition Assistance Program Formerly/colloquially known as Food Stamps, this program issues Electronic Benefit Transfers (also known as EBT cards) to use to buy groceries and food items.  They work much the same way as debit cards, and can only be used on specific food and foodrelated items.  Admission to the program occurs via http://www.dss.cahwnet.gov/foodstamps/entres/forms/english/DFA285A1.PDF application submission and an interview.
  The Food Bank of Yolo County has a number of programs to help alleviate hunger in the county.  You can visit http://www.foodbankyc.org/index.html their website for more information on distribution of food.
  Food Not Bombs can help you with a hot meal once a week in Central Park.  The Episcopal Church of St. Martin hosts hot meals to the community on for lunch on Saturday, as well as for dinner on Tuesdays and Thursdays.
  Cheap Food is a pretty nifty listing of good deals in town
  Free Food is even better.
  Free Beer

 Health Care
Termination of employment often results in termination of insurance benefits, so its good to know what other options you have available.

 MediCal
http://www.dhcs.ca.gov/services/medical/Pages/default.aspx MediCal is a statesponsored public health care program available to http://www.dhcs.ca.gov/services/medical/Pages/MediCalEligibility.aspx eligible individuals, notably caretakers (whether for children or dependent adults), pregnant women, and children.  DentiCal is available for dental care services.  http://www.dhcs.ca.gov/services/medical/Pages/MediCalApplications.aspx Download Application forms here.

 Community Clinic
CommuniCare Health Center is the local community clinic that serves low income and uninsured (or underinsured) individuals.  They accept various insurance plans, but also operate on a slidingscale payment system.

 Other Programs
  WIC Program WIC provides assistance to pregnant and breastfeeding women and infants
  QuickHealth inside the West Sacramento WalMart is a http://www.walmart.com/catalog/catalog.gsp?cat554492 nonemergency clinic designed to handle common ailments (like health screenings, wart removal, yeast infection treatment, etc) as well as general physicals.  They also offer package deals for related screenings and tests.  No insurance is needed.
  Planned Parenthood offers reproductive health education, gynecological exams, birth control, and abortion services

 Child Care
  http://headstart.seta.net/programoverview.htm HeadStart offers programs for http://headstart.seta.net/ehspregnantwomensvcs.htm expectant parents, infants/toddlers (http://headstart.seta.net/earlyheadstart.htm Early Head Start), and preschoolage children
  http://www.djusd.k12.ca.us/ChildrensCenter/ The Childrens Center is a statesubsidized preschool and childcare center with lowcost and free childcare.  The center is affiliated with the local school district.  For children ages 35.
  http://www.cdicdc.org/whiteboard/index.php CDC provides childcare and daycare for children of all ages in lowincome families
  The city also provides http://cityofdavis.org/cs/childcare/subsidy.cfm childcare subsidies to qualifying families

 LowIncome Housing
The city has a lowincome housing program for families meeting the requisite income and household size.  Essentially, qualifying households may receive http://en.wikipedia.org/wiki/Section_8_(housing) Section 8 housing vouchers to offset the cost of rent.  Section 8 vouchers are accepted at http://cityofdavis.org/housing/affordable/list.cfm a number of Davis housing complexes.  Unfortunately, there tend to be long waitlists for openings at these complexes, and it is not uncommon for the wait for an available apartment to last between six months and a year.

A more expedient option for individuals who rent housing may be to relocate to less expensive areas, especially Woodland, Dixon/Solano County, or in rural Yolo county.

 Finding Work
Check out Employment, Job Applications and Job Board

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine setboostatoms
      include 'moldy.h'
c
      integer*4 i, j, ii, ihe, ijk
      real*8 xij, yij, zij, rxij, ryij, rzij, r2
      real*8 xhea, yhea, zhea, rij
      integer*4 idya1, iikey

      ihenum = 0
      do i = 1, nn1
         if(atpe(i) .gt. -3.89)then
           ihenum = ihenum + 1
           iheid(ihenum) = i
         endif
      end do
      write(*,*) 'total number of atoms with Ep>-3.89eV is ', ihenum
c determine the surrounding atoms of these accelerating atoms
      do ihe = 1,ihenum
         ijk = iheid(ihe)
         xhea = x0(ijk)
         yhea = y0(ijk)
         zhea = z0(ijk)
c idya: the number of found near-by atoms around current He position
          idya1 = idya

c loop over all the atoms to find the atoms close to the current
c He position, distance less than 3 angstrom
          do i = 1, nn1
            xij = (x0(i)-xhea)*b0(1,1)
            yij = (y0(i)-yhea)*b0(2,2)
            zij = (z0(i)-zhea)*b0(3,3)
            if(xij.ge. b0(1,1)/2.0d0) xij = xij - b0(1,1)
            if(xij.lt.-b0(1,1)/2.0d0) xij = xij + b0(1,1)
            if(yij.ge. b0(2,2)/2.0d0) yij = yij - b0(2,2)
            if(yij.lt.-b0(2,2)/2.0d0) yij = yij + b0(2,2)
            if(zij.ge. b0(3,3)/2.0d0) zij = zij - b0(3,3)
            if(zij.lt.-b0(3,3)/2.0d0) zij = zij + b0(3,3)
            rij = sqrt(yij*yij + zij*zij)
c for loops, we consider the distance around one B
           if(abs(xij).le.0.5*b0(1,1) .and. rij .lt. 3.5d0)then
              iikey = 1
c idya_n(i): store the id of atoms near by He atoms
              do j = 1, idya1
                 iii = idya_n(j)
                 if(iii .eq. i)then
                   iikey = 0
                   exit
                 endif
              end do

              if(iikey .eq. 1)then
                idya = idya + 1
                idya_n(idya) = i
                xjp(idya) = x0(i)
                yjp(idya) = y0(i)
                zjp(idya) = z0(i)
              endif
            endif
         end do
      end do

      write(*,*) 'total accelerating atoms = ', idya

      do i = 1, idya
         ii = idya_n(i)
         x00 = xjp(i)*b0(1,1)
         y00 = yjp(i)*b0(2,2)
         z00 = zjp(i)*b0(3,3)
         write(*,*) x00, y00, z00, id(ii), ii
      end do

      open(unit=21,file='mdyn_region.psi',position='append')
      write(21,'(a)') 'ITEM: TIMESTEP'
      write(21,'(i9)') nst
      write(21,'(a)') 'ITEM: NUMBER OF ATOMS'
      write(21,'(i6)') idya
      write(21,'(a)') 'ITEM: BOX BOUNDS'
      write(21,'(2(f18.9))') -b0(1,1)/2.0d0, b0(1,1)/2.0d0
      write(21,'(2(f18.9))') -b0(2,2)/2.0d0, b0(2,2)/2.0d0
      write(21,'(2(f18.9))') -b0(3,3)/2.0d0, b0(3,3)/2.0d0
      write(21,'(a)') 'ITEM: ATOMS id x y z itype pe ke pre shr'
      do 12211 j = 1, idya
         iij = idya_n(j)
         xiij = x0(iij)*b0(1,1)
         yiij = y0(iij)*b0(2,2)
         ziij = z0(iij)*b0(3,3)
         write(21,'(i6,2x,3(f12.5),2x,i2,4(f12.5))')
     &   iij,xiij,yiij,ziij,id(iij),atpe(iij),ekin(iij),atpress(iij),
     &        atshear(iij)
12211 continue
      close(21)

      return
      end subroutine


text:
*
*
      subroutine formb3(nup)
*
*      compute formation of binaries due to three-body interactions
*      ------------------------------------------------------------
*      Goodman & Hut 1994, .....
*      -------------------------
*
*     abs(ikind(im)) - formb3 is the last routine in the interaction sequence
*     (intb3b3 - intb3f - formb3) therefore same objects could take place
*     already in interactions in intb3b3 or intb3f and have negative ikind
*
*
      include 'common.h'
*
      real*8 dt,co,vol,velave,smave,den,vb,f2,wrr,ptold,ptnew,
     &       xmm,xmm1,xmm2,age1,age2,sm1s,ehpot,ehkin,cor,timea,anb,
     &       body12,tform,tb3f,tb3b3,tcoll,tim,sm2s,ekicks1,ekicks2,
     &       ekicks3,sm1so,sm2so,zkold,zknew,ehmlevo,ehmlevn,ehbin3o,
     &       ehbin3n,ekickto,ekicktn,r1,r2,xlum1,xlum2,slum
*
      real*4 ran2
*
      integer nup,l,im1,im2,im3,ibin,lmin,lmax,n,inum,inumt,idb,
     &        id1,id2,is1,is2,is3,ik1,ik2,ik3,im,k1,k11,k2,k3,iexi,
     &        iexi1,iexi2,noibin,ik0
*
      common /timset/ tform(20),tb3f(20),tb3b3(20),tcoll(20)
*
      open(44,file=trim(datadir)//'/binarym.dat',access='append')
      open(43,file=trim(datadir)//'/starm.dat',access='append')

      n = nt
      sm2s = 0.d0
      f2 = 0.9d0
      ekisks1 = 0.d0
      ekicks2 = 0.d0
      ekicks3 = 0.d0
      pctonbu = rtidkg/rbar
      iexi = 0
      noibin = 0
*
*     Stodolkiewicz's density definition
*     cos = 4/3*pi*3^4.5/(4/3*pi)^3*f2*float(n)/log(gamma*n)*(2/3)^3
*
*     Henon's density definition
*     coh = 4/3*pi*3^4.5/(2*pi)^3*float(n)/log(gamma*n) = cos
*
      if(alphal.eq.-1.0d0) then
*
*     coefficient for single-mass case
*
        co = 2.36916d0*f2*float(nt)/log(gamma*nt)
      else
*
*     coefficient for multi-mass case
*
*     co = 3^6*f2*(m1*m2)^4*m3^2.5*nt0/(2^1.5*pi^2*((m1+m2)*(m1+m2+m3))^0.5*
*          <m>^4.5*(v^2)^4.5*log(gamma*nt0))
*
        co = 26.114565d0*f2*float(nt)/log(gamma*nt)
      endif
*
      ibin = 0
*
      do 10 l = 1,nup
*
         if(l.eq.1) then
           lmin = 1
         else
           lmin  = nzst(l-1) + 1
         endif
*
         lmax = nzst(l)
*
         if(ltwo(l).ne.0) then
           dt = tau*2.0d0**(-ltwo(l))
         else
           dt = tau
         endif
*
         tform(l) = tform(l) + dt*tscale0*log(gamma*nt00)
     &              /log(gamma*ntnew)*float(ntnew)/float(nt00)
         print*,'formb3- l,nup,tform(l) = ',l,nup,tform(l)
         call flush(6)
         tim = tform(l)
*
*      acomulate probability 'pbin3' for formation of 3-body binaries
*      find stars to form a binary
*
         k1 = lmin - 1
*
 20      continue
*
         k1 = k1 + 1
c         if(k1.gt.int(0.8*n)) go to 10
         if(k1.gt.lmax) go to 10
*
         im1 = iname(k1)
         if(abs(ikind(im1)).eq.2.or.names(im1).eq.0) go to 20
         if(abs(body(im1)).lt.tolm) go to 20
         if(nkick(im1).ge.1) go to 20
         if(r(k1).gt.1.d8) go to 20
         is1 = names(im1)
         call getLum(is1,slum)
         if(isNaN(slum)) then
           print*,'k1,im1,is1,lum = ',k1,im1,is1,slum
           go to 20
         endif
*
         k2 = k1
*
 30      continue
*
         k2 = k2 + 1
c         if(k2.gt.int(0.8*n)) go to 10
         if(k2.gt.lmax) go to 10
*
         im2 = iname(k2)
         if(abs(ikind(im2)).eq.2.or.names(im2).eq.0) go to 30
         if(abs(body(im2)).lt.tolm) go to 30
         if(nkick(im2).ge.1) go to 30
         if(r(k2).gt.1.d8) go to 30
         is1 = names(im2)
         call getLum(is1,slum)
         if(isNaN(slum)) then
           print*,'k2,im2,is2,lum = ',k2,im2,is1,slum
           go to 30
         endif
*
         k3 = k2
*
 40      continue
*
         k3 = k3 + 1
c         if(k3.gt.int(0.8*n)) go to 10
         if(k3.gt.lmax) go to 10
*
         im3 = iname(k3)
         if(abs(ikind(im3)).eq.2.or.names(im3).eq.0) go to 40
         if(abs(body(im3)).lt.tolm) go to 40
         if(nkick(im3).ge.1) go to 40
         if(r(k3).gt.1.d8) go to 40
         is1 = names(im3)
         call getLum(is1,slum)
         if(isNaN(slum)) then
           print*,'k3,im3,is3,lum = ',k3,im3,is1,slum
           go to 40
         endif
*
*       take care for multi-mass systems
*
         if(alphal.eq.-1.0d0) then
           if(k1.eq.1) then
             vol = co*r(k1+2)**3
           else
             vol = co*(r(k1+2)**3 - r(k1-1)**3)
           endif
           call veloci(k1,velave,smave)
           call densi(k1,den)
*
           pbin3 = vol*smave**5*den**3*dt/velave**4.5
*
         else
           call veloci(k1,velave,smave)
           call densi(k1,den)
cAdded DCH 28/8/6 to deal with a central triple of black holes
           if(k3.ge.3.and.k3.le.5) then
             is1 = names(im1)
             call get_ss_type(is1,ik1)
             if(ik1.ge.10.and.ik1.lt.15) then
               is2 = names(im2)
               call get_ss_type(is2,ik2)
               if(ik2.ge.10.and.ik2.lt.15) then
                 is3 = names(im3)
                 call get_ss_type(is3,ik3)
                 if(ik3.ge.10.and.ik3.lt.15) then
                   velave = body(im1)*(vr(im1)**2 + vt(im1)**2) +
     &                      body(im2)*(vr(im2)**2 + vt(im2)**2) +
     &                      body(im3)*(vr(im3)**2 + vt(im3)**2)
                   smave = body(im1) + body(im2) + body(im3)
                   velave = velave/smave
                   smave = smave/3.d0
                   den = 3.d0/(r(k3)**3 - r(k1)**3)
                   print*,'Three/or more inmost stars are WD-NS-BH',
     &                       body(im2),body(im2),body(im3),den,smave,
     &                       velave,co,dt,k3
                 endif
               endif
             endif
           endif
           if(body(im1)*body(im2)*body(im3).ne.0.d0) then
             xmm = (body(im1)*body(im2))**4*body(im3)**2.5
             xmm1 = body(im1) + body(im2)
             xmm2 = xmm1*(xmm1 + body(im3))
             xmm = xmm/(sqrt(xmm2)*(smave*velave)**4.5)
             pbin3 = co*xmm*dt*den**2
           else
             pbin3 = 0.d0
           endif
*
         endif
*
*      new check for binary formations starts from k3 + 1
*
         k11 = k1
         k1 = k3
         wrr = ran2(irun)
         if(wrr.gt.pbin3) go to 20
*
*         compute total potential before binary formation to deal with
*         potential changes due to substituting two single stars by binary
*
         if(ibin.eq.0.and.noibin.eq.0) then
           call energy(2)
           ptold = pot
           zkold = zkin
           ehmlevo = ehmlev
           ehbin3o = ehbin3
           ekickto = ekickt
         endif
*
*   evolve single stars to the time tim - time of binary formation
*
         call mloss_single(k11,tim,ekicks1,iexi)
*
         call mloss_single(k2,tim,ekicks2,iexi1)
*
         call mloss_single(k3,tim,ekicks3,iexi2)
*
*        check for obliterated stars
*
c         if(r(k11).gt.1.d8.or.r(k2).gt.1.d8.or.r(k3).gt.1.d8) then
c           print*,'formb3-obl  k1,k2,k3,im1,im2,im3,sm1,sm2,sm3 =',
c     &     k1,k2,k3,im1,im2,im3,body(im1)*zmbar,body(im2)*zmbar,
c     &     body(im3)*zmbra
c           go to 20
c         endif
*
         if(iexi.eq.0.or.iexi1.eq.0.or.iexi2.eq.0) then
           noibin = noibin + 1
           print*,'formb3-no  im1,im2,im3,iexi,iexi1 iexi2 = ',
     &             im1,im2,im3,iexi,iexi1,iexi2
           go to 20
         endif
*
c         tim = tim + 2.d-14
         pbin3 = 0.0d0
         nbin3 = nbin3 + 1
         ik0 = -1
*
         print*,'k1,velave,smave,den =',k11,velave,smave,den
         call veloc(k11,velave,smave)
         print*,'k1,velave,smave,den =',k11,velave,smave,den
*
*       define parameters of binary and compute kick velocities
*       for binary and field star
*
         vb = smave*velave
*
         ibin = ibin + 1
         ikind(im1) = 2
*
         id1 = names(im1)
         id2 = names(im2)
*
*      Update names and nameb and prepare the new binary for binary
*      evolution
*      nwhich contains the name of binary in the binary list.
*      if ibx=im1-nss0 is negative then binary index can be get from the
*      nwhich(im1)
*
*     set ibstra for newly formed binary. If ibstra(im1) or/and ibstra(im2)
*     is greater than 0 then ibstra(im1) is set to ibstra(im1) or/and
*     ibstra(im2)
*
         if(ibstra(im1).ge.ibstra(im2)) then
           ibstra(im1) = ibstra(im1)
           ibstra(im2) = 0
         else
           ibstra(im1) = ibstra(im2)
           ibstra(im2) = 0
         endif
*
         nameb(im1) = id1
         nameb(im2) = 0
         nwhich(im1) = nbin3
         nbinar(nbin3) = im1
         call get_mass(id1,sm1s)
         call get_mass(id2,sm2s)
         idb = nameb(im1)
         body12 = sm1s + sm2s
         call getAge(id1,age1)
         call getAge(id2,age2)
*
         print*,'k1,k2,k3,im1,im2,im3,id1,id2,idb,nameb,names1,',
     &   'names2,sm1b,sm2b,sm1s,sm2s body12,age1,age2 =',k11,k2,k3,im1,
     &   im2,im3,id1,id2,idb,nameb(im1),names(im1),names(im2),
     &   body(im1)*zmbar,body(im2)*zmbar,sm1s,sm2s,body12,age1,age2
*
         names(im1) = 0
         names(im2) = 0
*
*       change binary parameters in the case when one or two stars
*       change their masses because of stellar evolution (age1=age2)
*
         body(im1) = sm1s/zmbar
         body(im2) = sm2s/zmbar
         body12 = body(im1) + body(im2)
*
         bin(nbin3,3) = 0.5d0*body(im1)*body(im2)/vb
         bin(nbin3,4) = sqrt(ran2(irun))
         bin(nbin3,5) = vb
         bin(nbin3,6) = time
         bin(nbin3,7) = 1.0d8*float(id1) + float(id2)
*
         anb = bin(nbin3,3)*rbar/rtidkg/rsuntopc
*
      print*,'formb3-for im1,im2,id1,id2,idb,nb3,sm1,sm2,anb,ecc,tim=',
     &       im1,im2,id1,id2,idb,nbin3,sm1s,sm2s,anb,bin(nbin3,4),tim
*
*      create a new binary
*
         call create_binary(idb,id1,id2,anb,bin(nbin3,4),tim)
*
*       check if binary exist after formation
*
         call binary_exists(idb,iexist)
         if(iexist.eq.0) then
           print*,'formb3-noexist im1,idb,iexist',im1,idb,iexist
         endif
*
         call get_loid(idb,id1)
         call get_hiid(idb,id2)
         call get_loid_mass(idb,sm1so)
         call get_hiid_mass(idb,sm2so)
         call get_radius(id1,r1)
         call getLum(id1,xlum1)
         call get_radius(id2,r2)
         call getLum(id2,xlum2)
         call get_bs_type(idb,ik0n)
         call get_ss_type(id1,ik1n)
         call get_ss_type(id2,ik2n)
         ik0n = -1
         ik1 = -1
         ik2 = -1
*
         write(44,420) idb,id1,id2,im1,tim,sm1so,sm1s,sm2so,sm2s,
     &                 r1,r2,xlum1,xlum2,anb,bin(nbin3,4),r(k11),ik0,
     &                 ik0n,ik1,ik1n,ik2,ik2n,iexist,ibstra(im1),0
*
 420               format(1x,'#',4i7,1p12e16.8,9i4)
*
         bin(nbin3,1) = sm1so/zmbar
         bin(nbin3,2) = sm2so/zmbar
         bin(nbin3,7) = 1.0d8*float(id1) + float(id2)
c
      print*,'im1,im2,id1,id2,idb,ikind,sm1so,sm1s,sm2so,sm2,body(im1)'
     &,',body(im2),bin(nb,1),bin(nb,2),body12 =',im1,im2,id1,id2,idb,
     &ikind(im1),sm1so,sm1s,sm2so,sm2s,body(im1)*zmbar,body(im2)*zmbar,
     &sm1s,sm2s,body12*zmbar
*
*      store information about binary field star interaction
*
         iinte3(nbin3) = iinte3(nbin3) + 1
         inum = iinte3(nbin3)
*
         if(inum.gt.50) then
           iinte3(nbin3) = 1
           inum = 1
         endif
*
         r(k11) = (body(im1)*r(k11) + body(im2)*r(k2))/body12
         bin(nbin3,8) = r(k11)
*
         inumt = 50*(nbin3 - 1) + inum
         binin(inumt,1) = r(k11)/pctonbu
         binin(inumt,2) = vb
         binin(inumt,3) = body12*zmbar
         binin(inumt,4) = body(im3)*zmbar
         binin(inumt,5) = tim
         binin(inumt,6) = bin(nbin3,5)
         binin(inumt,7) = 0
*
*      determine energy generated during binary formation and split
*      this energy for the binary and the star
*
         print*,'vr(im1),vt(im1),vr(im2),vt(im2) = ',vr(im1),
     &          vt(im1),vr(im2),vt(im2)
         vb = 0.5d0*(body(im1)*(vr(im1)**2 + vt(im1)**2) +
     &        body(im2)*(vr(im2)**2 + vt(im2)**2) + smave*velave)
         body(im1) = body(im1) + body(im2)
c            vr(im1) = sqrt(smave*velave/(3.0d0*body12))
c            vt(im1) = sqrt(2.0d0*smave*velave/(3.0d0*body12))
         vr(im1) = sqrt(smave*velave/(3.0d0*body(im1)))
         vt(im1) = sqrt(2.0d0*smave*velave/(3.0d0*body(im1)))
c            nameb(im1) = nbin3
         r(k2) = 1.0d8 + float(im2)
         xmax(im2) = 1.01d0*r(k2)
*
         print*,'vb,smave,velave,body(im1),body(im2),body(im3),',
     &   'vr(im1),vt(im1),vr(im2),vt(im2) = ',vb,smave,velave,
     &   body(im1)*zmbar,body(im2)*zmbar,body(im3)*zmbar,vr(im1),
     &   vt(im1),vr(im2),vt(im2)
            vb = vb*body(im3)/(body(im1) + body(im3))
         call kick(k11,vb)
         write(6,*) 'kick-b-k1 k1,im1,debb,m12=',k11,im1,vb,
     &               body(im1)*zmbar
         vb = vb*body(im1)/body(im3)
         call kick(k3,vb)
         write(6,*) 'kick-s-k3  k3,im3,deb3,m3=',k3,im3,vb,
     &              body(im3)*zmbar
*
         ehbin3 = ehbin3 + bin(nbin3,5)
*
         go to 20
*
 10   continue
*
*      compute new potential and sort stars after removing from the
*      calculation star which is now one of the binary components (centre
*      of mass of newly formed binary is substitute for position of
*      another star which form binary)
*
      close(43)
      close(44)
*
      if(ibin.eq.0.and.noibin.eq.0) return
*
      call sort3(n,r,iname,inameo)
      nt = nt - ibin
      smt = smt
      nzst(nsuzon) = nt
      call coepot
      call energy(2)
      ptnew = pot
      zknew = zkin
      ehmlevn = ehmlev
      ehbin3n = ehbin3
      ekicktn = ekickt
      enepot = enepot + ptnew - ptold
      write(6,5432) zkold,zknew,ehmlevo,ehmlevn,ehbin3o,ehbin3n,
     &              ekickto,ekicktn,zknew-zkold,ehmlevn-ehmlevo,
     &              ehbin3n-ehbin3o,ekicktn-ekickto,(zknew-zkold)-
     &              (ehmlevn-ehmlevo+ehbin3n-ehbin3o+ekicktn-ekickto)
 5432 format(1x,'formb3  zkold,zknew,ehmlevo,ehmlevn,ekickto,ekicktn',
     &',ehbin3o,ehbin3n,Dzk,Dhm,Dhb3,Dhki,DDform =',1p13e14.6)
*
      etotn = zkin - pot + escsta - ehbin3 + enepot - ehb3b3 -
     &         ehmlev - ehcoll + eccoll - ekickt - enekin
      write(6,6543) etotn,zkin,pot,ehmlev,enepot,ehcoll,eccoll,
     &               ekickt,enekin,escsta,ehbin3,ehb3b3
 6543 format(1x,'formb3-f etotn,zkin,pot,ehmlev,enepot,ehcoll,',
     &'eccol,ekickt,enekin,escsta,ehbin3,ehb3b3 =',1p12e20.12)
*
      enekin = enekin + zknew - zkold
*
      print*,'formb3-fin  ibin,noibin,zko,zkn,poto,potn, = ',
     &  ibin,noibin,zkold,zknew,ptold,ptnew,enepot,enekin
      enekin = enekin - zknew + zkold
*
      return
      end
*
*

      logical function isnan(x)
      double precision x
      if (.not.x.le.0.d0.and..not.x.gt.0.d0) then
         isnan = .true.
      else
         isnan = .false.
      endif
      return
      end

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
   10 WRITE(*,*) 'SPAM'
      GO TO 10
      END

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
      subroutine print_mapinfo(lulog,map_info,n12)
      !
      ! print the content of the mysterious mapping information
      !
      implicit none

      integer, intent(in) ::
     &     lulog, n12, map_info(*)

      integer ::
     &     idx_minf, idx12, idx1, idx2, nsplit

      idx_minf = 0
      idx12 = 0
      do while(idx12.lt.n12)
        idx12 = idx12+1
        idx1 = 0
        idx2 = 0
        idx_minf = idx_minf+1
        nsplit = map_info(idx_minf)
        write(lulog,'("vtx",i4,": vtxs from op1: ",10i4)')
     &       idx12,map_info(idx_minf+1:idx_minf+nsplit)
     &
        idx_minf = idx_minf+nsplit

        idx_minf = idx_minf+1
        nsplit = map_info(idx_minf)
        write(lulog,'(3x,4x,2x,"vtxs from op2: ",10i4)')
     &       idx12,map_info(idx_minf+1:idx_minf+nsplit)
        idx_minf = idx_minf+nsplit
      end do

      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
c=======================================================================
      subroutine ptsurv
     &                (nrec,nsubject,nvar,p,maxni,datastr,
     &                 type,left,right,
     &                 x,subject,options,
     &                 priorv,prec,sb,tinvb,tauv,
     &                 mcmc,nsave,
     &                 acrate,randsave,thetasave,
     &                 cpar,beta,b,sigmab,sigmav,y,
     &                 iflagp,iflagr,parti,poss,typeint,whicho,whichn,
     &                 betac,linf,lsup,linf2,lsup2,
     &                 sigmainvb,sigmainvv,sigmainvvc,
     &                 sigmabc,sigmavc,ssb,ssbc,theta,thetac,vz,vzc,
     &                 workmr,workmr2,workmhp,
     &                 workmhr,workmhr2,workvp,workvr,
     &                 xtx,xty,yvec,yvecc,ztz,zty,
     &                 seed)

c=======================================================================
c     # 65 arguments
c
c     Subroutine `ptsurv' to run a Markov chain in a semiparametric
c     AFT effect model, using a Mixture of Multivariate Polya
c     trees prior for the distribution of the errors and the frailty
c     terms.
c
c     Copyright: Alejandro Jara and Tim Hanson, 2009-2010.
c
c     Version 1.0:
c
c     Last modification: 20-06-2007.
c
c      Alejandro Jara
c      Department of Statistics
c      Facultad de Matematicas
c      Pontificia Universidad Catolica de Chile
c      Casilla 306, Correo 22
c      Santiago
c      Chile
c      Voice: +56-2-3544506  URL  : http://www.mat.puc.cl/~ajara
c      Fax  : +56-2-3547729  Email: atjara@uc.cl
c
c=======================================================================
      implicit none

c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c++++ Data
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer maxni,nrec,nsubject,nvar,p
      integer datastr(nsubject,maxni+1)
      integer type(nrec,nvar)
      integer subject(nrec)
      real*8 left(nrec,nvar),right(nrec,nvar)
      real*8 x(nrec,p)

c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c++++ Options
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer frailty,latent,options(2)

c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+++++Prior
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer m1,m2
      real*8 prec(p*nvar,p*nvar),sb(p*nvar,2)
      real*8 nu0b,tinvb(nvar,nvar)
      real*8 tauv(nvar,2)
      real*8 priorv(3)

c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+++++MCMC
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer mcmc(5),nburn,nskip,nsave,ndisplay
      real*8 tune1,tune2

c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+++++Output
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      real*8 acrate(5)
      real*8 randsave(nsave,nsubject*nvar+2*nvar)
      real*8 thetasave(nsave,p*nvar+nvar+(nvar*(nvar+1))/2+2)

c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+++++Current values of the parameters
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      real*8 cpar(2),cparb,cparv,beta(p*nvar),b(nsubject,nvar)
      real*8 sigmab(nvar,nvar),sigmav(nvar,nvar)
      real*8 y(nrec,nvar)

c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c++++ External working space
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer iflagp(p*nvar)
      integer iflagr(nvar)
      integer parti(nvar)
      integer poss(nvar,p)
      integer typeint(nvar)
      integer whicho(nrec),whichn(nrec)

      real*8 betac(p*nvar)
      real*8 linf(nvar),lsup(nvar)
      real*8 linf2(nvar),lsup2(nvar)
      real*8 sigmainvb(nvar,nvar)
      real*8 sigmainvv(nvar,nvar)
      real*8 sigmainvvc(nvar,nvar)
      real*8 sigmabc(nvar,nvar)
      real*8 sigmavc(nvar,nvar)
      real*8 ssb(nsubject)
      real*8 ssbc(nsubject)
      real*8 theta(nvar)
      real*8 thetac(nvar)
      real*8 vz(nrec,nvar),vzc(nrec,nvar)
      real*8 workmr(nvar,nvar),workmr2(nvar,nvar)
      real*8 workmhp(p*nvar*(p*nvar+1)/2)
      real*8 workmhr(nvar*(nvar+1)/2),workmhr2(nvar*(nvar+1)/2)
      real*8 workvp(p*nvar)
      real*8 workvr(nvar)
      real*8 xtx(p*nvar,p*nvar)
      real*8 xty(p*nvar)
      real*8 yvec(nvar),yvecc(nvar)
      real*8 ztz(nvar,nvar)
      real*8 zty(nvar)

c+++++++++++++++++++++++++++++++++++++++++++++++++
c+++++Working space - RNG
c+++++++++++++++++++++++++++++++++++++++++++++++++
      integer seed1,seed2,seed(2)

c+++++++++++++++++++++++++++++++++++++++++++++++++
c++++ Internal working space
c+++++++++++++++++++++++++++++++++++++++++++++++++
      integer count
      integer dispcount
      integer i,j,k,l
      integer i1,j1
      integer ii,jj,kk,ll,mm
      integer iscan,isave
      integer marea
      parameter (marea=2**20)
      integer massi(marea)
      integer narea
      integer ni
      integer nscan
      integer nu
      integer skipcount
      integer sprint

      real*8 acrate2
      real*8 detlogb,detlogbc,detlogv,detlogvc
      real*8 logcgkn,logcgko
      real*8 loglikn,loglikn2,logliko,logliko2
      real*8 logpriorn,logprioro
      real*8 mass(marea)
      real*8 ratio
      real*8 sd,sdc
      real*8 ssevaln,ssevalo
      real*8 tmp1,tmp2

      integer maxp
      parameter(maxp=200)
      real*8 mumh(maxp),sigmamh(maxp,maxp)

      real*8 dlnrm,rtlnorm
      real runif

c+++++++++++++++++++++++++++++++++++++++++++++++++
c+++++CPU time
c+++++++++++++++++++++++++++++++++++++++++++++++++
      real*8 sec00,sec0,sec1,sec

c+++++++++++++++++++++++++++++++++++++++++++++++++
c++++ Initialize variables
c+++++++++++++++++++++++++++++++++++++++++++++++++

c+++++checking dimensions
      narea=2**nvar
      if(narea.gt.marea)then
         call rexit("increase dimension marea in subroutine ptsurv")
      end if

      if(nvar*p.gt.maxp)then
         call rexit("increase dimension maxp in subroutine ptsurv")
      end if

      do i=1,marea
         massi(i)=0
         mass(i)=0.d0
      end do

      do i=1,maxp
         do j=1,maxp
            sigmamh(i,j)=0.d0
         end do
         mumh(i)=0.d0
      end do

c++++ options
      frailty=options(1)
      latent=options(2)

c++++ mcmc, priors and "zipped"

      nburn=mcmc(1)
      nskip=mcmc(2)
      ndisplay=mcmc(3)
      m1=mcmc(4)
      m2=mcmc(5)

      nu0b=priorv(1)
      tune1=priorv(2)
      tune2=priorv(3)

      seed1=seed(1)
      seed2=seed(2)

      cparv=cpar(1)
      cparb=cpar(2)

c++++ set random number generator

      call setall(seed1,seed2)

c++++ set possition of covariates

      count=0
      do i=1,nvar
         do j=1,p
            count=count+1
            poss(i,j)=count
         end do
      end do

c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c++++ start the MCMC algorithm
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      isave=0
      skipcount=0
      dispcount=0
      nscan=nburn+(nskip+1)*nsave

      do i=1,nvar*p
         mumh(i)=beta(i)
         do j=1,nvar*p
            sigmamh(i,j)=0.d0
         end do
         sigmamh(i,i)=1.d0
      end do

      call cpu_time(sec0)
      sec00=0.d0

c++++ first evaluation of the likelihood

      detlogv=0.d0
      do i=1,nvar
         sigmainvv(i,i)=1.d0/sigmav(i,i)
         detlogv=detlogv+log(sigmav(i,i))
      end do

      call loglikpt_betacan_mz(m1,nvar,nrec,nsubject,subject,p,parti,
     &                         whicho,whichn,y,vz,cparv,detlogv,
     &                         linf,lsup,theta,sigmainvv,
     &                         x,beta,poss,b,logliko)

c      call dblepr("logliko",-1,logliko,1)
c      call dblepr("beta",-1,beta,p*nvar)
c      call dblepr("sigmav",-1,sigmav,nvar*nvar)
c      call intpr("frailty",-1,frailty,1)
c      call intpr("m1",-1,m1,1)
c      call intpr("m2",-1,m2,1)

      if(frailty.eq.1)then
         do i=1,nvar
            do j=1,nvar
               sigmainvb(i,j)=sigmab(i,j)
            end do
         end do
         call inversedet(sigmainvb,nvar,iflagr,detlogb)
      end if

      do iscan=1,nscan

c+++++++ check if the user has requested an interrupt
         call rchkusr()

c         call intpr("Step1",-1,latent,1)

c++++++++++++++++++++++++++++++++
c+++++++ latent variables     +++
c++++++++++++++++++++++++++++++++
         if(latent.eq.1)then

         acrate2=0.d0

         do i=1,nrec

c++++++++++ check if the user has requested an interrupt
            call rchkusr()

            do j=1,nvar
               yvec(j)=y(i,j)
               linf(j)=left(i,j)
               lsup(j)=right(i,j)
               linf2(j)=left(i,j)
               lsup2(j)=right(i,j)
               typeint(j)=type(i,j)
            end do

            do j=1,nvar
               tmp1=0.d0
               do k=1,p
                  l=poss(j,k)
                  tmp1=tmp1-x(i,k)*beta(l)
               end do
               theta(j)=tmp1+b(subject(i),j)
            end do

c++++++++++ generating the candidate value

            call rtmvn(nvar,theta,sigmav,linf,lsup,workmr,workmr2,
     &                 typeint,yvecc,logcgko)

            call dtmvn(nvar,theta,sigmav,linf,lsup,workmr,workmr2,
     &                 typeint,yvecc,yvec,logcgkn)

c++++++++++ evaluating likelihood

            call loglikpt_ycan_mz(i,m1,nvar,nrec,nsubject,subject,
     &                            p,parti,whicho,whichn,yvec,zty,
     &                            vz,cparv,detlogv,linf,lsup,theta,
     &                            sigmainvv,x,beta,poss,b,logliko)

            call loglikpt_ycan_mz(i,m1,nvar,nrec,nsubject,subject,
     &                            p,parti,whicho,whichn,yvecc,zty,
     &                            vz,cparv,detlogv,linf,lsup,theta,
     &                            sigmainvv,x,beta,poss,b,loglikn)

c++++++++++ mh step

            ratio=loglikn-logliko+logcgkn-logcgko

            if(log(dble(runif())).lt.ratio)then
               acrate2=acrate2+1.d0
               do j=1,nvar
                  y(i,j)=yvecc(j)
                  vz(i,j)=zty(j)
               end do
            end if
         end do
         acrate(1)=acrate(1)+acrate2/dble(nrec)

c+++++++ evaluating the current value of the likelihood

         call loglikpt_betacan_mz(m1,nvar,nrec,nsubject,subject,p,parti,
     &                         whicho,whichn,y,vz,cparv,detlogv,
     &                         linf,lsup,theta,sigmainvv,
     &                         x,beta,poss,b,logliko)

         end if

c         call dblepr("y",-1,y,nsubject*nvar)

c++++++++++++++++++++++++++++++++
c+++++++ fixed effects        +++
c++++++++++++++++++++++++++++++++

         do i=1,p*nvar
            do j=1,p*nvar
               xtx(i,j)=prec(i,j)
            end do
            xty(i)=sb(i,1)
         end do

         do ii=1,nrec
            do jj=1,nvar
               yvec(jj)=y(ii,jj)-b(subject(ii),jj)
            end do

            do jj=1,nvar
               do kk=1,p
                  i1=poss(jj,kk)
                  do ll=1,nvar
                     do mm=1,p
                        j1=poss(ll,mm)
                        xtx(i1,j1)=xtx(i1,j1)+sigmainvv(jj,ll)*
     &                             x(ii,kk)*x(ii,mm)
                     end do
                     xty(i1)=xty(i1)-sigmainvv(jj,ll)*yvec(ll)*x(ii,kk)
                  end do
               end do
            end do
         end do
         call inverse(xtx,p*nvar,iflagp)

         if(iscan.lt.nburn/2)then
            do i=1,p*nvar
               do j=1,p*nvar
                  xtx(i,j)=0.01d0*xtx(i,j)
               end do
            end do
            call rmvnorm(p*nvar,beta,xtx,workmhp,workvp,betac)

          else
            ratio=dble(runif())
            if(ratio.le.0.25d0)then
               do i=1,p*nvar
                  do j=1,p*nvar
                      xtx(i,j)=(5.4264d0/dble(p*nvar))*sigmamh(i,j)
                  end do
               end do
               call rmvnorm(p*nvar,beta,xtx,workmhp,workvp,betac)
             else if(ratio.le.0.5d0)then
               do i=1,p*nvar
                  do j=1,p*nvar
                      xtx(i,j)=(0.001d0)*xtx(i,j)
                  end do
               end do
               call rmvnorm(p*nvar,beta,xtx,workmhp,workvp,betac)
             else if(ratio.le.0.75d0)then
               do i=1,p*nvar
                  do j=1,p*nvar
                      xtx(i,j)=(0.001d0)*xtx(i,j)
                  end do
               end do
               call rmvnorm(p*nvar,beta,xtx,workmhp,workvp,betac)
             else
               do i=1,p*nvar
                  do j=1,p*nvar
                     xtx(i,j)=xtx(i,j)
                  end do
               end do
               call rmvnorm(p*nvar,beta,xtx,workmhp,workvp,betac)
            end if
         end if

c+++++++ prior ratio
         logprioro=0.d0
         logpriorn=0.d0

         do i=1,p*nvar
            do j=1,p*nvar
               logpriorn=logpriorn+(betac(i)-sb(i,2))*
     &                   prec(i,j)       *
     &                  (betac(j)-sb(j,2))

               logprioro=logprioro+(beta(i)-sb(i,2))*
     &                   prec(i,j)      *
     &                   (beta(j)-sb(j,2))

            end do
         end do

         logpriorn=-0.5d0*logpriorn
         logprioro=-0.5d0*logprioro

c+++++++ evaluating likelihood for betac

         call loglikpt_betacan_mz(m1,nvar,nrec,nsubject,subject,p,parti,
     &                         whicho,whichn,y,vzc,cparv,detlogv,
     &                         linf,lsup,theta,sigmainvv,
     &                         x,betac,poss,b,loglikn)

c         call dblepr("logliko",-1,logliko,1)
c         call dblepr("loglikn",-1,loglikn,1)

c+++++++ acceptance step

         ratio=loglikn-logliko+logpriorn-logprioro

         if(log(dble(runif())).lt.ratio)then
            acrate(2)=acrate(2)+1.d0
            do i=1,p*nvar
               beta(i)=betac(i)
            end do
            logliko=loglikn
            do i=1,nrec
               do j=1,nvar
                  vz(i,j)=vzc(i,j)
               end do
            end do
         end if

c+++++++ adapting the parameters for the MH algorithm

         do i=1,p*nvar
            do j=1,p*nvar
               tmp1=(beta(i)-mumh(i))*(beta(j)-mumh(j))

               if(tmp1.eq.0.d0)then
                  tmp2=sigmamh(i,j)/dble(iscan)
                else
                  tmp2=
     &            (tmp1-
     &             (dble(iscan+1)/dble(iscan))*sigmamh(i,j)
     &            )/dble(iscan+1)
               end if
               sigmamh(i,j)=sigmamh(i,j)+tmp2
            end do
         end do

         do i=1,p*nvar
            tmp1=(beta(i)-mumh(i))/dble(iscan+1)
            mumh(i)=mumh(i)+tmp1
         end do

         call dblepr("beta",-1,beta,p*nvar)

c         call dblepr("betac",-1,betac,p*nvar)
c         call dblepr("mumh",-1,mumh,p*nvar)
c         do i=1,p*nvar
c            do j=i,p*nvar
c                call dblepr("sigmamh",-1,sigmamh(i,j),1)
c            end do
c         end do

c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+++++++ Updating sigmav using a MH step                    +++
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

c+++++++ check if the user has requested an interrupt
         call rchkusr()

c+++++++ generating the candidate value/
c+++++++ evaluating the candidate generating kernel/
c+++++++ evaluate log-prior

         logcgkn=0.d0
         logcgko=0.d0
         detlogvc=0.d0
         logpriorn=0.d0
         logprioro=0.d0

         do i=1,nvar
            sd=sigmainvv(i,i)
            sdc=rtlnorm(log(sd),tune1*0.01d0,0,0,.true.,.true.)
            sigmainvvc(i,i)=sdc

            sigmavc(i,i)=1.d0/sigmainvvc(i,i)
            detlogvc=detlogvc+log(sigmavc(i,i))

            logcgkn=logcgkn+dlnrm(sd,log(sdc),tune1*0.01d0,1)
            logcgko=logcgko+dlnrm(sdc,log(sd ),tune1*0.01d0,1)

            call dgamma2(sigmainvvc(i,i),0.5d0*tauv(i,1),
     &                   0.5d0*tauv(i,2),tmp1)
            logpriorn=logpriorn+tmp1

            call dgamma2(sigmainvv(i,i) ,0.5d0*tauv(i,1),
     &                   0.5d0*tauv(i,2),tmp1)
            logprioro=logprioro+tmp1
         end do

c+++++++ evaluating likelihood

         call loglikpt_betacan_mz(m1,nvar,nrec,nsubject,subject,p,parti,
     &                            whicho,whichn,y,vzc,cparv,detlogvc,
     &                            linf,lsup,yvec,sigmainvvc,
     &                            x,beta,poss,b,loglikn)

c+++++++ acceptance step

         ratio=loglikn-logliko+logcgkn-logcgko+
     &         logpriorn-logprioro

         if(log(dble(runif())).lt.ratio)then
            do i=1,nvar
               sigmav(i,i)=sigmavc(i,i)
               sigmainvv(i,i)=sigmainvvc(i,i)
            end do
            detlogv=detlogvc
            do i=1,nrec
               do j=1,nvar
                  vz(i,j)=vzc(i,j)
               end do
            end do
            logliko=loglikn
            acrate(3)=acrate(3)+1.d0
         end if

         call dblepr("sigmav",-1,sigmav,nvar*nvar)

c++++++++++++++++++++++++++++++++
c+++++++ frailty terms        +++
c++++++++++++++++++++++++++++++++

         if(frailty.eq.1)then
            acrate2=0.d0
            do i=1,nsubject

c+++++++++++++ check if the user has requested an interrupt
               call rchkusr()

               do j=1,nvar
                  theta(j)=b(i,j)
               end do

c+++++++++++++ generating a candidate
               do j=1,nvar
                  tmp1=0.d0
                  do k=1,nvar
                     ztz(j,k)=sigmainvb(j,k)
                  end do
                  zty(j)=0.d0
               end do

               ni=datastr(i,1)

               do j=1,ni
                  do k=1,nvar
                     tmp1=0.d0
                     do l=1,p
                        jj=poss(j,k)
                        tmp1=tmp1-x(datastr(i,j+1),l)*beta(jj)
                     end do
                     thetac(k)=tmp1
                  end do

                  do k=1,nvar
                     ztz(k,k)=ztz(k,k)+sigmainvv(k,k)
                     zty(k)=zty(k)+(y(datastr(i,j+1),k)-thetac(k))*
     &                      sigmainvv(k,k)
                  end do
               end do

               call inverse(ztz,nvar,iflagr)

               call rmvnorm(nvar,theta,ztz,workmhr,workvr,thetac)

c+++++++++++++ evaluating the prior

               logprioro=0.d0
               logpriorn=0.d0

               call condptpriorell(i,m2,nrec,nsubject,nvar,theta,ssb,
     &                          sigmainvb,detlogb,cparb,
     &                          whicho,whichn,
     &                          logprioro,ssevalo)

               call condptpriorell(i,m2,nrec,nsubject,nvar,thetac,ssb,
     &                          sigmainvb,detlogb,cparb,
     &                          whicho,whichn,
     &                          logpriorn,ssevaln)

c+++++++++++++ evaluating likelihood

               do j=1,nvar
                  b(i,j)=theta(j)
               end do

               call loglikpt_bcan_mz(m1,nvar,nrec,nsubject,subject,p,
     &                         parti,whicho,whichn,y,vzc,cparv,detlogv,
     &                         linf,lsup,yvec,sigmainvv,
     &                         x,beta,poss,b,logliko)

               do j=1,nvar
                  b(i,j)=thetac(j)
               end do

               call loglikpt_bcan_mz(m1,nvar,nrec,nsubject,subject,p,
     &                         parti,whicho,whichn,y,vzc,cparv,detlogv,
     &                         linf,lsup,yvec,sigmainvv,
     &                         x,beta,poss,b,loglikn)

c+++++++++++++ mh step

               ratio=loglikn-logliko+
     &               logcgkn-logcgko+
     &               logpriorn-logprioro

               if(log(dble(runif())).lt.ratio)then
                  acrate2=acrate2+1.d0
                  do j=1,nvar
                     b(i,j)=thetac(j)
                     vz(i,j)=vzc(i,j)
                  end do
                  ssb(i)=ssevaln
                  logliko=loglikn
                 else
                  do j=1,nvar
                     b(i,j)=theta(j)
                  end do
               end if
            end do
            acrate(4)=acrate(4)+acrate2/dble(nsubject)
         end if

c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+++++++ Updating sigmab using a MH step                    +++
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

         if(frailty.eq.1)then

c++++++++++ check if the user has requested an interrupt
            call rchkusr()

c++++++++++ generating the candidate value

            nu=int((dble(nsubject))*tune2)

            do i=1,nvar
               do j=1,nvar
                  sigmabc(i,j)=dble(nu-nvar-1)*sigmab(i,j)
               end do
            end do

            call riwishart(nvar,nu,sigmabc,sigmainvvc,workmr,workvr,
     &                     workmhr,workmhr2,iflagr)

c++++++++++ evaluating the candidate generating kernel

            do i=1,nvar
               do j=1,nvar
                  ztz(i,j)=dble(nu-nvar-1)*sigmab(i,j)
               end do
            end do

            call diwishart(nvar,nu,sigmabc,ztz,workmr,workmr2,workvr,
     &                  iflagr,logcgko)

            do i=1,nvar
               do j=1,nvar
                  ztz(i,j)=dble(nu-nvar-1)*sigmabc(i,j)
               end do
            end do

            call diwishart(nvar,nu,sigmab,ztz,workmr,workmr2,workvr,
     &                  iflagr,logcgkn)

c++++++++++ evaluating the prior

            call diwishart(nvar,int(nu0b),sigmabc,tinvb,workmr,
     &                  workmr2,workvr,iflagr,logpriorn)

            call diwishart(nvar,int(nu0b),sigmab,tinvb,workmr,
     &                  workmr2,workvr,iflagr,logprioro)

c++++++++++ evaluating likelihood

            call loglik_ell(m2,nrec,nsubject,nvar,b,ssb,
     &                      sigmainvb,detlogb,cparb,
     &                      whicho,whichn,
     &                      logliko2)

            do i=1,nvar
               do j=1,nvar
                  sigmainvvc(i,j)=sigmabc(i,j)
               end do
            end do

            call inversedet(sigmainvvc,nvar,iflagr,detlogbc)

            call loglik_ell(m2,nrec,nsubject,nvar,b,ssbc,
     &                      sigmainvvc,detlogbc,cparb,
     &                      whicho,whichn,
     &                      loglikn2)

c++++++++++ acceptance step

            ratio=loglikn2-logliko2+logcgkn-logcgko+
     &            logpriorn-logprioro

            if(log(dble(runif())).lt.ratio)then
               do i=1,nvar
                  do j=1,nvar
                     sigmab(i,j)=sigmabc(i,j)
                     sigmainvb(i,j)=sigmainvvc(i,j)
                  end do
               end do
               detlogb=detlogbc
               do i=1,nsubject
                  ssb(i)=ssbc(i)
               end do
               acrate(5)=acrate(5)+1.d0
            end if

         end if

c++++++++++++++++++++++++++++++++++
c+++++++ save samples
c++++++++++++++++++++++++++++++++++

         cpar(1)=cparv
         cpar(2)=cparb

         if(iscan.gt.nburn)then
            skipcount=skipcount+1
            if(skipcount.gt.nskip)then
               isave=isave+1
               dispcount=dispcount+1

c+++++++++++++ regression coefficient information

               do i=1,p*nvar
                  thetasave(isave,i)=beta(i)
               end do

c+++++++++++++ baseline covariance for the errors

               do i=1,nvar
                  thetasave(isave,p*nvar+i)=sigmav(i,i)
               end do

c+++++++++++++ baseline covariance for the frailty

               if(frailty.eq.1)then
                  k=0
                  do i=1,nvar
                     do j=i,nvar
                        k=k+1
                        thetasave(isave,p*nvar+nvar+k)=sigmab(i,j)
                     end do
                  end do
                else
                  k=nvar*(nvar+1)/2
               end if

c+++++++++++++ precision parameter information

               thetasave(isave,p*nvar+nvar+k+1)=cparv
               if(frailty.eq.1)then
                  thetasave(isave,p*nvar+nvar+k+2)=cparb
               end if

c+++++++++++++ frailty

               if(frailty.eq.1)then
                  k=0
                  do i=1,nsubject
                     do j=1,nvar
                        k=k+1
                        randsave(isave,k)=b(i,j)
                     end do
                  end do
                 else
                  k=nsubject*nvar
               end if

c+++++++++++++ error predictive information

               call sampredptun(marea,nvar,nrec,parti,m1,mass,massi,
     &                          iflagr,typeint,whichn,whicho,vz,
     &                          cparv,thetac,linf,lsup,workmr,
     &                          workmhr,workvr,theta,1)

               do i=1,nvar
                  tmp1=0.d0
                  theta(i)=sqrt(sigmav(i,i))*theta(i)
                  k=k+1
                  randsave(isave,k)=theta(i)
               end do

c+++++++++++++ frailty predictive information
               if(frailty.eq.1)then
                  do i=1,nvar
                     thetac(i)=0.d0
                  end do

                  call sampredellpt(nrec,nsubject,nvar,m2,ssb,cparb,
     &                              whicho,whichn,workmr,workmhr,workvr,
     &                              thetac,sigmab,theta)
                  do i=1,nvar
                     k=k+1
                     randsave(isave,k)=theta(i)
                  end do
               end if

c+++++++++++++ print
               skipcount = 0
               if(dispcount.ge.ndisplay)then
                  call cpu_time(sec1)
                  sec00=sec00+(sec1-sec0)
                  sec=sec00
                  sec0=sec1
                  tmp1=sprint(isave,nsave,sec)
                  dispcount=0
               end if
            end if
         end if

      end do

      do i=1,5
         acrate(i)=acrate(i)/dble(nscan)
      end do

      return

      end

c=======================================================================
      subroutine loglikpt_ycan_mz(ind,m,nvar,nrec,nsubject,subject,
     &                            p,parti,
     &                            whicho,whichn,yc,tvec,vz,cpar,detlogl,
     &                            linf,lsup,theta,sigmainv,
     &                            x,beta,poss,b,loglikc)
c=======================================================================
c     This subroutine evaluate the log-likelihood for the candidate
c     value of the survival time in a marginal Multivariate PT.
c
c     Alejandro Jara, 2006
c=======================================================================
      implicit none

c-----Input
      integer ind,m,nvar,nrec,nsubject,p
      integer parti(nvar),poss(nvar,p)
      integer subject(nrec)
      integer whicho(nrec),whichn(nrec)
      real*8 beta(nvar*p),b(nsubject,nvar),x(nrec,p)
      real*8 yc(nvar),tvec(nvar),vz(nrec,nvar),cpar,detlogl
      real*8 linf(nvar),lsup(nvar)
      real*8 theta(nvar),sigmainv(nvar,nvar)

c-----Output
      real*8 loglikc

c-----Working
      integer countero,countern,final
      integer i,j,je2,k,k1,k2,l,nint,ok
      real*8 dnrm,invcdfnorm,prob,quan,tmp1

c-----Routine

      loglikc=0.d0

      do i=1,nvar
         tmp1=0.d0
         do j=1,p
            k=poss(i,j)
            tmp1=tmp1-x(ind,j)*beta(k)
         end do
         theta(i)=tmp1+b(subject(ind),i)
      end do

      do i=1,nvar
         tmp1=0.d0
         do j=1,nvar
            tmp1=tmp1+sigmainv(i,j)*(yc(j)-theta(j))
         end do
         tvec(i)=tmp1
      end do

      nint=2
      prob=1.d0/dble(nint)
      quan=invcdfnorm(prob,0.d0,1.d0,1,0)

      countero=0

      do j=1,nvar
         if(tvec(j).le.quan)then
            linf(j)=-999999.d0
            lsup(j)=quan
            parti(j)=1
          else
            linf(j)=quan
            lsup(j)= 999999.d0
            parti(j)=2
         end if
      end do

      do l=1,nrec
         if(l.ne.ind)then
            final=1
            do j=1,nvar
               if(vz(l,j).gt.lsup(j).or.vz(l,j).lt.linf(j))then
                  final=0
               end if
c               if(y(l,j).gt.lsup2(j).or.y(l,j).lt.linf2(j))then
c                  final=0
c               end if
            end do

            if(final.eq.1)then
               countero=countero+1
               whicho(countero)=l
            end if
         end if
      end do

      if(countero.eq.0) go to 1

      ok=1
      j=2
      do while(ok.eq.1.and.j.le.m)
         nint=2**j
         je2=j**2
         prob=1.d0/dble(nint)

         do k=1,nvar
            k1=2*(parti(k)-1)+1
            k2=2*(parti(k)-1)+2
            quan=invcdfnorm(dble(k1)*prob,0.d0,1.d0,1,0)

            if(tvec(k).le.quan)then
              parti(k)=k1
              lsup(k)=quan
             else
              parti(k)=k2
              linf(k)=quan
            end if
         end do

         countern=0
         do l=1,countero
            final=1
            do k=1,nvar
               if(vz(whicho(l),k).gt.lsup(k).or.
     &            vz(whicho(l),k).lt.linf(k)    )then
                  final=0
               end if
            end do

            if(final.eq.1)then
              countern=countern+1
              whichn(countern)=whicho(l)
            end if
         end do

         loglikc=loglikc+
     &           log((2.d0**nvar)*cpar*dble(je2)+
     &               dble((2.d0**nvar)*countero))-
     &           log((2.d0**nvar)*cpar*dble(je2)+dble(i-1))

         if(countern.eq.0)then
            ok=0
          else
            countero=countern
            do l=1,countern
               whicho(l)=whichn(l)
            end do
            j=j+1
         end if
      end do

1     continue

      loglikc=loglikc-0.5d0*detlogl
      do j=1,nvar
         loglikc=loglikc+dnrm(tvec(j),0.d0, 1.d0, 1)
      end do

      return
      end

c=======================================================================
      subroutine loglikpt_betacan_mz(m,nvar,nrec,nsubject,subject,
     &                               p,parti,whicho,whichn,y,vzc,cpar,
     &                               detlogl,linf,lsup,theta,sigmainv,
     &                               x,betac,poss,b,loglikc)
c=======================================================================
c     This subroutine evaluate the log-likelihood for the candidate
c     value of the regression coefficients using a marginal
c     Multivariate PT in a median zero AFT regression model.
c
c     Alejandro Jara, 2006
c=======================================================================
      implicit none

c-----Input
      integer m,nvar,nrec,nsubject,p
      integer parti(nvar),poss(nvar,p)
      integer subject(nrec)
      integer whicho(nrec),whichn(nrec)
      real*8 betac(nvar*p),b(nsubject,nvar),x(nrec,p)
      real*8 y(nrec,nvar),vzc(nrec,nvar),cpar,detlogl
      real*8 linf(nvar),lsup(nvar)
      real*8 theta(nvar),sigmainv(nvar,nvar)

c-----Output
      real*8 loglikc

c-----Working
      integer countero,countern,final
      integer i,j,je2,k,k1,k2,l,nint,ok
      real*8 dnrm,invcdfnorm,prob,quan,tmp1

c-----Routine

      loglikc=0.d0

      do i=1,nrec
         do j=1,nvar
            tmp1=0.d0
            do k=1,p
               l=poss(j,k)
               tmp1=tmp1-x(i,k)*betac(l)
            end do
            theta(j)=tmp1+b(subject(i),j)
         end do

         do j=1,nvar
            tmp1=0.d0
            do k=1,nvar
               tmp1=tmp1+sigmainv(j,k)*(y(i,k)-theta(k))
            end do
            vzc(i,j)=tmp1
         end do

c+++++++ check if the user has requested an interrupt
         call rchkusr()

c+++++++ first subject
         if(i.eq.1)then
            loglikc=-0.5d0*detlogl
            do j=1,nvar
               loglikc=loglikc+dnrm(vzc(i,j),0.d0, 1.d0, 1)
            end do

c+++++++ following subjects
          else

            nint=2
            prob=1.d0/dble(nint)
            quan=invcdfnorm(prob,0.d0,1.d0,1,0)

            countero=0

            do j=1,nvar
               if(vzc(i,j).le.quan)then
                  linf(j)=-999999.d0
                  lsup(j)=quan
                  parti(j)=1
                else
                  linf(j)=quan
                  lsup(j)= 999999.d0
                  parti(j)=2
               end if
            end do

            do l=1,i-1
               final=1
               do j=1,nvar
                  if(vzc(l,j).gt.lsup(j).or.vzc(l,j).lt.linf(j))then
                    final=0
                  end if
               end do

               if(final.eq.1)then
                  countero=countero+1
                  whicho(countero)=l
               end if
            end do

            if(countero.eq.0) go to 1

            ok=1
            j=2
            do while(ok.eq.1.and.j.le.m)
               nint=2**j
               je2=j**2
               prob=1.d0/dble(nint)

               do k=1,nvar
                  k1=2*(parti(k)-1)+1
                  k2=2*(parti(k)-1)+2
                  quan=invcdfnorm(dble(k1)*prob,0.d0,1.d0,1,0)

                  if(vzc(i,k).le.quan)then
                    parti(k)=k1
                    lsup(k)=quan
                   else
                    parti(k)=k2
                    linf(k)=quan
                  end if
               end do

               countern=0
               do l=1,countero
                  final=1
                  do k=1,nvar
                     if(vzc(whicho(l),k).gt.lsup(k).or.
     &                  vzc(whicho(l),k).lt.linf(k)    )then
                        final=0
                     end if
                  end do

                  if(final.eq.1)then
                    countern=countern+1
                    whichn(countern)=whicho(l)
                  end if
               end do

               loglikc=loglikc+
     &            log((2.d0**nvar)*cpar*dble(je2)+
     &                 dble((2.d0**nvar)*countero))-
     &            log((2.d0**nvar)*cpar*dble(je2)+dble(i-1))

               if(countern.eq.0)then
                  ok=0
                else
                  countero=countern
                  do l=1,countern
                     whicho(l)=whichn(l)
                  end do
                  j=j+1
               end if
            end do

1           continue

            loglikc=loglikc-0.5d0*detlogl
            do j=1,nvar
               loglikc=loglikc+dnrm(vzc(i,j),0.d0, 1.d0, 1)
            end do
         end if
      end do

      return
      end

c=======================================================================
      subroutine loglikpt_bcan_mz(m,nvar,nrec,nsubject,subject,p,parti,
     &                            whicho,whichn,y,vzc,cpar,detlogl,
     &                            linf,lsup,muc,sigmainv,
     &                            x,beta,poss,b,loglikc)
c=======================================================================
c     This subroutine evaluate the log-likelihood for the candidate
c     value of the frailty term in a marginal Multivariate PT.
c
c     Alejandro Jara, 2006
c=======================================================================
      implicit none

c-----Input
      integer m,nvar,nrec,nsubject,p
      integer parti(nvar),poss(nvar,p)
      integer subject(nrec)
      integer whicho(nrec),whichn(nrec)
      real*8 beta(nvar*p),b(nsubject,nvar),x(nrec,p)
      real*8 y(nrec,nvar),vzc(nrec,nvar),cpar,detlogl
      real*8 linf(nvar),lsup(nvar)
      real*8 muc(nvar),sigmainv(nvar,nvar)

c-----Output
      real*8 loglikc

c-----Working
      integer countero,countern,final
      integer i,j,je2,k,k1,k2,l,nint,ok
      real*8 dnrm,invcdfnorm,prob,quan,tmp1

c-----Routine

      loglikc=0.d0

      do i=1,nrec

         do j=1,nvar
            tmp1=0.d0
            do k=1,p
               l=poss(j,k)
               tmp1=tmp1-x(i,k)*beta(l)
            end do
            muc(j)=tmp1+b(subject(i),j)
         end do

         do j=1,nvar
            tmp1=0.d0
            do k=1,nvar
               tmp1=tmp1+sigmainv(j,k)*(y(i,k)-muc(k))
            end do
            vzc(i,j)=tmp1
         end do

c+++++++ check if the user has requested an interrupt
         call rchkusr()

c+++++++ first subject
         if(i.eq.1)then
            loglikc=-0.5d0*detlogl
            do j=1,nvar
               loglikc=loglikc+dnrm(vzc(i,j),0.d0, 1.d0, 1)
            end do

c+++++++ following subjects
          else

            nint=2
            prob=1.d0/dble(nint)
            quan=invcdfnorm(prob,0.d0,1.d0,1,0)

            countero=0

            do j=1,nvar
               if(vzc(i,j).le.quan)then
                  linf(j)=-999999.d0
                  lsup(j)=quan
                  parti(j)=1
                else
                  linf(j)=quan
                  lsup(j)= 999999.d0
                  parti(j)=2
               end if
            end do

            do l=1,i-1
               final=1
               do j=1,nvar
                  if(vzc(l,j).gt.lsup(j).or.vzc(l,j).lt.linf(j))then
                    final=0
                  end if
               end do

               if(final.eq.1)then
                  countero=countero+1
                  whicho(countero)=l
               end if
            end do

            if(countero.eq.0) go to 1

            ok=1
            j=2
            do while(ok.eq.1.and.j.le.m)
               nint=2**j
               je2=j**2
               prob=1.d0/dble(nint)

               do k=1,nvar
                  k1=2*(parti(k)-1)+1
                  k2=2*(parti(k)-1)+2
                  quan=invcdfnorm(dble(k1)*prob,0.d0,1.d0,1,0)

                  if(vzc(i,k).le.quan)then
                    parti(k)=k1
                    lsup(k)=quan
                   else
                    parti(k)=k2
                    linf(k)=quan
                  end if
               end do

               countern=0
               do l=1,countero
                  final=1
                  do k=1,nvar
                     if(vzc(whicho(l),k).gt.lsup(k).or.
     &                  vzc(whicho(l),k).lt.linf(k)    )then
                        final=0
                     end if
                  end do

                  if(final.eq.1)then
                    countern=countern+1
                    whichn(countern)=whicho(l)
                  end if
               end do

               loglikc=loglikc+
     &            log((2.d0**nvar)*cpar*dble(je2)+
     &                dble((2.d0**nvar)*countero))-
     &            log((2.d0**nvar)*cpar*dble(je2)+dble(i-1))

               if(countern.eq.0)then
                  ok=0
                else
                  countero=countern
                  do l=1,countern
                     whicho(l)=whichn(l)
                  end do
                  j=j+1
               end if
            end do

1           continue

            loglikc=loglikc-0.5d0*detlogl
            do j=1,nvar
               loglikc=loglikc+dnrm(vzc(i,j),0.d0, 1.d0, 1)
            end do
         end if
      end do

      return
      end

max_stars_repo_name: kdevick/densitymediation

max_issues_repo_name: kdevick/densitymediation

max_forks_repo_name: kdevick/densitymediation


text:
no toad sexing in Davis
by Users/DanMasiel

{{{One fine evening, Larry and I were bored, broke and didnt want to do our homework. We decided it would be
interesting to attempt to penetrate the mysterious depths of the buildings on campus. Before class that day I
noticed an unlocked ventilation grate that seemed to lead to a tunnel. We set out from the Malcom dormitory
and headed towards Olson Hall. The unlocked grate outside seemed to lead to noting but a dark hole, so we
opted to explore the interior of the building instead.

Once inside we walked downstairs to the basement. We soon discovered a door that hadnt been closed completely.
Inside was a dimly lit concrete room full of pipes, wires, and boxes. At first the room didnt appear to be
that large, but after further inspection we discovered the entrance to a long tunnel at the back of the room.
We had set out, bored and not really expecting to find anything interesting. The sight of what appeared to be
a 100 meter long tunnel left us dumbfounded. The tunnel was approximately a meter wide and three meters tall
with pipes running along the ceiling. It was well lit and we immediately walked through it. On the other side
was a huge room. Half of the room was filled with boxes of old student records enclosed by a huge wire cage.
We wandered around the room for a while and found an elevator.

We got into the elevator and noticed that there were 9 floors in the building. Upon realizing that we were
in the basement of Sproul Hall, we were overcome with a feeling of immense satisfaction and a growing sense
that nothing could stop us on this adventure. Larry pressed the 9 button and we soon found ourselves in a
hallway full of offices. We wandered around a looked out the windows. We then entered the stairwell, and
before us stood something we could not refuse. A small metal staircase stretched out over the spiraling
stairwell and ended at a hatch to the roof. There appeared to be a lock hanging from the handle of the
hatch, but I decided to check it anyway. To my amazement the lock was not closed and the hatch was easily
opened.

Larry and I climbed out the hatch and emerged near a bunch of communications equipment on roof. The view
was incredible. Essentially all of Davis was laid out before us. After enjoying the view for a while we
climbed back through the hatch and walked back to the elevator. We pressed the basement button but
received no response from the machine. Thinking we had accomplished quite enough for one night, we took
the elevator to the bottom floor and walked victoriously back to the dorms.}}}

20041209 17:47:24 nbsp Awesome story.  Ive been through the SproulOlson underground passage, but I have yet to climb to the roof.  Users/MikeIvanov

20041209 22:10:54 nbsp Awesome,  I too have explored the tunnels under sproul  though Im prepared to deny that upon questioning by disciplinary authorities (lets apply the disclaimer in the emosnail userinfo to this).  You really lucked out with all the doors unlocked.  There also is a tunnel between Kerr Hall Kerr and Wellman Hall Wellman which Ive seen on floor schematics but never explored.  Users/KrisFricke

20041209 23:21:41 nbsp Kris:  The tunnel between Kerr and Wellman is far less exciting, simply a passage between the basements.  It is usually locked from the Wellman side, but you can always get from Kerr to Wellman.  Professors use this tunnel all the time in the rainy season.  Users/BrentLaabs

20041209 23:46:13 nbsp I wish I had a tunnel to campus to use when it rained  Users/KrisFricke

20041210 02:17:30 nbsp I hear that the tunnel from Kerr to Wellman is useful if you want to sneak into Wellman late at night to watch DVDs on the big screen with auditorium seating.  Of course, this is only good if you can get into Kerr.  Users/ZachSaul

20041210 07:36:27 nbsp And whats even better is, you could sneak into wellman late at night and play Halo 2 on the projection screens with several xboxes networked on the campus LAN  Users/DanMasiel

20060224 19:09:10 nbsp I havent seen the tunnel, but on my Birthday 2004, the same day Gary Snyder was retiring, I had a similar adventure on the roof of Sproul.  It is UCDs own Stairway to Heaven  Users/EricHolmes

20060518 23:10:17 nbsp Anyone know if this is still possible? Im pretty sure this trick doesnt work anymore, but when i was a freshman, Kerr had a trick to get in after hours. If you pressed the handicap button on the east entrance, the door didnt open, but it would unlock long enough to pull the door opened and walk in. If this still worked as well as always having the tunnel entrance opened from the kerr side, thatd be a good way to get into wellman after hours. Users/CarlosBarahona

20060519 18:24:24 nbsp wait where does the no toad sexing come in? Users/KrisFricke

20060729 17:35:10 nbsp No Toad Sexing: http://www.everything2.com/index.pl?node_id1078092 Users/PenguiN42

20061104 16:20:29 nbsp Having finally made it in to it last night, I can say that what struck me the most was how hot it is down there. Users/JosephBleckman

20080617 15:25:33 nbsp Toad sexing is also forbidden by commandment at MIT: http://hacks.mit.edu/by_year/1992/cathedral_7/father_tool_description.html Users/IDoNotExist

20100803 01:04:16 nbsp I discovered the Kerr/Wellman tunnel during my first year as an undergraduate, and later, as a math grad student, discovered to my delight that the key to the front door of Kerr also opened the door of said tunnel. Users/BarnabasTruman


text:
c $Id: comwid.f,v 1.15 2007/01/30 14:50:24 bleicher Exp $
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     Revision : 1.0
c
c     common block for the tabulated branching ratios
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c number of points for the interpolation
      integer widnsp
c the branching ratios are interpolated with high precision in the range
c from mintab to maxtab1 and from maxtab1 to maxtab2 with a smaller precision

c version number of the table
      integer tabver
      real*8 mintab, maxtab1, maxtab2

c set the default values
      parameter (widnsp=120)
      parameter (mintab=0.1d0)
      parameter (maxtab1=5.0d0)
      parameter (maxtab2=50.0d0)
c increase this parameter, if you make changes, which require a new table
      parameter (tabver=9)

c tabulated x-values (i.e. sqrt(s) of the collision)
      real*8 tabx (1:widnsp)
c tabulated y-values (i.e. the branching ratios) and the second
c derivatives of the function.
c full baryon ratio
      real*8 fbtaby (1:widnsp,minbar:maxbar,1:2)
c partial baryon ratios
      real*8 pbtaby (1:widnsp,1:2,minbar:maxbar,0:maxbra)
c full meson ratio
      real*8 fmtaby (1:widnsp,minmes:maxmes,1:2)
c partial meson ratios
      real*8 pmtaby (1:widnsp,1:2,minmes:maxmes,0:maxbrm)

c Breit-Wigner norms
c norm of Breit-Wigner with mass dependent widths baryons/mesons
        real*8 bwbarnorm(minbar:maxbar),bwmesnorm(minmes:maxmes)

c tabulated fppfit()
c tabulated x-values (i.e. sqrt(s) of the collision)
      real*8 tabxnd (1:widnsp)
c 2-resonance channels
c                      x     deriv ND N*..D*
      real*8 frrtaby(1:widnsp,1:2,1:2,2:maxdel)

c this flag indicates the progress of tabulating the function
      integer wtabflg
c name of file containing the tables
      character*77 tabname

      common /decaywidth/ tabx,fbtaby,pbtaby,fmtaby,pmtaby,wtabflg
        common /brwignorm/ bwbarnorm,bwmesnorm
        common /xsections/ tabxnd,frrtaby
      common /tabnames/ tabname

max_stars_repo_name: NRSummerfield/hic

max_issues_repo_name: NRSummerfield/hic

max_forks_repo_name: NRSummerfield/hic


text:
!< WenOOF test UI: definition of common User Interface (UI) for WenOOF tests.
module wenoof_test_ui
!< WenOOF test UI: definition of common User Interface (UI) for WenOOF tests.

use flap, only : command_line_interface
use penf, only : FR_P, I_P, R_P

implicit none
private
public :: test_ui

character(99), parameter :: interpolators(1:4) = ["JS  ", &
                                                  "M-JS", &
                                                  "M-Z ", &
                                                  "Z   "] !< List of available interpolators.

type :: test_ui
  !< Class to handle test(s) User Interface (UI).
  type(command_line_interface) :: cli                                  !< Command line interface handler.
  integer(I_P)                 :: error=0                              !< Error handler.
  character(99)                :: interpolator_type='JS'               !< Interpolator used.
  character(99)                :: output_bname='unset'                 !< Output files basename.
  character(99)                :: output_dir=''                        !< Output directory.
  integer(I_P)                 :: pn_number                            !< Number of different points-number tested.
  integer(I_P), allocatable    :: points_number(:)                     !< Points number used to discretize the domain.
  integer(I_P)                 :: S_number                             !< Number of different stencils tested.
  integer(I_P), allocatable    :: S(:)                                 !< Stencils used.
  real(R_P)                    :: eps                                  !< Small epsilon to avoid zero-division.
  real(R_P)                    :: x_target                             !< Interpolation target coordinate.
  logical                      :: interpolate=.false.                  !< Flag for activating interpolation.
  logical                      :: reconstruct=.false.                  !< Flag for activating reconstruction.
  logical                      :: errors_analysis=.false.              !< Flag for activating errors analysis.
  logical                      :: plots=.false.                        !< Flag for activating plots saving.
  logical                      :: results=.false.                      !< Flag for activating results saving.
  logical                      :: verbose=.false.                      !< Flag for activating verbose output.
  contains
    ! public methods
    procedure, pass(self) :: get               !< Get user options.
    procedure, pass(self) :: loop_interpolator !< Loop over available interpolators.
endtype test_ui

contains
  ! public methods
  subroutine get(self)
  !< Get user options.
  class(test_ui), intent(inout) :: self !< Test UI.

  call set_cli
  call parse_cli
  contains
    subroutine set_cli()
    !< Set Command Line Interface.

    associate(cli => self%cli)
      call cli%init(progname    = 'WenOOF Test',                                    &
                    authors     = 'Fortran-FOSS-Programmers',                       &
                    license     = 'GNU GPLv3',                                      &
                    description = 'Test WenOOF library on function reconstruction', &
                    examples    = ["$EXECUTABLE --interpolate -i JS --results",     &
                                   "$EXECUTABLE --interpolate -i JS-Z -r     ",     &
                                   "$EXECUTABLE --reconstruct -i JS-M        ",     &
                                   "$EXECUTABLE --reconstruct -i all -p -r   "])

      call cli%add(switch='--interpolate', help='Perform interpolation', required=.false., act='store_true', def='.false.')
      call cli%add(switch='--reconstruct', help='Perform reconstruction', required=.false., act='store_true', def='.false.')
      call cli%add(switch='--x_target', switch_ab='-x', help='WENO interpolation target point coordinate', &
                   required=.false., def='0.0', act='store')
      call cli%add(switch='--interpolator', switch_ab='-i', help='WENO interpolator type', required=.false., &
                   def='JS', act='store', choices='all,JS,M-JS,M-Z,Z')
      call cli%add(switch='--points_number', switch_ab='-pn', nargs='+', help='Number of points used to discretize the domain', &
                   required=.false., act='store', def='50 100')
      call cli%add(switch='--stencils', switch_ab='-s', nargs='+', help='Stencils dimensions (and number)', required=.false., &
                   act='store', def='2 3 4 5 6 7 8 9', choices='2, 3, 4, 5, 6, 7, 8, 9')
      call cli%add(switch='--eps', help='Small epsilon to avoid zero-division', required=.false., act='store', def='1.e-6')
      call cli%add(switch='--output_dir', help='Output directory', required=.false., act='store', def='./')
      call cli%add(switch='--results', switch_ab='-r', help='Save results', required=.false., act='store_true', def='.false.')
      call cli%add(switch='--plots', switch_ab='-p', help='Save plots', required=.false., act='store_true', def='.false.')
      call cli%add(switch='--output', help='Output files basename', required=.false., act='store', def='output')
      call cli%add(switch='--errors_analysis', help='Peform errors analysis', required=.false., act='store_true', def='.false.')
      call cli%add(switch='--verbose', help='Verbose output', required=.false., act='store_true', def='.false.')
    endassociate
    endsubroutine set_cli

    subroutine parse_cli()
    !< Parse Command Line Interface and check its validity.

    call self%cli%parse(error=self%error) ; if (self%error/=0) stop
    call self%cli%get(switch='--interpolate', val=self%interpolate, error=self%error) ; if (self%error/=0) stop
    call self%cli%get(switch='--reconstruct', val=self%reconstruct, error=self%error) ; if (self%error/=0) stop
    call self%cli%get(switch='-x', val=self%x_target, error=self%error) ; if (self%error/=0) stop
    call self%cli%get(switch='-i', val=self%interpolator_type, error=self%error) ; if (self%error/=0) stop
    call self%cli%get_varying(switch='-pn', val=self%points_number, error=self%error) ; if (self%error/=0) stop
    call self%cli%get_varying(switch='-s', val=self%S, error=self%error) ; if (self%error/=0) stop
    call self%cli%get(switch='--eps', val=self%eps, error=self%error) ; if (self%error/=0) stop
    call self%cli%get(switch='--output_dir', val=self%output_dir, error=self%error) ; if (self%error/=0) stop
    call self%cli%get(switch='-r', val=self%results, error=self%error) ; if (self%error/=0) stop
    call self%cli%get(switch='-p', val=self%plots, error=self%error) ; if (self%error/=0) stop
    call self%cli%get(switch='--output', val=self%output_bname, error=self%error) ; if (self%error/=0) stop
    call self%cli%get(switch='--errors_analysis', val=self%errors_analysis, error=self%error) ; if (self%error/=0) stop
    call self%cli%get(switch='--verbose', val=self%verbose, error=self%error) ; if (self%error/=0) stop

    if ((.not.self%interpolate).and.(.not.self%reconstruct)) then
       ! activate both for coverage tests
       self%interpolate = .true.
       self%reconstruct = .true.
    endif

    self%pn_number = size(self%points_number, dim=1)
    self%S_number = size(self%S, dim=1)
    endsubroutine parse_cli
  endsubroutine get

  function loop_interpolator(self, interpolator) result(again)
  !< Loop over available interpolators.
  class(test_ui), intent(in)  :: self          !< Test UI.
  character(99),  intent(out) :: interpolator  !< Interpolator name.
  logical                     :: again         !< Flag continuing the loop.
  integer(I_P), save          :: i = 0         !< Counter.

  again = .false.
  if (i==0) then
    i = 1
    interpolator = trim(adjustl(interpolators(i)))
    again = .true.
  elseif (i<size(interpolators, dim=1)) then
    i = i + 1
    interpolator = trim(adjustl(interpolators(i)))
    again = .true.
  else
    i = 0
    interpolator = ''
    again = .false.
  endif
  endfunction loop_interpolator
endmodule wenoof_test_ui

max_stars_repo_name: Fortran-FOSS-Programmers/WenOOF

max_issues_repo_name: Fortran-FOSS-Programmers/WenOOF

max_forks_repo_name: Fortran-FOSS-Programmers/WenOOF


Module halo_serial_module

  Use, Intrinsic :: iso_fortran_env, Only :  wp => real64

  Use halo_setter_base_module, Only : halo_setter_base_class, halo_setter_base_data_class

  Implicit None

  Type, Public, Extends( halo_setter_base_data_class ) :: halo_serial_data
  End type halo_serial_data

  ! Mostly dummy type for the moment, eventually will contain things like MPI communicator
  Type, Public, Extends( halo_setter_base_class ) :: halo_serial_setter
     Private
   Contains
     Procedure, Public :: fill => halo_fill
  End Type halo_serial_setter

  Private

Contains

  Subroutine halo_fill( H, halo_width, hdlb, gin, hout )

    Class( halo_serial_setter ),                                 Intent( In    ) :: H
    Integer   ,                                                  Intent( In    ) :: halo_width
    Integer   , Dimension( 1:3 ),                                Intent( In    ) :: hdlb
    Real( wp ), Dimension( 0:, 0:, 0: )                        , Intent( In    ) :: gin
    Real( wp ), Dimension( hdlb( 1 ):, hdlb( 2 ):, hdlb( 3 ): ), Intent(   Out ) :: hout

    Integer, Dimension( 1:3 ) :: gub, gn
    Integer, Dimension( 1:3 ) :: hub, hlb

    Integer, Dimension( 1:3 ) :: hi, gi

    Integer :: hi1, hi2, hi3

    Logical, Parameter :: debug = .True.

    ! Bounds on grid without halo - we know the lower bound is zero from declaration
    gub = Ubound( gin )

    ! Bounds on grid with halo
    ! Have to be careful as HOUT might be declared bigger than needed
    hlb =     - halo_width
    hub = gub + halo_width

    ! Size
    gn = gub + 1

    ! Fill the central section
    hout( 0:gub( 1 ), 0:gub( 2 ), 0:gub( 3 ) ) = gin

    ! Need to finish actually filling halos!!!
    ! First all the stuff below the bottom face of the non-haloed grid
    Do hi3 = hlb( 3 ), -1
       Do hi2 = hlb( 2 ), hub( 2 )
          Do hi1 = hlb( 1 ), hub( 1 )
             hi = [ hi1, hi2, hi3 ]
             gi = Modulo( hi, gn )
             hout( hi( 1 ), hi( 2 ), hi( 3 ) ) = gin( gi( 1 ), gi( 2 ), gi( 3 ) )
          End Do
       End Do
    End Do

    ! Now above the top face
    Do hi3 = gub( 3 ) + 1, hub( 3 )
       Do hi2 = hlb( 2 ), hub( 2 )
          Do hi1 = hlb( 1 ), hub( 1 )
             hi = [ hi1, hi2, hi3 ]
             gi = Modulo( hi, gn )
             hout( hi( 1 ), hi( 2 ), hi( 3 ) ) = gin( gi( 1 ), gi( 2 ), gi( 3 ) )
          End Do
       End Do
    End Do

    ! Now the front face
    Do hi3 = 0, gub( 3 )
       Do hi2 = hlb( 2 ), -1
          Do hi1 = hlb( 1 ), hub( 1 )
             hi = [ hi1, hi2, hi3 ]
             gi = Modulo( hi, gn )
             hout( hi( 1 ), hi( 2 ), hi( 3 ) ) = gin( gi( 1 ), gi( 2 ), gi( 3 ) )
          End Do
       End Do
    End Do

    ! Back face
    Do hi3 = 0, gub( 3 )
       Do hi2 = gub( 2 ) + 1, hub( 2 )
          Do hi1 = hlb( 1 ), hub( 1 )
             hi = [ hi1, hi2, hi3 ]
             gi = Modulo( hi, gn )
             hout( hi( 1 ), hi( 2 ), hi( 3 ) ) = gin( gi( 1 ), gi( 2 ), gi( 3 ) )
          End Do
       End Do
    End Do

    ! Left Face
    Do hi3 = 0, gub( 3 )
       Do hi2 = 0, gub( 2 )
          Do hi1 = hlb( 1 ), -1
             hi = [ hi1, hi2, hi3 ]
             gi = Modulo( hi, gn )
             hout( hi( 1 ), hi( 2 ), hi( 3 ) ) = gin( gi( 1 ), gi( 2 ), gi( 3 ) )
          End Do
       End Do
    End Do

    ! Right Face
    Do hi3 = 0, gub( 3 )
       Do hi2 = 0, gub( 2 )
          Do hi1 = gub( 1 ) + 1, hub( 1 )
             hi = [ hi1, hi2, hi3 ]
             gi = Modulo( hi, gn )
             hout( hi( 1 ), hi( 2 ), hi( 3 ) ) = gin( gi( 1 ), gi( 2 ), gi( 3 ) )
          End Do
       End Do
    End Do

    ! Sanity check on fiddly routine
    If( debug ) Then
       Do hi3 = hlb( 3 ), hub( 3 )
          Do hi2 = hlb( 2 ), hub( 2 )
             Do hi1 = hlb( 1 ), hub( 1 )
                hi = [ hi1, hi2, hi3 ]
                gi = Modulo( hi, gn )
                If( Abs( hout( hi( 1 ), hi( 2 ), hi( 3 ) ) - &
                     gin( gi( 1 ), gi( 2 ), gi( 3 ) ) ) > Tiny( hout ) ) Then
                   Write( *, '( a, 1x, 2( 3( i5, 1x ), g20.12, 5x ) )'  ) &
                        'halo filling problem',                           &
                        hi, hout( hi( 1 ), hi( 2 ), hi( 3 ) ),            &
                        gi,  gin( gi( 1 ), gi( 2 ), gi( 3 ) )
                End If
             End Do
          End Do
       End Do
    End If

  End Subroutine halo_fill

End Module halo_serial_module


!> Simple MPI program to test the timing API when not all ranks
!! enter a given timed region
PROGRAM timer_test
  use dl_timer_constants_mod
  use dl_timer
  use mpi
  implicit none

  integer :: time0
  integer(i_def64), parameter :: nstep = 5

  real(r_def) :: mysum

  integer numtasks, rank, ierr, j

  !--------------------------------------------------------------
  ! Initialisation

  ! Initialize the MPI library:
  call MPI_INIT(ierr)
  if (ierr .ne. MPI_SUCCESS) then
     print *,'Error starting MPI program. Terminating.'
     call MPI_ABORT(MPI_COMM_WORLD, 1, ierr)
  end if

  ! Get the number of processors this job is using
  call MPI_COMM_SIZE(MPI_COMM_WORLD, numtasks, ierr)

  ! Get the rank of the processor this thread is running on
  call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)

  ! Initialise the timer library (must be done *after* MPI has been
  ! initialised)
  call timer_init()

  mysum = 0.0d0

  !--------------------------------------------------------------
  ! Time-stepping

  call timer_start(time0, label='Time-stepping', num_repeats=nstep)

  ! Don't do time-stepping on odd ranks so that not all PEs enter the
  ! 'fake section' region
  if (mod(rank,2) == 0) call step()

  call timer_stop(time0)

  !---------------------------------------------------------------
  ! Finalise

  ! Output the result to prevent the compiler optimising-out the loop
  write (*,*) 'Fake checksum = ', mysum

  call timer_report()

  ! Tell the MPI library to release all resources it is using:
  call MPI_FINALIZE(ierr)

contains

  !--------------------------------------------------------------
  subroutine step()
    use timer_test_utils_mod, only: fsleep
    implicit none
    ! Locals
    integer :: istep, time1, nloops, ret

    ! Artificially create a load imbalance
    nloops = 200*(rank+1)

    do istep = 1, nstep

       call timer_start(time1, label='Fake section')
       do j = 1, nloops
          mysum = mysum + sqrt(5.0d0*istep*istep)
       end do
       ret = fsleep(1)
       call timer_stop(time1)

    end do

  end subroutine step

END PROGRAM timer_test


text:
great_ne(dd1,v1).
great_ne(aa1,y1).
great_ne(aa1,i1).
great_ne(kk1,hh1).
great_ne(k1,x1).
great_ne(kk1,bb1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


        !COMPILER-GENERATED INTERFACE MODULE: Fri Mar  6 15:34:10 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE DIAG2__genmod
          INTERFACE
            SUBROUTINE DIAG2(SUM,AVN,AVS)
              REAL(KIND=8) :: SUM(112)
              REAL(KIND=8) :: AVN
              REAL(KIND=8) :: AVS
            END SUBROUTINE DIAG2
          END INTERFACE
        END MODULE DIAG2__genmod


text: #redirect AT&T Wireless


text:
MODULE m_forceb8
  !-------------------
  ! Implements the surface contribution to the force. Madsen Eq.(B8)
  !
  ! FZJ 15/3-01 GMadsen
  !---------------------------------
CONTAINS
  SUBROUTINE force_b8(atoms,ecwk,stars, sym,cell, jspin, force,f_b8)

    USE m_constants, ONLY : tpi_const
    USE m_sphbes
    USE m_stern
    USE m_types
    IMPLICIT NONE
    TYPE(t_sym),INTENT(IN)       :: sym
    TYPE(t_stars),INTENT(IN)     :: stars
    TYPE(t_cell),INTENT(IN)      :: cell
    TYPE(t_atoms),INTENT(IN)     :: atoms
    !     ..
    !     .. Arguments
    INTEGER, INTENT (IN) :: jspin
    COMPLEX, INTENT (IN) :: ecwk(stars%ng3)
    COMPLEX, INTENT (INOUT) :: f_b8(3,atoms%ntype)
    real,intent(inout)      :: force(:,:,:)
    !     ..
    !     .. Local Variables
    INTEGER g(3),nst,stg(3,sym%nop),ia,istr,i,j,jj,jneq
    REAL    fj(0:atoms%lmaxd),rotkzz(3),rstg(3,sym%nop)
    REAL    frmt,gl,pha,s
    COMPLEX taup(sym%nop),factor,fact,fstar(3),fsur2(3)
    !

    ia=1
    DO jneq=1,atoms%ntype
       frmt = 2.0*tpi_const*atoms%rmt(jneq)**2
       fsur2(1:3) = cmplx(0.0,0.0)

       ! skip G=(0,0,0) no contribution to ekin
       DO istr=2,stars%ng3_fft

          g(:)     = stars%kv3(:,istr)
          fstar(:) = cmplx(0.0,0.0)
          CALL stern(sym,cell,g, nst,stg,taup,gl,rstg)

          CALL sphbes(atoms%lmax(jneq),atoms%rmt(jneq)*gl,fj)
          fact = ecwk(istr) * fj(1) / gl

          DO jj=1,nst
             pha=(atoms%taual(1,ia)*stg(1,jj)+atoms%taual(2,ia)*stg(2,jj)&
                  +atoms%taual(3,ia)*stg(3,jj))*tpi_const
             !
             !cc   swaped sin and cos because there's an i in the equation
             !
             factor = fact * cmplx(-sin(pha),cos(pha)) * taup(jj)
             DO i=1,3
                fstar(i) = fstar(i) + factor*rstg(i,jj)
             ENDDO
          ENDDO
          DO i=1,3
             fsur2(i)=fsur2(i)+fstar(i)*frmt
          ENDDO
       ENDDO
       DO i=1,3
          f_b8(i,jneq) = f_b8(i,jneq) + fsur2(i)
          force(i,jneq,jspin) = force(i,jneq,jspin) + real(fsur2(i))
       ENDDO
       ia=ia+atoms%neq(jneq)
    ENDDO

  END SUBROUTINE force_b8
END MODULE m_forceb8

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR