text:
!
!  getinf_el.f
!
!
!  Created by Xingshi on 1/26/08.
!  Copyright 2008 __MyCompanyName__. All rights reserved.
!

       subroutine getinf_el_3d(finf, x, xna, nn, nsd, ne, nen, ien, maxconn)
	   implicit none
	   integer :: ninf, nn, nsd, ne, nen, maxconn ! fluid variable from main
	   integer ien(nen,ne) ! transfer matrix connectted local information to global
	   integer inf(maxconn) ! subroutine variable to store the index of element which contains the solid node
	   real(8) x(nsd), xna(nsd, nn) ! x - coordinates of solid node   xna - coordinates of fluid nodes in global view
           real(8) xn(nsd,nen)
           integer nnum
	   real(8) length_max(nsd)  ! the minimun rectangular to contian the element
	   real(8) length_min(nsd)
	   real(8) xna_p1(nsd)      ! coordinates of two neighbour points in the element
	   real(8) xna_p2(nsd)
	   real(8) xna_p3(nsd)
	   real(8) vector_0(nsd)
	   real(8) vector_1(nsd)
	   real(8) vector_2(nsd)
	   real(8) dot00
	   real(8) dot01
	   real(8) dot02
	   real(8) dot11
	   real(8) dot12
	   real(8) invDenom
	   real(8) sign_1
	   real(8) sign_2
	   real(8) vector_fe1(nsd)
	   real(8) vector_fe2(nsd)
	   real(8) vector_p(nsd)
	   real(8) sign_fe
	   integer finf  ! the final result we want to get
	   integer ie  ! loop variables
	   integer inl
	   integer isd
	   integer i

       ninf = 0
     !  write(*,*) 'nen', nen
     !  write(*,*) '1000 el in get ', xna(:,124)
	   do ie = 1, ne
	      do inl = 1, nen

                     nnum=ien(inl,ie)
              ! write(*,*) 'nnum=  ', nnum
              ! write(*,*) 'ie=    ',ie
              ! write(*,*) 'inl=   ',inl
		         do isd = 1, nsd
		!	write(*,*) '485 x', xna(1,485)
			 xn(isd, inl) = xna(isd, nnum)
                   !      write(*,*) 'ne=',ie, 'nnum=', xna(2,nnum)
			 end do
	      end do

		  do isd = 1, nsd
		  length_max(isd) = maxval(xn(isd,1:nen))
		  length_min(isd) = minval(xn(isd,1:nen))
		!  write(*,*) 'xna = ', xn(isd,1), xn(isd,2), xn(isd,3), xn(isd,4)
                  end do

		  if (nsd==2) then
		     if (x(1)>length_min(1) .and. x(1)<length_max(1) .and. &
                         x(2)>length_min(2) .and. x(2)<length_max(2)) then
				 ninf = ninf +1
		         inf(ninf) = ie
                         !write(*,*) 'finf',ie
                        ! write(*,*) 'x-y  ', length_min(2), length_max(2)
			  end if
		   end if

		  if (nsd==3) then
		     if (x(1)>length_min(1) .and. x(1)<length_max(1) &
                         .and. x(2)>length_min(2) .and. x(2)<length_max(2) .and. &
				 x(3)>length_min(3) .and. x(3)<length_max(3)) then
				 ninf = ninf +1
		         inf(ninf) = ie

			  end if
		   end if
	  !write(*,*) 'x min' ,length_min(1)
         ! write(*,*) 'x max' , length_max(1)
         ! write(*,*) 'x 4 nodal4', xn(1,1), xn(1,2), xn(1,3), xn(1,4)
          ! write(*,*) 'ie= ', ie
           end do
         !  write(*,*) 'ninf= ', ninf

         !  write(*,*) 'inf= ' , inf(1)

	  if(ninf == 0) goto 2000

	   if (nsd==2 .and. nen==3) then
	      do i= 1, ninf
				   xna_p1(1:nsd) = xna(1:nsd, ien(1,inf(i)))
				   xna_p2(1:nsd) = xna(1:nsd, ien(2,inf(i)))
				   xna_p3(1:nsd) = xna(1:nsd, ien(3,inf(i)))

			 vector_0(1:nsd) = xna_p3(1:nsd)-xna_p1(1:nsd)
			 vector_1(1:nsd) = xna_p2(1:nsd)-xna_p1(1:nsd)
			 vector_2(1:nsd) = x(1:nsd)-xna_p1(1:nsd)

			 dot00=vector_0(1)*vector_0(1)+vector_0(2)*vector_0(2)
			 dot01=vector_0(1)*vector_1(1)+vector_0(2)*vector_1(2)
			 dot02=vector_0(1)*vector_2(1)+vector_0(2)*vector_2(2)
			 dot11=vector_1(1)*vector_1(1)+vector_1(2)*vector_1(2)
			 dot12=vector_1(1)*vector_2(1)+vector_1(2)*vector_2(2)
			 invDenom=1/(dot00*dot11-dot01*dot01)
			 sign_1=(dot11*dot02-dot01*dot12)*invDenom
			 sign_2=(dot00*dot12-dot01*dot02)*invDenom

			 if(sign_1>0.0 .and. sign_2>0.0 .and. (sign_1+sign_2)<1.0) then
			 finf = inf(i)
			 goto 2000
			 end if
		  end do
           end if
          if (nsd==2 .and. nen==4) then
                 ! write(*,*) 'I am in the right loop'

                 ! write(*,*) 'get  x=' ,x(:)
                  do i= 1, ninf
                 ! write(*,*) 'xna= ',i, xna(1:nsd,ien(1:nen,inf(i)))
                  vector_fe1(1:nsd)=xna(1:nsd,ien(2,inf(i)))-xna(1:nsd,ien(1,inf(i)))
                  vector_fe2(1:nsd)=xna(1:nsd,ien(3,inf(i)))-xna(1:nsd,ien(1,inf(i)))
                  vector_p(1:nsd)=x(1:nsd)-xna(1:nsd,ien(1,inf(i)))
                  sign_1=vector_fe1(1)*vector_fe2(2)-vector_fe1(2)*vector_fe2(1)
                  sign_2=vector_fe1(1)*vector_p(2)-vector_fe1(2)*vector_p(1)
                  sign_fe=sign_2*sign_1
                 ! write(*,*) '1-2', sign_1, sign_2
                  if (sign_fe<0)  then
                  go to 2100
                  end if
                  !1-2
                  vector_fe1(1:nsd)=xna(1:nsd,ien(3,inf(i)))-xna(1:nsd,ien(2,inf(i)))
                  vector_fe2(1:nsd)=xna(1:nsd,ien(4,inf(i)))-xna(1:nsd,ien(2,inf(i)))
                  vector_p(1:nsd)=x(1:nsd)-xna(1:nsd,ien(2,inf(i)))
                  sign_1=vector_fe1(1)*vector_fe2(2)-vector_fe1(2)*vector_fe2(1)
                  sign_2=vector_fe1(1)*vector_p(2)-vector_fe1(2)*vector_p(1)
                  sign_fe=sign_2*sign_1
                 ! write(*,*) '2-3', sign_1, sign_2
                  if (sign_fe<0)  then
                  go to 2100
                  end if
                  !2-3

                  vector_fe1(1:nsd)=xna(1:nsd,ien(4,inf(i)))-xna(1:nsd,ien(3,inf(i)))
                  vector_fe2(1:nsd)=xna(1:nsd,ien(1,inf(i)))-xna(1:nsd,ien(3,inf(i)))
                  vector_p(1:nsd)=x(1:nsd)-xna(1:nsd,ien(3,inf(i)))
                  sign_1=vector_fe1(1)*vector_fe2(2)-vector_fe1(2)*vector_fe2(1)
                  sign_2=vector_fe1(1)*vector_p(2)-vector_fe1(2)*vector_p(1)
                  sign_fe=sign_2*sign_1
                 ! write(*,*)'3-4',sign_1,sign_2
                  if (sign_fe<0)  then
                  go to 2100
                  end if
                  !3-4

                  vector_fe1(1:nsd)=xna(1:nsd,ien(1,inf(i)))-xna(1:nsd,ien(4,inf(i)))
                  vector_fe2(1:nsd)=xna(1:nsd,ien(2,inf(i)))-xna(1:nsd,ien(4,inf(i)))
                  vector_p(1:nsd)=x(1:nsd)-xna(1:nsd,ien(4,inf(i)))
                  sign_1=vector_fe1(1)*vector_fe2(2)-vector_fe1(2)*vector_fe2(1)
                  sign_2=vector_fe1(1)*vector_p(2)-vector_fe1(2)*vector_p(1)
                  sign_fe=sign_2*sign_1
                 ! write(*,*) '4-1',sign_1,sign_2
                  if (sign_fe<0)  then
                  go to 2100
                  end if
                  !4-1
                  finf=inf(i)
                  go to 2000
2100              continue
                  end do
	  end if

	  if (nsd==3) then
!          write(*,*) 'i am here'
          call search_3d(finf, x, xna, nn, nsd, ne, nen, ien, inf, ninf,maxconn) ! call search_3d for 3d
          !cases
         end if
		2000		continue
	!	write (*,*) 'Find the elemet in fluid domain as:',  finf
       return
       end

max_stars_repo_name: biofluids/IFEM-archive

max_issues_repo_name: biofluids/IFEM-archive

max_forks_repo_name: biofluids/IFEM-archive


text:
      FUNCTION CHABIMF(ALO,AUP,M1,M2,STARMASS)

c     Returns number of stars in range from M1 to M2 for Chabrier (2001) IMF
c     Returns the mass in stars from M1 to M2
c     ALO and AUP are the lower and upper cutoffs of the IMF, which
c     is normalized to 1 Mo over this range

      implicit real (A-H,O-Z)
      real*8 m1,m2,starmass

      common /fmlog/xml,sigmac,anormlog
      common /fmlin/anormlin

*******************************************************
*                 - Adjustable parameters -
*
*    Lower and upper masses for the entire mass
*    integration of the IMF: alo, aup
*    Index "alpha=1+x" in the power-law part > 1 Msol
*           dn/dm propto m**(-alpha)
      alpha=2.3
*******************************************************

      xml=-1.1
      sigmac=0.69
      call chabint(alo,aup,alpha)    ! calculates the normalization
      al1=1.-alpha
      al2=2.-alpha
      anorm=anormlin/log(10.)
      v1log=anormlog*exp(-xml**2/2./sigmac**2)
c      print*,'normlog, normlin',anormlog,anormlin
c      print*,'norm @ 1 Msol',v1log,anormlin
      if(v1log.ne.anormlin) then
                            print*,'normalization pb in routine chabimf!'
                            stop
      endif

	if (m2.le.m1) then
		chabimf=0.
		starmass=0.
		return
	endif

      IF(m1.LE.1.0) THEN
            if(m2.LE.1.0) then
                vnum=chabnum(m1,m2)     ! number-density
                vmas=chabmas(m1,m2)     ! mass-density
            else
                alim=1.0
                vnum=chabnum(m1,dble(alim))+anorm/al1*(m2**al1-1.0)
                vmas=chabmas(m1,dble(alim))+anorm/al2*(m2**al2-1.0)
            endif
      ELSE
                vnum=anorm/al1*(m2**al1-m1**al1)
                vmas=anorm/al2*(m2**al2-m1**al2)
      ENDIF

      chabimf=vnum
      starmass=vmas

      return
      end

      SUBROUTINE CHABINT(alo,aup,alpha)
**************************************************************************
**    chabmas calculates the normalization such that
**        [int m.(dn/dm)dm = 1]  from m=alo to m=aup
**************************************************************************
      implicit real (A-H,O-Z)
      common /fmlog/xml,sigmac,anormlog
      common /fmlin/anormlin

      pi=4.*atan(1.)
      cte2=exp(-xml**2/sigmac**2/2.)
      al2=2.-alpha
      clin=cte2/al2*(aup**al2-1.)
      if(aup.LE.1.) clin=0.

      sigmacp=sigmac*log(10.)
      xmlp=xml*log(10.)
      cte1=exp(xmlp+sigmacp**2/2.)
      x1=-(xmlp+sigmacp**2)/sqrt(2.)/sigmacp
      xlo=(log10(alo)*log(10.)-(xmlp+sigmacp**2))/sqrt(2.)/sigmacp
      clog=cte1*sqrt(2.*pi)/2.*sigmacp
     !       *  (erf(x1) - erf(xlo))

      anorm=1./(clog+clin)
      anormlog=anorm*log(10.)     ! cte de la partie lognormal
      anormlin=anormlog*cte2         ! cte de la partie loi de puissance

      return
      end

      FUNCTION CHABNUM(xmin,xmax)
**************************************************************************
**                         *    IMF Chabrier   *
**    chabnum returns [int (dn/dm)dm]  from m=xmin to m=xmax
*************************************************************************
      implicit real (A-H,O-Z)
	real*8 xmin,xmax
c	real*8 y1,y2,de
	include 'chabimf.dec'
      common /fmlog/xml,sigmac,anormlog

c	check mass range
	if (xmax.le.xmin) then
		chabnum=0.
		return
	else
		pi=4.*atan(1.)
c		Modified by GBA, April 19, 2003
c		chabnum=anormlog * sqrt(2.*pi)/2.*sigmac
c    !		*( erf((log10(xmax)-xml)/sqrt(2.)/sigmac )
c    !		-  erf((log10(xmin)-xml)/sqrt(2.)/sigmac ) )
		y1 =  (dlog10(xmax*1.D00)-xml)/sqrt(2.)/sigmac
		y2 =  (dlog10(xmin*1.D00)-xml)/sqrt(2.)/sigmac
		de =  erf(y1) - erf(y2)
		chabnum=anormlog * sqrt(2.*pi)/2.*sigmac * de
		return
	endif
	end

      FUNCTION CHABMAS(xmin,xmax)
**************************************************************************
**                         *    IMF Chabrier   *
**    chabmas returns [int m.(dn/dm)dm]  from m=xmin to m=xmax
*************************************************************************
      implicit real (A-H,O-Z)
	real*8 xmin,xmax
c	real*8 y1,y2,de
	include 'chabimf.dec'
      common /fmlog/xml,sigmac,anormlog

c	Check mass range
	if (xmax.le.xmin) then
		chabmas=0.
		return
	else
		pi=4.*atan(1.)
		sigmacp=log(10.)*sigmac
		xmlp=xml*log(10.)
		anormp=anormlog/log(10.)
		cte=exp(xmlp+sigmacp**2/2.)
c		Modified by GBA, April 19, 2003
c		x1=(log10(xmin)*log(10.)-(xmlp+sigmacp**2))/sqrt(2.)/sigmacp
c		x2=(log10(xmax)*log(10.)-(xmlp+sigmacp**2))/sqrt(2.)/sigmacp
c		chabmas=anormp*cte*sqrt(2.*pi)/2.*sigmacp * (erf(x2) - erf(x1))
		y1=(dlog10(xmin*1.D00)*dlog(10.D00)-(xmlp+sigmacp**2))/sqrt(2.)/sigmacp
		y2=(dlog10(xmax*1.D00)*dlog(10.D00)-(xmlp+sigmacp**2))/sqrt(2.)/sigmacp
		de=erf(y2) - erf(y1)
		chabmas=anormp * cte * sqrt(2.*pi)/2.*sigmacp * de
		return
	endif
	end

ext: f

max_stars_repo_name: jchavesmontero/galaxpy

max_issues_repo_name: jchavesmontero/galaxpy

max_forks_repo_name: jchavesmontero/galaxpy


text:
If you would like to ask a computer question, or answer one, you may do so here. Also any information relevant to computers that can be considered helpful can go here.

 Talk Topics to date
Click to jump to the section
#A Laptop Coolers  are they necessary?
#B Proxies to circumvent blocked websites
#C Removal of Malware, security software, reformatting
#D Spilled liquid onto laptop, slightly malfunctioning
#E Refurbished computer equipment: good buy, or too risky?
#F Sluggish laptop  worth the memory upgrade?
#G Hardware questions  power to Usbs?

  If you have other questions questions, ask the wiki and post below!
  For serious issues, and specific questions, it may be worth your while to contact Computer repair a professional.

Anchor(A)20090519 15:34:10 nbsp Laptop coolers: gimmick or gottahave? Users/TheAmazingLarry

  20090519 15:39:08 nbsp Depends on the machine and how you plan on using it. Using a laptop on your bed for example often causes it to suck in dust, lint, and possibly hair potentially causing it to overheat.  Some performance intensive applications such as gaming can cause high levels of heat over prolonged periods of time which also has the potential to cause overheating.  Most often, and for most applications a cooler is completely useless.  If you MUST work on your bed, a lap desk, or large hard cover book can raise the computer up and allow for better ventilation.  On the other hand, some models (Such as older dells with air intake on the bottom) are notorious for overheating even under the best conditions and would benefit from a cooling stand.   Users/MasonMurray

  20090519 15:48:24 nbsp A metal cookie sheet can also work well. Users/NickSchmalenberger
   ...or pretty much anything that allows air to circulate under the laptop, right?  I mean, the fans in the laptop coolers are pretty much there just to make them look hightech, and not because the fans add much cooling power, right?

  I had a Mac whose fan came on all of the time, sometimes just while browsing the internet or watching a video.  Using a laptop pad (without its own fan) to elevate the laptop definitely reduced the amount of time the fan needed to come on.  Best case: lengthen life of harddrive.  Worst case: less time listening to annoying fan noise.  I use a Targus coolpad.   Users/CovertProfessor

  Exactly as Mason said, it depends on what you have and how youre using it.  I wouldnt use anything like a cookie sheet or anything in direct contact under the laptop (unless you get elaborate with it), but even simply propping a book under the back of it to give it an inch of air likely helps.  Lowend coolers are gimmicky, the fans tend to suck.  Theyre helpful but not much more than other stuff you can do.  There actually are nice coolers that work great out there (and really are useful), but it again depends on what youre doing to make the laptop run that hot.  users/edwins ES

  Slightly vulgar/NSFW, but you really gotta be careful Larry. This has was a huge issue when the new Macbook Pros http://www.pennyarcade.com/comic/2009/1/30/ came out.

  20090520 14:21:27 nbsp If anyone is having issues with cooling or fans being too noisey in a laptop, it is very likely that your laptop needs to be cleaned. The heatsink for the CPU in a laptop is quite small. Dust can easily and quickly build up around the heatsink and eventually restrict all airflow leaving the CPU to stagnate in its own hot air. If you are having heating issues, I would not reccomend putting it on a cookie sheet or anything like that as a serious solution. If your laptop is experiencing this type of issue, it is only going to get worse and as time goes by, your CPU is spending a lot of time either over heating or lowering its clock speed to avoid over heating. If you have a CPU with more than one core, which is almost for sure if your laptop is less that 3 years old, then you should use a program called CoreTemp to see if your CPU is over heating or not. Core temp is available here: http://www.alcpu.com/CoreTemp/ Once you have CoreTemp downlaoded and running, you need to stress your computer to find out if it over heats when it is being used or not. To do this, you simply need to sun a stress test while you run CoreTemp. A good stress test can be found here: http://www.pcworld.com/downloads/file/fid,7146order,1page,1c,systemresourcestuneup/description.html

 If you do not have any knowledge about this sort of thing, do not try to clean the fan yourself. You can try Dust off or sucking air through the fan hole with the a vaccum hose but dont try to take anything apart!

 This is a very common issue with pretty much all laptops and its not especially difficult or expensive to fix. Users/RealComputers
     I have an old can of compressed air at home, but it has some additive to keep youngsters from huffing it.  Can I still spray that into my laptop to get the dust out? Or do I need to buy the powderfree canned air?
        You can use any kind of compressed air. Yes some laptop fans are louder than others, but when a CPU is overheating, the CPU fan will spin faster and make more noise as it tries to compensate. This would be quite noticable and much louder than the fan would normally be regardless of how noisey it origionally was.  Users/RealComputers
  Brian, it seems to me that your answer overlooks the wide variety of laptops that are out there.  Some laptop fans are noisier than others right out of the box.  Some fans come on frequently, some almost never.  I had a Powerbook whose fan came on all the time, and all the Powerbooks of that generation were like that; my current MBP fan hardly ever comes on.  And I just got out of a meeting with someone who had an HP laptop that sounded like a jet engine.  In the worst of days, my Powerbook never sounded like that.   In short, noisy fan doesnt necessarily mean that there is dust, and fan coming on could just be a product of a particular model together with certain types of use.  Users/CovertProfessor

   What you just stated is why people should follow the instructions above if they have a noisy fan and suspect their computer is overheating. Doing so will allow them to determine if there is a problem or not. Users/RealComputers
    That site appears to be Windows only.   Users/CovertProfessor

   20090521 00:01:49 nbsp Wiping dust out with a tissue or paper towel can often work better than canned air, an important part is getting access to where the dust is and giving it room to leave. I wouldnt just blow air into a laptop without opening it because it might blow the dust in farther. If I have access to an air compressor I definitely use that too. Users/NickSchmalenberger
      The fan inside a laptop is sealed in a small metal chamber. Dust can not make it into other parts of the laptop from the fan vents.  Users/RealComputers
 Again you are making a generalization about the hardware arrangement of laptops, which in my experience is very diverse. If the fan is sealed with only one opening then blowing air in could be the only way to get dust out. Many laptops are not like this however. For example, in http://www.flickr.com/photos/macpolen/2256203086/ this picture I found in a google image search for ibook fan, it is clearly open on top so even when the keyboard is on dust will probably go through that opening. Whether or not the fan blows in or out also makes a difference and this too will vary between laptops.  Users/NickSchmalenberger
   That picture is of a laptop that is not assembled .. And yes I am generalizing, however all laptops fans are sealed and dust can not get inside, when the laptop is assembled. Users/RealComputers

It mostly depends on how youre planning to use it. In my opinion a general rule of thumb is. If its feeling really hot its better to err on the side of safety and spend the small amount of money to get the cooler. Better safe than sorry. Curtis Computing Solutions James Curtis Computing Solutions

Anchor(B)20090519 17:21:32 nbsp Hi, Hope this a considered a computer question: Does anyone know of a proxy IP/site that will support viewing vids on YouTube without the annoying notice that says I need to change my java script settings? I am in China and YouTube is blocked right now along with a bunch of other sites....Thanks! Users/jsbmeb
  There are lots of ways proxies can work. What I have done is ssh to another server I have and my ssh client program (openssh or putty) can be a socks proxy server on the localhost interface that will send the traffic through the ssh connection. Youtube does seem to work through the ssh socks proxy.  Users/NickSchmalenberger
  Reporters without Borders has a number of resources, including a http://www.rsf.org/rubrique.php3?id_rubrique542 Handbook for bloggers and cyberdissidents, and http://www.rsf.org/article.php3?id_article15013 Technical ways to get around censorship.  There is quite a bit of information out there, and if you need somebody to email you anything about how to get around it, Im sure somebody here would... although it might be illegal for you to receive it.   Users/JabberWokky Evan JabberWokky Edwards

Update: I was able to get through from what Nick told me. Thanks! jsbmeb

Anchor(C)20090521 11:27:44 nbsp Ive got an old laptop, an IBM Thinkpad T42. Its been a pretty stable commuting workhorse for me.  Recently, I slipped up and took in a ton of cruel malware. (Was streaming episodes of The Office online, followed a bad link, and left my computer for an hour since I was at work  had a zillion things popped up when I came back).  Ive exhausted myself playing with the typical Adaware/Spybot/HijackThis/Combofix, and for the first time, Im stumped. Usually I can fix it up after a few days of annoying reboots and scans and digging through my registry.  Im sure if I keep going specifically with HJT and Combofix I could clear my computer, but it seems that some of them just keep slipping back, and Im losing the willpower to devote an entire weekend to sitting in front of it.  Ive got a backup external harddrive, but Im concerned about compromising it, or my current flash drive, so I bought a new 8gig flash drive (I guess so I can keep it quarantined).  Im planning to copy over my pictures and documents, and simply reformat the laptop.  Its been over two years since I reformatted it, so its about due anyway.  However, Im concerned over how effective this will be.  If I really got something nasty, its possible (and has been known) to come back, but I doubt I have a boot sector virus, just good ole fashioned shiny new trojans and malware. Any advice? Any knowledge, anyone?  My biggest concern is that I havent backed up my documents and most importantly my pictures in a few months, so I really dont want to lose those.  As far as I can tell, Im not having virus trouble, just some really sneaky trojan/malware tidbits strongly hidden.  Users/EdWins ES
  You seem tech savy enough so, heres what you want to do: Boot to safe mode and copy your data to your flash drive. Do your reinstall, then install your usual barrage of security softwares.  Scan the flash drive before copying the data over, preferably in safe mode.  You want to work in safe mode because it loads minimal drivers and boot loaders, so it bypasses may viriis attack vectors.  It may be possible to clean your machine however. I would recomend a temp file cleaner such as CCleaner.  Malicious programs downloaded from the internet often hide in your temp files.  IE does not wipe these when it is done (One of the many reasons firefox is more secure). Users/MasonMurray
  As an aside, thats the exact model of laptop I use for casual and onsite work. They are darn near indestructible, and have decent strength in their wireless card range. Users/MasonMurray
     Re temp cleaner, I use http://www.technibble.com/repairtooloftheweekatfcleaner/ ATF Cleaner. And yes, the T42 is an incredibly well built machine, I love the durability. Ive got a stronger battery in it too, it outlasts newer laptops (even with the same battery) by hours.  Users/EdWins ES
      This is getting ridiculous! I ALSO have the extended life battery XD!  But if you are using all of those tools, and you are not able to get it booting properly, you need a reinstall.  Back up your most precious data, and hit the nuke button.  Users/MasonMurray
  20090522 08:23:33 nbsp Did you try checking your start up programs to see if it can be disabled? Go to start>run then type in mscofig and go to the start up tab. Run is located under all programs>accesories in vista. See if that helps with your pop ups. If that dosent work, then what Mason said is true and you will have to reinstall.  Users/RealComputers
     Yes, that was among the first things done. As well as doing multiple scans and digging through the reg in safe mode.  Its not really an issue with popups, that was just probably when I first got all this junk. I cant recall the malware names since its been a little bit, but the problem is two or three key trojans which not only hide in the registry, but also end up helping install other things. So I can remove everything, and eventually theyll sneak back in.  Apparently theyre among the most disingenuous and compromising. Ill be ok with reformatting, just wasnt sure if it would work out great or not. Ill just try what Mason said, do as much of it in safe mode as I can and be sure to scan through the flash drive before bringing that stuff back. users/edwins ES
       If you can find the executable file names or some other key information then you can probably find solutions online. Otherwise if you are just going to back up your data and the reformat then you dont really need to do anything in safe mode .. but it wouldnt hurt either. Also, make sure you disable all non microsoft services in the services tab on msconfig if you havent already. I usually just reformat when I get a virus just to make sure its totally dead. Users/RealComputers
  Id reformat, too.  For future reference, though, Ill add one more tool to the iteration you mentioned above: http://technet.microsoft.com/enus/sysinternals/bb963902.aspx Sysinternals Autoruns.  It reveals the same sorts of things as HJT, but not just related to IE.  Ive found it very helpful in locating where malware is hiding its .dlls and .exes. Users/MattJurach
 Final update on my situation: During the above thread, I was waiting for the mailman to deliver my USB stick and copies of the XP/install CDs.  Once they all arrived, the situation had worsened. I hadnt touched the laptop for a week, on first boot tried to log in and got a weird loop, even in safemode (would type in password, hit enter, logging in...logging out, need to reenter password).  After a few attempts, the computer refused to load! It would load up almost as normal (again, even in safemode) but by the time you get to the normal login screen, the screen would go black.  Stuck!   I booted from the XP disks, tried to repair the install, checked a few things out (userinit.exe), and basically reinstalled windows over the current install (without formatting).  Was able to finally log in.  Crappily enough, the computer now refused to recognize my USB flash drives.  When plugged in, Id see the Xyx Flash Stick Brand Blah is connected! but couldnt access it (or any of my drives).  Couldnt figure this one out either. Luckily, I had backed up some stuff more recently than I thought, so I ended up burning 6 or 7 CDs full of data/duplicates.  Then I reformatted, reinstalled, spent forever patching  it up (didnt have an original driver disk for this laptop), and finally, it was as good as new. Reinstalled a few essential programs, including an antivirus program and some antimalware stuff, transferred all my data back, etc.  It was a hell of an annoying trip. TL;DR  booted from XP cd, tried to repair install, burned data to disks, reformatted/reinstalled, good as new.  I really want to just buy a new laptop, but if I wait a year maybe a training grant will cover me for it <.<   users/edwins ES

Anchor(D)20090525 14:47:36 nbsp I have an older Dell laptop onto which someone dumped a cup of water.  This happened about a month ago, and Im pretty sure it has completely dried out now. But it still has two problems: (1) sometimes when it starts up, instead of loading Windows, it shows a black screen with white text saying Exception... and a bunch of codes.  And, (2) when it does load Windows, it is often the case that user input is slightly screwy  for example, when you click a bookmark folder in Firefox, all the bookmarks in that folder are loaded into a single page at once.  Or, when you click on something on the desktop, it selects multiple things instead of just one.  It is as though the shift or ctrl keys are always down (but note that this only happens some of the time that the computer is on, not always).  What could be the problem here?  And what can I do to fix it? Users/TheAmazingLarry
   Well... thats a bad situation. Id like to first say that fixing a computer can be alot like fixing a car.  In that the objective in fixing it is not making it new again, it is getting it to work.  That said, the main problem SEEMS to be that it is reading one of the keys as being depressed.  If you replace the keyboard it may fix that problem.  The problem with it booting could be due to water damage to the motherboard or the hard drive, although the hard drive is much less likely in that if it were damaged, the computer should never load at all, or would crash often due to lookup errors.  The motherboard problem would be around 300375 for a professional to replace whatever capacitors/chips were overloaded.  All of THAT said, I cant guarantee anything because I have not done a diagnostic on it, so it could just be a couple of loose connectors.  Users/MasonMurray
   Sounds like you need a new touchpad.Users/RealComputers
    Thanks for your assessments of this.  It looks like Im screwed on this computer because its too old and junky to put much money into repair.  The strange input behavior was from a mouse and not the touchpad, so I think thats not the problem.  The part of the problem that confuses me is that the computer sometimes works fine, but at other time exhibits the symptoms above...Users/TheAmazingLarry
     Well, the touchpad may be giving input regardless, unless you have some method for turning it off and it still does that. If you cant just press a button or flip a switch to turn it off, you can try going into control panel (assuming youre using windows) and then into the mouse options, and there should be something in there that you can disable. If you are getting screwy input from a touchpad, then disabling it using either software or hardware techniques may work...Oh, and if you dont mind having the touchpad be easily usable you could just open it up and unplug the pad from the motherboard. Of course, I would advise trying the software method first.Users/JoePomidor
  I feel like Im stating the obvious here, but have you tried using an external keyboard and/or a new external mouse to see what happens?  You should be able to get them through freecycle if you dont have any handy.   Users/CovertProfessor
    Yes, I have tried the external keyboard (but not a different mouse), but that brings up another symptom that I should have explained above.  In fact, we keep the external keyboard plugged into the laptop all the time.  If we unplug it, the laptop keyboard starts exhibiting even stranger behavior.  I tried unplugging it yesterday, for example, and the laptop keyboard started behaving as though the Fn key was being held down (so M, J, K,L, U, I, O, and P would return 0  9 rather than their alphabetic character). Users/TheAmazingLarry
     You can disable the touchpad by right clicking on computer under start, going to properties and then device manager. From there you can disable all your system components. At least this way you can be sure your touchpad is not the issue.  Users/RealComputers
  I will try disabling the touchpad to see if thats the source of the problem.  We usually do get stuck using it sometimes though, because this laptop has only two USB ports, so whenever the mouse and the wifi dongle are plugged in, the mouse has to come out to let other devices (th USB printer or camera) be used.Users/TheAmazingLarry

Anchor(E)20090616 08:04:40 nbsp Advice on buying refurb computers? Users/TheAmazingLarry
  Expanding on the above comment: Im now thinking of buying two items: a new laptop (see comment above) and a carGPS unit.  I see a lot of deals on refurbished items, but Im always a little cautious because I dont know if refurbished means refurbished like new or shorter lifespan than a new product.  Any advice on this?  My gut says that the quality of the refurb is directly proportional to the reputation of the merchant, but maybe theres more to it than that?
 I assume that a refurbished item is used but basically works. Also, items with moving parts like a hard drive or floppy drive are likely to have a shorter lifespan. Laptops have so many tiny screws and delicate cables that the more its taken apart the more chance there is for more problems. Ive replaced laptop motherboards and lcd inverter boards successfully, but Id rather not touch it if theres not an actual problem. A good technician will be well organized about the screws and stuff, but sometimes the phone rings and you drop something. Refurbished could mean something was broken but it was replaced, so then its questionable where the replacement part came from (ebay?) and what got lost or broken in the process of repairing the broken thing. Users/NickSchmalenberger
 Adding to Nicks comment, it really depends on your vendor. Apple refurbs, for example, are indistinguishable from their new computers and carry the same warranty. Some random thing at Frys.... not so much. Users/WilliamLewis

I once received a refurb Kensington trackball because the one I had bought was defective.  Except the trackball had had its finish sanded off!  So I called and complained that I had bought a new one and this was what they gave me.  So they sent me another.  Same problem.  Then one day I came home and found two new ones to replace each of the two refurbs!  Apples refurbs carry the same warranty as new products.  They take returns or defective products, replace defective parts, test them to the same standards as their new products, and ship them.  Im not positive that they will look as nice as new ones  it may be that they dont replace casings with scratches and such.  I would avoid the ones at Frys.  Seriously.  They are known for repackaging defective products in new shrink wrap.  Users/IDoNotExist

   I never buy refurbished hard or optical drives.  Nothing wrong with other components.  Note that there appears to be no quality or procedural standards for what constitutes a refurbished component.  And nobody actually disassembles a hard drive and replaces heads, bearings or platters.  The very best you can hope for is that the unit will be tested.  Worst and most probable case would be a quick repartition and a new antistatic bag.  Users/JimStewart

   Ive had Western Digital send me refurbished hard drives.  The one they sent me failed in only a few weeks.  Id avoid those... Users/IDoNotExist

Anchor(F)20090714 18:22:28 nbsp My laptop(hp ze2000) which is a few years old is very sluggish at times.  In response to a virus invasion I rebooted the hard drive, but that didnt seem to help.  I am a photographer and use Photoshop, but I store all my photos on a portable drive.  I was thinking of getting more RAM, I currently have 384 MB, but a friend of mine said that he did this and it didnt speed up his computer at all.  I also frequently get the message that my virtual memory is low and my computer is increasing the size of it (whatever that means).  What Im wondering is, is it worth it to diagnose and/or upgrade my current computer or is it time to just buy a new one?  Also, if I do get a new one, are macs better? Users/MRathje
  Which operating system are you using? Users/JasonAller

Windows Me Users/MRathje
  RAM is your physical memory and is faster than virtual memory. If your virtual memory is low your physical memory probably is also and increasing the RAM should definitely help. If it all gets taken up by viruses and stuff though its still going to be bad, so you have to be careful. By rebooting the hard drive, do you mean you erased it and reinstalled Windows? This is the simplest and surest way to get rid of viruses.  Users/NickSchmalenberger

Yes, I erased it and reinstalled Windows so I think that got rid of the viruses.  I have AVG antivirus on it as well.Users/MRathje
    I did a quick google search and came up with http://www.memorystock.com/memory/HewlettPackardPavilionze2000.html this page. It says that you can upgrade your HewlettPackard Pavilion ze2000 Laptop to up to a maximum of 2 GB Memory, the system has 2 sockets to install Memory, already with 512 MB  256 MB (removable) + 256 MB (removable) standard Memory installed.  Id recommend adding more memory; this website offers a 1 gb memory card for about $47.   What you should do, in my opinion, is contact HP or browse their their site, and find out the exact specifications of the type of memory you can put in.  384 is really low; if you were to purchase a new laptop, theres no way itd come with under 1 gigs worth of it.  Assuming youre using new software, like photoshop and other memory heavy programs, you should see a huge improvement.  Definitely browse around and price shop, but a memory upgrade is a cheaper alternative than a total upgrade (which might be an option too). users/edwins ES
      To go down the memory line further, the system requirements for http://www.adobe.com/products/photoshop/photoshop/systemreqs/ adobe photoshop CS4 states 512MB of RAM (1GB recommended).  JoePomidor linked you some some good stuff below in his comment, so again, a memory upgrade is a much cheaper option than a total upgrade, and will almost assuredly help you out.  users/edwins ES

   Cant say for sure but I would bet your computer has become sluggish because its overheating. The fan probably needs to be cleaned for dust. No virus can survive a reformat so if you did that and still had issues than its probably heat related.  Real Computers

  You computer supplements the real RAM with virtual RAM, which is on your hard drive. If you dont have enough RAM, then it will try to improve its performance by increasing the amount of virtual RAM it is using. If you dont mind gambling about $40 on trying to keep this computer going for another little while, you could check out http://www.newegg.com/Product/ProductList.aspx?SubmitENE&N2010170381%201052907864&name200Pin%20DDR%20SODIMM Newegg, get your RAM up to a good amount, and see if that helps. Another thing, I dont know to what degree this is true, but one of my friends once told me that using an irregular amount of RAM (not some exponent of 2, i.e. 256, 512, 1024) can make some computers do funny things. That could be untrue, but given the idiosyncrasies of modern  operating systems, who knows? Oh, also, if you want to try to remove some of the dust and whatnot that is probably clogging up your fan assembly, heres the http://www.frctech.com/docs/manuals/SG/hpcompaq/MSG%20HP%20Pavilion%20zd2000%20%20Presario%20M2000.pdf maintenance guide (fan stuff is on page 144). It looks like a bit of disassembly is required, so if you dont feel up to that then I would take it somewhere to be looked at by someone who does.Users/JoePomidor
      Not matching RAM in some systems or not having certain key values might mean youre not using all of the memory that is physically in the system, but I think (not having used Windows in about 11 years), that the OS should report how much it can actually see and use.  Are you checking your memory by looking in your settings, or are you assuming your memory amount because you know you put those memory cards in?.  Some video cards also borrow system memory, reducing how much you have to actually use to run programs.  jw
       You can increase the size of your virtual memory by right clicking on the hard drive, then properties, then the quota tab. However, vistual memory is very slow. I would still say it is a heat related issue if your computer used to perform much faster even with a relatively low amount of ram installed, and you have also reinstalled your OS.  Real Computers
        Virtual memory (the address used by a program can be directed to a different physical address) and paging (using mass storage to store the contents of system memory) are two distinct things. wl

Are you sure you did a clean install, or did you do a repair? Before you upgrade your memory I would suggest running these two tools. http://www.filehippo.com/download_ccleaner/ and http://www.filehippo.com/download_defraggler/. Download and run those and see if they help (make sure you run the registry cleaner portion of CCleaner). If that does not help try opening the programs that you normally run when it is sluggish and then open Windows Task Manager (Ctrl + Alt + Del). Once that opens go the the Performance tab. There should be a bar that says Memory. This is telling you how much memory your system is currently using. If the bar is high (probably more than half way) then there is a good chance that upgrading memory could help you with your problem. However, if the bar is remaining relatively low and your computer is still sluggish then adding more ram will probably not help you out much. Of course Ram is cheap so you could always try to add more ram at this point and hops it works. If after adding more ram your computer is still slow I would recommend you take it into a professional or look into buying a new PC. Users/jlc1988

Anchor(G)20101224 21:17:40 nbsp Hardware question  Users/mathewblake I have a really stupid question but its got me baffled. I went and did a super quick clean on the fan assembly. Really easy yea thats what i thought. Now from what i can tell the panels for the usbs and the monitor cord, audio cables (front and back) have no power. but when i plug in a flash drive its blinks once for a quick second and then all dark. on top of that the ethernet cable has a green light that still comes on. i have checked all the cables making sure none of them came undone in the process of the cleaning but hell i used canned air for 90% and qtips for the other 10%  so if you guys have an idea of what it could be that would be awesome.Users/mathewblake

  Did you unplug anything, by any chance?  If so, did you plug everything back in in the correct place?   Users/IDoNotExist
   Open it back up and look around the fan area.  Chances are that you loosened or unplugged a cable between the panel USBs and the motherboard.  Usually the cable is silver, about the diameter of a pencil and has a small rectangular connector on the end that mates with the motherboard.  Plug it back in (it should only go in one way) and you should be good to go.  users/jimstewart

yes i did unplug the fans there happens to be two other than the one in the power box  and im sure i plugged those back in the correct spot  i know this because the come on when i turn the computer on. I have no idea what it could be i checked the plugs on everything. But like i said above the ethernet light comes on and just the monitor mouse and keyboard dont work along with the flash on the any portable drive i plug in. Users/mathewblake

20111106 06:52:44 nbsp I have a question: After I turned my laptop off last night, all the desktop icons on my macbook have vanished except the harddrive icon when I turned it on this morning. I can find my desktop items if I go to the folder in Finder. My iTunes account also started as brandnew until I used some prompts. I would also get the pinwheel of death a lot when I was using safari until I started backing up with time machine, after I did, it got better. I was also told that the pinwheeel of death can mean my harddrive is failing. Is this weird issue with the icons vanishing also a symptom? (I know my macbook is getting old)
What do you think? Any help is much appreciated! I think Im going to go to the Mac group myg at the seniot center next week if needed....
MacBook2.1 Version: 10.5.8 Processor Name:	Intel Core 2 Duo Processor Speed:	2 GHz (Snow leopard)
  Users/jsbmeb
  The standard (easy) solutions are 1) restart, and if that doesnt work, 2) In the Finder, go to View > Arrange by > (pick anything).  If neither of those works, let me know and we can try some more complicated solutions.  I wouldnt worry about the pinwheels of death too much  we all get them from time to time.   Users/CovertProfessor
  Ok, so tried that. Not working still. I also just discovered that my chinese characters have all turned to boxes. Its like my mac had a giant fart and cant recover without more help than I can give, so I made an appt. at the apple store to see what the heck happened. I have a feeling Snow Leopard will have to reinstalled, but hopefully it will just be a couple clicks by a genius and then will be good. Thanks Professor :)

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*----------------------------------------------------------------------*
      subroutine frm_sched(xret,flist,depend_info,idxselect,nselect,
     &         init,no_skip,op_info,str_info,strmap_info,orb_info)
*----------------------------------------------------------------------*
!>     wrapper for actual scheduler routines
!!
!!     default scheduler is frm_sched2
*----------------------------------------------------------------------*
      implicit none

      include 'routes.h'
      include 'contr_times.h'

      include 'stdunit.h'
      include 'opdim.h'
      include 'mdef_operator_info.h'
      include 'def_graph.h'
      include 'def_strinf.h'
      include 'def_strmapinf.h'
      include 'def_orbinf.h'
      include 'def_contraction.h'
      include 'def_formula_item.h'
      include 'def_dependency_info.h'

      real(8), intent(out) ::
     &     xret(*)
      integer, intent(in) ::
     &     nselect, idxselect
      logical, intent(in) ::
     &     init, no_skip
      type(formula_item), intent(inout) ::
     &     flist
      type(operator_info) ::
     &     op_info
      type(strinf) ::
     &     str_info
      type(strmapinf) ::
     &     strmap_info
      type(orbinf) ::
     &     orb_info
      type(dependency_info) ::
     &     depend_info

      real(8) ::
     &     cpu, sys, wall, cpu0, sys0, wall0
c stat
      real(8) ::
     &     avg, sig
c stat

      call atim_csw(cpu0,sys0,wall0)
      cnt_op1op2(1:3) = 0d0
      cnt_dloop(1:3) = 0d0
      cnt_kernel(1:3) = 0d0
      cnt_coll1(1:2) = 0d0
      cnt_coll2(1:2) = 0d0
      cnt_dgemm(1:2) = 0d0
      cnt_scatt(1:2) = 0d0
      cnt_reo(1:2) = 0d0
      cnt_rd(1:3) = 0d0
      cnt_wr(1:3) = 0d0
      cnt_maxscr = 0
      cnt_used_reo = .false.
c dbg
      cnt_test(1:12) = 0d0
c dbg
c stat
      mm_call   = 0
      mm_dim1   = 0
      mm_dim1sq = 0
      mm_dim2   = 0
      mm_dim2sq = 0
      mm_cnt    = 0
      mm_cntsq  = 0
c stat
c more stat
      st_total = 0
      st_ex1c_z= 0
      st_ex1a_z= 0
      st_ex2c_z= 0
      st_ex2a_z= 0
      st_cntc_z= 0
      st_cnta_z= 0
      st_noop1_reo= 0
      st_noop2_reo= 0
      st_noop12_reo= 0
c more stat

      select case (irt_sched)
      case (0) !only 0 if requested by the schedule argument
        if (.not.init) call quit(1,'frm_sched',
     &                      'init option not (yet) for this scheduler')
        call frm_sched1(xret,flist,depend_info,idxselect,nselect,
     &         op_info,str_info,strmap_info,orb_info)
      case (1) ! default
        call frm_sched2(xret,flist,depend_info,idxselect,nselect,
     &         init,no_skip,op_info,str_info,strmap_info,orb_info)
      case default
        call quit(1,'frm_sched','illegal route')
      end select

      call atim_csw(cpu,sys,wall)
      if (.false..and.iprlvl.ge.10) then
        call prtim(lulog,'contraction kernel',
     &       cnt_kernel(1),cnt_kernel(2),-1d0)
c dbg
        call prtim(lulog,'check kernel     ',
     &       cnt_test(1),cnt_test(2),-1d0)
        call prtim(lulog,'     in setup    ',
     &       cnt_test(3),cnt_test(4),-1d0)
        call prtim(lulog,'     inner loop ?',
     &       cnt_test(11),cnt_test(12),-1d0)
        call prtim(lulog,'     inner loop  ',
     &       cnt_test(9),cnt_test(10),-1d0)
        call prtim(lulog,'    middle loop  ',
     &       cnt_test(7)-cnt_test(9),cnt_test(8)-cnt_test(10),-1d0)
        call prtim(lulog,'     outer loop  ',
     &       cnt_test(5)-cnt_test(7),cnt_test(6)-cnt_test(8),-1d0)
c dbg
        call prtim(lulog,'     in collect 1',
     &       cnt_coll1(1),cnt_coll1(2),-1d0)
        call prtim(lulog,'     in collect 2',
     &       cnt_coll2(1),cnt_coll2(2),-1d0)
        call prtim(lulog,'         in dgemm',
     &       cnt_dgemm(1),cnt_dgemm(2),-1d0)
        call prtim(lulog,'       in scatter',
     &       cnt_scatt(1),cnt_scatt(2),-1d0)
        if (cnt_used_reo)
     &       call prtim(lulog,'additional reord.',
     &       cnt_reo(1),cnt_reo(2),-1d0)
        call prtim(lulog,'overhead contraction 1',
     &       cnt_dloop(1)-cnt_kernel(1)-cnt_reo(1),
     &       cnt_dloop(2)-cnt_kernel(2)-cnt_reo(2),
     &       -1d0)
        call prtim(lulog,'IO read',
     &      cnt_rd(1),cnt_rd(2),-1d0)
        call prtim(lulog,'IO write',
     &      cnt_wr(1),cnt_wr(2),-1d0)
        call prtim(lulog,'overhead contraction 2',
     &       cnt_op1op2(1)-cnt_dloop(1)-cnt_rd(1)-cnt_wr(1),
     &       cnt_op1op2(2)-cnt_dloop(2)-cnt_rd(2)-cnt_wr(2),
     &       -1d0)
        call prtim(lulog,'overhead scheduler',
     &       cpu-cpu0-cnt_op1op2(1),
     &       sys-sys0-cnt_op1op2(2),
     &       -1d0)
c stat
        write(lulog,*) 'calls to dgemm: ',mm_call
        avg = dble(mm_dim1)/dble(mm_call)
        sig = sqrt(abs(avg*avg - dble(mm_dim1sq)/dble(mm_call)))
        write(lulog,'(2x,a,2f12.2)') ' avg dim1, sigma:',avg,sig
        avg = dble(mm_dim2)/dble(mm_call)
        sig = sqrt(abs(avg*avg - dble(mm_dim2sq)/dble(mm_call)))
        write(lulog,'(2x,a,2f12.2)') ' avg dim2, sigma:',avg,sig
        avg = dble(mm_cnt)/dble(mm_call)
        sig = sqrt(abs(avg*avg - dble(mm_cntsq)/dble(mm_call)))
        write(lulog,'(2x,a,2f12.2)') ' avg cnt , sigma:',avg,sig
        write(lulog,'(/2x,a,i10,a,f12.3,a)')
     &       'max. scratch: ',cnt_maxscr,' = ',
     &       dble(cnt_maxscr)*8d0/1024d0**3d0,' Gbytes'
c stat

      end if
      if (iprlvl.ge.1)
     &     call prtim(lulog,'formula evaluation',
     &     cpu-cpu0,sys-sys0,wall-wall0)
c more stat
       write(lulog,*) 'total: ',st_total
       write(lulog,*) 'EX1 C        ',st_ex1c_z
       write(lulog,*) 'EX1 A        ',st_ex1a_z
       write(lulog,*) 'EX2 C        ',st_ex2c_z
       write(lulog,*) 'EX2 A        ',st_ex2a_z
       write(lulog,*) 'CNT C        ',st_cntc_z
       write(lulog,*) 'CNT A        ',st_cnta_z
       write(lulog,*) 'NO OP1 REO   ',st_noop1_reo
       write(lulog,*) 'NO OP2 REO   ',st_noop2_reo
       write(lulog,*) 'NO FINAL REO ',st_noop12_reo
c more stat

      return
      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:
!--------------------------------------------------------------------------------
!M+
! NAME:
!       Endian_Utility
!
! PURPOSE:
!       Module containing functions to byte-swap intrinsic data types.
!
! CATEGORY:
!       Utility
!
! LANGUAGE:
!       Fortran-95
!
! CALLING SEQUENCE:
!       USE Endian_Utility
!
! MODULES:
!       Type_Kinds:   Module to hold specification kinds for variable
!                     declaration.
!
! CONTAINS:
!       Big_Endian:   Logical function that returns .TRUE. if platform
!                     is big endian.
!
!       Swap_Endian:  Function that byte-swaps input arguments.
!
! INCLUDE FILES:
!       None.
!
! EXTERNALS:
!       None.
!
! COMMON BLOCKS:
!       None.
!
! FILES ACCESSED:
!       None.
!
! SIDE EFFECTS:
!       None.
!
! RESTRICTIONS:
!       None.
!
! CREATION HISTORY:
!         Written by:   Paul van Delst, CIMSS/SSEC, 17-Mar-2000
!                       paul.vandelst@ssec.wisc.edu
!
!  Copyright (C) 2000 Paul van Delst
!
!M-
!--------------------------------------------------------------------------------

MODULE Endian_Utility

  ! ----------
  ! Module use
  ! ----------

  USE Type_Kinds

  ! ---------------------------
  ! Disable all implicit typing
  ! ---------------------------

  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------

  PRIVATE
  PUBLIC :: Big_Endian
  PUBLIC :: Swap_Endian

  ! ------------------
  ! Overload interface
  ! ------------------

  INTERFACE Swap_Endian
    MODULE PROCEDURE Swap_Short_Integer
    MODULE PROCEDURE Swap_Long_Integer
    MODULE PROCEDURE Swap_LLong_Integer
    MODULE PROCEDURE Swap_Single_Float
    MODULE PROCEDURE Swap_Double_Float
    MODULE PROCEDURE Swap_Single_Complex
    MODULE PROCEDURE Swap_Double_Complex
  END INTERFACE Swap_Endian

CONTAINS

  FUNCTION Big_Endian()

    ! ---------------
    ! Local variables
    ! ---------------

    INTEGER( Short ) :: Source = 1_Short

    ! ------------
    ! The function
    ! ------------

    LOGICAL :: Big_Endian

    ! ----------
    ! Intrinsics
    ! ----------

    INTRINSIC TRANSFER, ICHAR

    ! ----------------------------------
    ! Initialise result to little-endian
    ! ----------------------------------

    Big_Endian = .FALSE.

    ! ------------------------------------------------------------
    ! Test for "endian-ness".
    !
    ! TRANSFER( source, 'a' ) returns a result with the physical
    !   representation of the number 1, i.e. an integer, but
    !   interpreted as a character (the type of 'a' - a character,
    !   not the value, is what is important).
    !
    ! IACHAR returns the position of a character in the ASCII
    !   collating sequence associated with the kind type parameter
    !   of the character.
    ! ------------------------------------------------------------

    IF ( IACHAR( TRANSFER( Source, 'a' ) ) == 0 ) Big_Endian = .TRUE.

  END FUNCTION Big_Endian

  ELEMENTAL FUNCTION Swap_Short_Integer ( Input ) RESULT ( Output )

    ! -------------------
    ! Argument and result
    ! -------------------

    INTEGER( Short ), INTENT( IN ) :: Input
    INTEGER( Short )               :: Output

    ! ----------------
    ! Local parameters
    ! ----------------

    INTEGER, PARAMETER :: N = n_Bytes_Short

    ! ---------------
    ! Local variables
    ! ---------------

    INTEGER( Byte ), DIMENSION( N ) :: Byte_Equivalent

    ! ------------------------------------------------
    ! Byte swap the data. The extra step in the middle
    ! is necessary for those compilers that can't
    ! handle a negative strided input to TRANSFER
    ! ------------------------------------------------

    Byte_Equivalent = TRANSFER( Input, Byte_Equivalent )
    Byte_Equivalent = Byte_Equivalent( N:1:-1 )
    Output          = TRANSFER( Byte_Equivalent, Output )

  END FUNCTION Swap_Short_Integer

  ELEMENTAL FUNCTION Swap_Long_Integer ( Input ) RESULT ( Output )

    ! -------------------
    ! Argument and result
    ! -------------------

    INTEGER( Long ), INTENT( IN ) :: Input
    INTEGER( Long )               :: Output

    ! ----------------
    ! Local parameters
    ! ----------------

    INTEGER, PARAMETER :: N = n_Bytes_Long

    ! ---------------
    ! Local variables
    ! ---------------

    INTEGER( Byte ), DIMENSION( N ) :: Byte_Equivalent

    ! ------------------------------------------------
    ! Byte swap the data. The extra step in the middle
    ! is necessary for those compilers that can't
    ! handle a negative strided input to TRANSFER
    ! ------------------------------------------------

    Byte_Equivalent = TRANSFER( Input, Byte_Equivalent )
    Byte_Equivalent = Byte_Equivalent( N:1:-1 )
    Output          = TRANSFER( Byte_Equivalent, Output )

  END FUNCTION Swap_Long_Integer

  ELEMENTAL FUNCTION Swap_LLong_Integer ( Input ) RESULT ( Output )

    ! -------------------
    ! Argument and result
    ! -------------------

    INTEGER( LLong ), INTENT( IN ) :: Input
    INTEGER( LLong )               :: Output

    ! ----------------
    ! Local parameters
    ! ----------------

    INTEGER, PARAMETER :: N = n_Bytes_LLong

    ! ---------------
    ! Local variables
    ! ---------------

    INTEGER( Byte ), DIMENSION( N ) :: Byte_Equivalent

    ! ------------------------------------------------
    ! Byte swap the data. The extra step in the middle
    ! is necessary for those compilers that can't
    ! handle a negative strided input to TRANSFER
    ! ------------------------------------------------

    Byte_Equivalent = TRANSFER( Input, Byte_Equivalent )
    Byte_Equivalent = Byte_Equivalent( N:1:-1 )
    Output          = TRANSFER( Byte_Equivalent, Output )

  END FUNCTION Swap_LLong_Integer

  ELEMENTAL FUNCTION Swap_Single_Float ( Input ) RESULT ( Output )

    ! -------------------
    ! Argument and result
    ! -------------------

    REAL( Single ), INTENT( IN ) :: Input
    REAL( Single )               :: Output

    ! ----------------
    ! Local parameters
    ! ----------------

    INTEGER, PARAMETER :: N = n_Bytes_Single

    ! ---------------
    ! Local variables
    ! ---------------

    INTEGER( Byte ), DIMENSION( N ) :: Byte_Equivalent

    ! ------------------------------------------------
    ! Byte swap the data. The extra step in the middle
    ! is necessary for those compilers that can't
    ! handle a negative strided input to TRANSFER
    ! ------------------------------------------------

    Byte_Equivalent = TRANSFER( Input, Byte_Equivalent )
    Byte_Equivalent = Byte_Equivalent( N:1:-1 )
    Output          = TRANSFER( Byte_Equivalent, Output )

  END FUNCTION Swap_Single_Float

  ELEMENTAL FUNCTION Swap_Double_Float ( Input ) RESULT ( Output )

    ! -------------------
    ! Argument and result
    ! -------------------

    REAL( Double ), INTENT( IN ) :: Input
    REAL( Double )               :: Output

    ! ----------------
    ! Local parameters
    ! ----------------

    INTEGER, PARAMETER :: N = n_Bytes_Double

    ! ---------------
    ! Local variables
    ! ---------------

    INTEGER( Byte ), DIMENSION( N ) :: Byte_Equivalent

    ! ------------------------------------------------
    ! Byte swap the data. The extra step in the middle
    ! is necessary for those compilers that can't
    ! handle a negative strided input to TRANSFER
    ! ------------------------------------------------

    Byte_Equivalent = TRANSFER( Input, Byte_Equivalent )
    Byte_Equivalent = Byte_Equivalent( N:1:-1 )
    Output          = TRANSFER( Byte_Equivalent, Output )

  END FUNCTION Swap_Double_Float

  ELEMENTAL FUNCTION Swap_Single_Complex ( Input ) RESULT ( Output )

    ! -------------------
    ! Argument and result
    ! -------------------

    COMPLEX( Single ), INTENT( IN ) :: Input
    COMPLEX( Single )               :: Output

    ! ------------------
    ! Byte-swap the data
    ! ------------------

    Output = CMPLX( Swap_Endian( REAL( Input,          Single ) ), &   ! Real
                    Swap_Endian( REAL( AIMAG( Input ), Single ) ), &   ! Imaginary
                    Single )

  END FUNCTION Swap_Single_Complex

  ELEMENTAL FUNCTION Swap_Double_Complex ( Input ) RESULT ( Output )

    ! -------------------
    ! Argument and result
    ! -------------------

    COMPLEX( Double ), INTENT( IN ) :: Input
    COMPLEX( Double )               :: Output

    ! ------------------
    ! Byte-swap the data
    ! ------------------

    Output = CMPLX( Swap_Endian( REAL( Input,          Double ) ), &   ! Real
                    Swap_Endian( REAL( AIMAG( Input ), Double ) ), &   ! Imaginary
                    Double )

  END FUNCTION Swap_Double_Complex

END MODULE Endian_Utility

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
*
* This file contains a test program and functions TVTL (trivariate normal
* and t), BVTL (bivariate t), BVND (bivariate normal), STUDNT (univariate
* t), PHID (univariate normal), plus some support functions.
* The file is self contained and should compile without errors on (77)
* standard Fortran compilers. The test program demonstrates the use of
* TVTL for computing trivariate distribution values 20 test problems
* with NU = 0 (normal case), 3, 6, 9, and 12.
*
* The software is based on work described in the paper
*  "Numerical Computation of Rectangular Bivariate and Trivariate Normal
*    and t Probabilities", by
*
*          Alan Genz
*          Department of Mathematics
*          Washington State University
*          Pullman, WA 99164-3113
*          Email : alangenz@wsu.edu
*
      PROGRAM TVTST
      INTEGER I, J, NU, NT
      PARAMETER ( NT = 20 )
      DOUBLE PRECISION TVTL, LIMIT(3,NT), SIGMA(3,NT), EPS, V
      DATA ( LIMIT(I,1), I = 1, 3 ), ( SIGMA(I,1), I = 1, 3 )
     &     / .5D0, .5D0, .8D0, .1D0, .6D0, .8D0 /
      DATA ( LIMIT(I,2), I = 1, 3 ), ( SIGMA(I,2), I = 1, 3 )
     &     / -2.5D0, .5D0, .8D0, .1D0, -.6D0, -.8D0 /
      DATA ( LIMIT(I,3), I = 1, 3 ), ( SIGMA(I,3), I = 1, 3 )
     &     / 1.5D0, .5D0, .8D0, .1D0, .6D0, .8D0 /
      DATA ( LIMIT(I,4), I = 1, 3 ), ( SIGMA(I,4), I = 1, 3 )
     &     / .5D0, .5D0, .8D0, .1D0, -.6D0, -.8D0 /
      DATA ( LIMIT(I,5), I = 1, 3 ), ( SIGMA(I,5), I = 1, 3 )
     &     / .5D0, .5D0, .8D0, .1D0, -.5D0, .5D0 /
      DATA ( LIMIT(I,6), I = 1, 3 ), ( SIGMA(I,6), I = 1, 3 )
     &     / -1.5D0, .5D0, .8D0, .1D0, -.5D0, .5D0 /
      DATA ( LIMIT(I,7), I = 1, 3 ), ( SIGMA(I,7), I = 1, 3 )
     &     / 1.5D0, .5D0, .8D0, .1D0, .5D0, -.5D0 /
      DATA ( LIMIT(I,8), I = 1, 3 ), ( SIGMA(I,8), I = 1, 3 )
     &     / -.5D0, 1D0, 1.2D0, -.4D0, .2D0, .7D0 /
      DATA ( LIMIT(I,9), I = 1, 3 ), ( SIGMA(I,9), I = 1, 3 )
     &     / 1D0, 1D0, 2D0, .4D0, .8D0, .8D0 /
      DATA ( LIMIT(I,10), I = 1, 3 ), ( SIGMA(I,10), I = 1, 3 )
     &     / 1D0, 2D0, 1D0, .4D0, .8D0, .8D0 /
      DATA ( LIMIT(I,11), I = 1, 3 ), ( SIGMA(I,11), I = 1, 3 )
     &     / -2D0, -2D0, -2D0, .4D0, .8D0, .8D0 /
      DATA ( LIMIT(I,12), I = 1, 3 ), ( SIGMA(I,12), I = 1, 3 )
     *     / 1D0, 2D0, 3D0, -.998D0, -0.248D0, 0.248D0 /
      DATA ( LIMIT(I,13), I = 1, 3 ), ( SIGMA(I,13), I = 1, 3 )
     *     / -1D0, 2D0, 3D0, .25D0, 0.25D0, 0.25D0 /
      DATA ( LIMIT(I,14), I = 1, 3 ), ( SIGMA(I,14), I = 1, 3 )
     *     /  1D0, 1D0, 3D0, .998D0, 0.2482D0, 0.2487D0 /
      DATA ( LIMIT(I,15), I = 1, 3 ), ( SIGMA(I,15), I = 1, 3 )
     *     /  1D0, 1D0, 3D0, .998D0, 0.5D0, 0.5D0 /
      DATA ( LIMIT(I,16), I = 1, 3 ), ( SIGMA(I,16), I = 1, 3 )
     *     /  1D0, 1D0, 3D0, .99D0, 0.99D0, 0.99D0 /
      DATA ( LIMIT(I,17), I = 1, 3 ), ( SIGMA(I,17), I = 1, 3 )
     *     /  1D0, 2D0, 3D0, -1D0, -.99D0, .99D0 /
      DATA ( LIMIT(I,18), I = 1, 3 ), ( SIGMA(I,18), I = 1, 3 )
     *     /  1D0, 2D0, 3D0, 1D0, -.99D0, -.99D0 /
      DATA ( LIMIT(I,19), I = 1, 3 ), ( SIGMA(I,19), I = 1, 3 )
     *     /  1D0, -1D0, 1D0, .998D0, -0.2482D0, -0.2482D0 /
      DATA ( LIMIT(I,NT), I = 1, 3 ), ( SIGMA(I,NT), I = 1, 3 )
     *     /  1D0, -1D0, 2D0, .99992D0, 0.64627D0, 0.63975D0 /
      EPS = 1D-6
      PRINT '(''      Trivariate t Test with EPS ='', E10.1)', EPS
      DO NU = 0, 12, 3
         PRINT '(''NU   B1   B2   B3    R21      R31      R32    TVT'')'
         DO J = 1, NT
            V = TVTL( NU, LIMIT(1,J), SIGMA(1,J), EPS )
            PRINT '(I2,3F5.1,3F9.5,F13.10)', NU,
     &           ( LIMIT(I,J), I = 1, 3 ), ( SIGMA(I,J), I = 1, 3 ), V
         END DO
      END DO
      END
*
      DOUBLE PRECISION FUNCTION TVTL( NU, H, R, EPSI )
*
*     A function for computing trivariate normal and t-probabilities.
*     This function uses algorithms developed from the ideas
*     described in the papers:
*       R.L. Plackett, Biometrika 41(1954), pp. 351-360.
*       Z. Drezner, Math. Comp. 62(1994), pp. 289-294.
*     with adaptive integration from (0,0,1) to (0,0,r23) to R.
*
*      Calculate the probability that X(I) < H(I), for I = 1,2,3
*    NU   INTEGER degrees of freedom; use NU = 0 for normal cases.
*    H    REAL array of uppoer limits for probability distribution
*    R    REAL array of three correlation coefficients, R should
*         contain the lower left portion of the correlation matrix r.
*         R should contains the values r21, r31, r23 in that order.
*   EPSI  REAL required absolute accuracy; maximum accuracy for most
*          computations is approximately 1D-14
*
*    The software is based on work described in the paper
*     "Numerical Computation of Rectangular Bivariate and Trivariate
*      Normal and t Probabilities", by the code author:
*
*       Alan Genz
*       Department of Mathematics
*       Washington State University
*       Pullman, WA 99164-3113
*       Email : alangenz@wsu.edu
*
*
* Copyright (C) 2013, Alan Genz,  All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided the following conditions are met:
*   1. Redistributions of source code must retain the above copyright
*      notice, this list of conditions and the following disclaimer.
*   2. Redistributions in binary form must reproduce the above copyright
*      notice, this list of conditions and the following disclaimer in
*      the documentation and/or other materials provided with the
*      distribution.
*   3. The contributor name(s) may not be used to endorse or promote
*      products derived from this software without specific prior
*      written permission.
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
* OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
* TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF USE
* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
      EXTERNAL TVTMFN
      INTEGER NU, NUC
      DOUBLE PRECISION H(3), H1, H2, H3, R(3), R12, R13, R23, EPSI
      DOUBLE PRECISION ONE, ZRO, EPS, ZROS(3), HS(3), TVT
      DOUBLE PRECISION RUA, RUB, AR, RUC, PT, BVTL, PHID, ADONET
      PARAMETER ( ZRO = 0, ONE = 1 )
      COMMON /TVTMBK/ H1, H2, H3, R23, RUA, RUB, AR, RUC, NUC
      EPS = MAX( 1D-14, EPSI )
      PT = ASIN(ONE)
      NUC = NU
      H1 = H(1)
      H2 = H(2)
      H3 = H(3)
      R12 = R(1)
      R13 = R(2)
      R23 = R(3)
*
*     Sort R's and check for special cases
*
      IF ( ABS(R12) .GT. ABS(R13) ) THEN
         H2 = H3
         H3 = H(2)
         R12 = R13
         R13 = R(1)
      END IF
      IF ( ABS(R13) .GT. ABS(R23) ) THEN
         H1 = H2
         H2 = H(1)
         R23 = R13
         R13 = R(3)
      END IF
      TVT = 0
      IF ( ABS(H1) + ABS(H2) + ABS(H3) .LT. EPS ) THEN
         TVT = ( 1 + ( ASIN(R12) + ASIN(R13) + ASIN(R23) )/PT )/8
      ELSE IF ( NU .LT. 1 .AND. ABS(R12) + ABS(R13) .LT. EPS ) THEN
         TVT = PHID(H1)*BVTL( NU, H2, H3, R23 )
      ELSE IF ( NU .LT. 1 .AND. ABS(R13) + ABS(R23) .LT. EPS ) THEN
         TVT = PHID(H3)*BVTL( NU, H1, H2, R12 )
      ELSE IF ( NU .LT. 1 .AND. ABS(R12) + ABS(R23) .LT. EPS ) THEN
         TVT = PHID(H2)*BVTL( NU, H1, H3, R13 )
      ELSE IF ( 1 - R23 .LT. EPS ) THEN
         TVT = BVTL( NU, H1, MIN( H2, H3 ), R12 )
      ELSE IF ( R23 + 1 .LT. EPS ) THEN
         IF  ( H2 .GT. -H3 )
     &        TVT = BVTL( NU, H1, H2, R12 ) - BVTL( NU, H1, -H3, R12 )
      ELSE
*
*        Compute singular TVT value
*
         IF ( NU .LT. 1 ) THEN
            TVT = BVTL( NU, H2, H3, R23 )*PHID(H1)
         ELSE IF ( R23 .GE. 0 ) THEN
            TVT = BVTL( NU, H1, MIN( H2, H3 ), ZRO )
         ELSE IF ( H2 .GT. -H3 ) THEN
            TVT = BVTL( NU, H1, H2, ZRO ) - BVTL( NU, H1, -H3, ZRO )
         END IF
*
*        Use numerical integration to compute probability
*
*
         RUA = ASIN( R12 )
         RUB = ASIN( R13 )
         AR = ASIN( R23)
         RUC = SIGN( PT, AR ) - AR
         TVT = TVT + ADONET( TVTMFN, ZRO, ONE, EPS )/( 4*PT )
      END IF
      TVTL = MAX( ZRO, MIN( TVT, ONE ) )
      END
*
      DOUBLE PRECISION FUNCTION TVTMFN( X )
*
*     Computes Plackett formula integrands
*
      INTEGER NU
      DOUBLE PRECISION X, H1, H2, H3, R23, RUA, RUB, AR, RUC
      DOUBLE PRECISION R12, RR2, R13, RR3, R, RR, ZRO, PNTGND
      PARAMETER ( ZRO = 0 )
      COMMON /TVTMBK/ H1, H2, H3, R23, RUA, RUB, AR, RUC, NU
      TVTMFN = 0
      CALL SINCS( RUA*X, R12, RR2 )
      CALL SINCS( RUB*X, R13, RR3 )
      IF ( ABS(RUA) .GT. 0 )
     &     TVTMFN = TVTMFN + RUA*PNTGND( NU, H1,H2,H3, R13,R23,R12,RR2 )
      IF ( ABS(RUB) .GT. 0 )
     &     TVTMFN = TVTMFN + RUB*PNTGND( NU, H1,H3,H2, R12,R23,R13,RR3 )
      IF ( NU .GT. 0 ) THEN
         CALL SINCS( AR + RUC*X, R, RR )
         TVTMFN = TVTMFN - RUC*PNTGND( NU, H2, H3, H1, ZRO, ZRO, R, RR )
      END IF
      END
*
      SUBROUTINE SINCS( X, SX, CS )
*
*     Computes SIN(X), COS(X)^2, with series approx. for |X| near PI/2
*
      DOUBLE PRECISION X, SX, CS, PT, EE
      PARAMETER ( PT = 1.57079632679489661923132169163975D0 )
      EE = ( PT - ABS(X) )**2
      IF ( EE .LT. 5D-5 ) THEN
         SX = SIGN( 1 - EE*( 1 - EE/12 )/2, X )
         CS = EE*( 1 - EE*( 1 - 2*EE/15 )/3 )
      ELSE
         SX = SIN(X)
         CS = 1 - SX*SX
      END IF
      END
*
      DOUBLE PRECISION FUNCTION PNTGND( NU, BA, BB, BC, RA, RB, R, RR )
*
*     Computes Plackett formula integrand
*
      INTEGER NU
      DOUBLE PRECISION BA, BB, BC, RA, RB, R, RR
      DOUBLE PRECISION DT, FT, BT, PHID, STUDNT
      PNTGND = 0
      DT = RR*( RR - ( RA - RB )**2 - 2*RA*RB*( 1 - R ) )
      IF ( DT .GT. 0 ) THEN
         BT = ( BC*RR + BA*( R*RB - RA ) + BB*( R*RA -RB ) )/SQRT(DT)
         FT = ( BA - R*BB )**2/RR + BB*BB
         IF ( NU .LT. 1 ) THEN
            IF ( BT .GT. -10 .AND. FT .LT. 100 ) THEN
               PNTGND = EXP( -FT/2 )
               IF ( BT .LT. 10 ) PNTGND = PNTGND*PHID(BT)
            END IF
         ELSE
            FT = SQRT( 1 + FT/NU )
            PNTGND = STUDNT( NU, BT/FT )/FT**NU
         END IF
      END IF
      END
*
      DOUBLE PRECISION FUNCTION ADONET( F, A, B, TOL )
*
*     One Dimensional Globally Adaptive Integration Function
*
      EXTERNAL F
      DOUBLE PRECISION F, A, B, TOL
      INTEGER NL, I, IM, IP
      PARAMETER ( NL = 100 )
      DOUBLE PRECISION EI(NL), AI(NL), BI(NL), FI(NL), FIN, ERR, KRNRDT
      COMMON /ABLK/ ERR, IM
      AI(1) = A
      BI(1) = B
      ERR = 1
      IP = 1
      IM = 1
      DO WHILE ( 4*ERR .GT. TOL .AND. IM .LT. NL )
         IM = IM + 1
         BI(IM) = BI(IP)
         AI(IM) = ( AI(IP) + BI(IP) )/2
         BI(IP) = AI(IM)
         FI(IP) = KRNRDT( AI(IP), BI(IP), F, EI(IP) )
         FI(IM) = KRNRDT( AI(IM), BI(IM), F, EI(IM) )
         ERR = 0
         FIN = 0
         DO I = 1, IM
            IF ( EI(I) .GT. EI(IP) ) IP = I
            FIN = FIN + FI(I)
            ERR = ERR + EI(I)**2
         END DO
         ERR = SQRT( ERR )
      END DO
      ADONET = FIN
      END
*
      DOUBLE PRECISION FUNCTION KRNRDT( A, B, F, ERR )
*
*     Kronrod Rule
*
      DOUBLE PRECISION A, B, ERR, T, CEN, F, FC, WID, RESG, RESK
*
*        The abscissae and weights are given for the interval (-1,1);
*        only positive abscissae and corresponding weights are given.
*
*        XGK    - abscissae of the 2N+1-point Kronrod rule:
*                 XGK(2), XGK(4), ...  N-point Gauss rule abscissae;
*                 XGK(1), XGK(3), ...  optimally added abscissae.
*        WGK    - weights of the 2N+1-point Kronrod rule.
*        WG     - weights of the N-point Gauss rule.
*
      INTEGER J, N
      PARAMETER ( N = 11 )
      DOUBLE PRECISION WG(0:(N+1)/2), WGK(0:N), XGK(0:N)
      SAVE WG, WGK, XGK
      DATA WG( 0)/ 0.2729250867779007D+00/
      DATA WG( 1)/ 0.5566856711617449D-01/
      DATA WG( 2)/ 0.1255803694649048D+00/
      DATA WG( 3)/ 0.1862902109277352D+00/
      DATA WG( 4)/ 0.2331937645919914D+00/
      DATA WG( 5)/ 0.2628045445102478D+00/
*
      DATA XGK( 0)/ 0.0000000000000000D+00/
      DATA XGK( 1)/ 0.9963696138895427D+00/
      DATA XGK( 2)/ 0.9782286581460570D+00/
      DATA XGK( 3)/ 0.9416771085780681D+00/
      DATA XGK( 4)/ 0.8870625997680953D+00/
      DATA XGK( 5)/ 0.8160574566562211D+00/
      DATA XGK( 6)/ 0.7301520055740492D+00/
      DATA XGK( 7)/ 0.6305995201619651D+00/
      DATA XGK( 8)/ 0.5190961292068118D+00/
      DATA XGK( 9)/ 0.3979441409523776D+00/
      DATA XGK(10)/ 0.2695431559523450D+00/
      DATA XGK(11)/ 0.1361130007993617D+00/
*
      DATA WGK( 0)/ 0.1365777947111183D+00/
      DATA WGK( 1)/ 0.9765441045961290D-02/
      DATA WGK( 2)/ 0.2715655468210443D-01/
      DATA WGK( 3)/ 0.4582937856442671D-01/
      DATA WGK( 4)/ 0.6309742475037484D-01/
      DATA WGK( 5)/ 0.7866457193222764D-01/
      DATA WGK( 6)/ 0.9295309859690074D-01/
      DATA WGK( 7)/ 0.1058720744813894D+00/
      DATA WGK( 8)/ 0.1167395024610472D+00/
      DATA WGK( 9)/ 0.1251587991003195D+00/
      DATA WGK(10)/ 0.1312806842298057D+00/
      DATA WGK(11)/ 0.1351935727998845D+00/
*
*           Major variables
*
*           CEN  - mid point of the interval
*           WID  - half-length of the interval
*           RESG - result of the N-point Gauss formula
*           RESK - result of the 2N+1-point Kronrod formula
*
*           Compute the 2N+1-point Kronrod approximation to
*            the integral, and estimate the absolute error.
*
      WID = ( B - A )/2
      CEN = ( B + A )/2
      FC = F(CEN)
      RESG = FC*WG(0)
      RESK = FC*WGK(0)
      DO J = 1, N
         T = WID*XGK(J)
         FC = F( CEN - T ) + F( CEN + T )
         RESK = RESK + WGK(J)*FC
         IF( MOD( J, 2 ) .EQ. 0 ) RESG = RESG + WG(J/2)*FC
      END DO
      KRNRDT = WID*RESK
      ERR = ABS( WID*( RESK - RESG ) )
      END
*
      DOUBLE PRECISION FUNCTION STUDNT( NU, T )
*
*     Student t Distribution Function
*
*                       T
*         STUDNT = C   I  ( 1 + y*y/NU )**( -(NU+1)/2 ) dy
*                   NU -INF
*
      INTEGER NU, J
      DOUBLE PRECISION T, ZRO, ONE, PI, PHID
      DOUBLE PRECISION CSSTHE, SNTHE, POLYN, TT, TS, RN
      PARAMETER ( ZRO = 0, ONE = 1 )
      PI = ACOS(-ONE)
      IF ( NU .LT. 1 ) THEN
         STUDNT = PHID( T )
      ELSE IF ( NU .EQ. 1 ) THEN
         STUDNT = ( 1 + 2*ATAN(T)/PI )/2
      ELSE IF ( NU .EQ. 2 ) THEN
         STUDNT = ( 1 + T/SQRT( 2 + T*T ))/2
      ELSE
         TT = T*T
         CSSTHE = 1/( 1 + TT/NU )
         POLYN = 1
         DO J = NU-2, 2, -2
            POLYN = 1 + ( J - 1 )*CSSTHE*POLYN/J
         END DO
         IF ( MOD( NU, 2 ) .EQ. 1 ) THEN
            RN = NU
            TS = T/SQRT(RN)
            STUDNT = ( 1 + 2*( ATAN(TS) + TS*CSSTHE*POLYN )/PI )/2
         ELSE
            SNTHE = T/SQRT( NU + TT )
            STUDNT = ( 1 + SNTHE*POLYN )/2
         END IF
         STUDNT = MAX( ZRO, MIN( STUDNT, ONE ) )
      ENDIF
      END
*
      DOUBLE PRECISION FUNCTION BVTL( NU, DH, DK, R )
*
*     A function for computing bivariate t probabilities.
*
*       Alan Genz
*       Department of Mathematics
*       Washington State University
*       Pullman, WA 99164-3113
*       Email : alangenz@wsu.edu
*
*    This function is based on the method described by
*        Dunnett, C.W. and M. Sobel, (1954),
*        A bivariate generalization of Student's t-distribution
*        with tables for certain special cases,
*        Biometrika 41, pp. 153-169.
*
* BVTL - calculate the probability that X < DH and Y < DK.
*
* parameters
*
*   NU number of degrees of freedom
*   DH 1st lower integration limit
*   DK 2nd lower integration limit
*   R   correlation coefficient
*
      INTEGER NU, J, HS, KS
      DOUBLE PRECISION DH, DK, R
      DOUBLE PRECISION TPI, PI, ORS, HRK, KRH, BVT, SNU, BVND, STUDNT
      DOUBLE PRECISION GMPH, GMPK, XNKH, XNHK, QHRK, HKN, HPK, HKRN
      DOUBLE PRECISION BTNCKH, BTNCHK, BTPDKH, BTPDHK, ONE, EPS
      PARAMETER ( ONE = 1, EPS = 1D-15 )
      IF ( NU .LT. 1 ) THEN
         BVTL = BVND( -DH, -DK, R )
      ELSE IF ( 1 - R .LE. EPS ) THEN
            BVTL = STUDNT( NU, MIN( DH, DK ) )
      ELSE IF ( R + 1  .LE. EPS ) THEN
         IF ( DH .GT. -DK )  THEN
            BVTL = STUDNT( NU, DH ) - STUDNT( NU, -DK )
         ELSE
            BVTL = 0
         END IF
      ELSE
         PI = ACOS(-ONE)
         TPI = 2*PI
         SNU = NU
         SNU = SQRT(SNU)
         ORS = 1 - R*R
         HRK = DH - R*DK
         KRH = DK - R*DH
         IF ( ABS(HRK) + ORS .GT. 0 ) THEN
            XNHK = HRK**2/( HRK**2 + ORS*( NU + DK**2 ) )
            XNKH = KRH**2/( KRH**2 + ORS*( NU + DH**2 ) )
         ELSE
            XNHK = 0
            XNKH = 0
         END IF
         HS = SIGN( ONE, DH - R*DK )
         KS = SIGN( ONE, DK - R*DH )
         IF ( MOD( NU, 2 ) .EQ. 0 ) THEN
            BVT = ATAN2( SQRT(ORS), -R )/TPI
            GMPH = DH/SQRT( 16*( NU + DH**2 ) )
            GMPK = DK/SQRT( 16*( NU + DK**2 ) )
            BTNCKH = 2*ATAN2( SQRT( XNKH ), SQRT( 1 - XNKH ) )/PI
            BTPDKH = 2*SQRT( XNKH*( 1 - XNKH ) )/PI
            BTNCHK = 2*ATAN2( SQRT( XNHK ), SQRT( 1 - XNHK ) )/PI
            BTPDHK = 2*SQRT( XNHK*( 1 - XNHK ) )/PI
            DO J = 1, NU/2
               BVT = BVT + GMPH*( 1 + KS*BTNCKH )
               BVT = BVT + GMPK*( 1 + HS*BTNCHK )
               BTNCKH = BTNCKH + BTPDKH
               BTPDKH = 2*J*BTPDKH*( 1 - XNKH )/( 2*J + 1 )
               BTNCHK = BTNCHK + BTPDHK
               BTPDHK = 2*J*BTPDHK*( 1 - XNHK )/( 2*J + 1 )
               GMPH = GMPH*( 2*J - 1 )/( 2*J*( 1 + DH**2/NU ) )
               GMPK = GMPK*( 2*J - 1 )/( 2*J*( 1 + DK**2/NU ) )
            END DO
         ELSE
            QHRK = SQRT( DH**2 + DK**2 - 2*R*DH*DK + NU*ORS )
            HKRN = DH*DK + R*NU
            HKN = DH*DK - NU
            HPK = DH + DK
            BVT = ATAN2( -SNU*( HKN*QHRK + HPK*HKRN ),
     &                          HKN*HKRN-NU*HPK*QHRK )/TPI
            IF ( BVT .LT. -EPS ) BVT = BVT + 1
            GMPH = DH/( TPI*SNU*( 1 + DH**2/NU ) )
            GMPK = DK/( TPI*SNU*( 1 + DK**2/NU ) )
            BTNCKH = SQRT( XNKH )
            BTPDKH = BTNCKH
            BTNCHK = SQRT( XNHK )
            BTPDHK = BTNCHK
            DO J = 1, ( NU - 1 )/2
               BVT = BVT + GMPH*( 1 + KS*BTNCKH )
               BVT = BVT + GMPK*( 1 + HS*BTNCHK )
               BTPDKH = ( 2*J - 1 )*BTPDKH*( 1 - XNKH )/( 2*J )
               BTNCKH = BTNCKH + BTPDKH
               BTPDHK = ( 2*J - 1 )*BTPDHK*( 1 - XNHK )/( 2*J )
               BTNCHK = BTNCHK + BTPDHK
               GMPH = 2*J*GMPH/( ( 2*J + 1 )*( 1 + DH**2/NU ) )
               GMPK = 2*J*GMPK/( ( 2*J + 1 )*( 1 + DK**2/NU ) )
            END DO
         END IF
         BVTL = BVT
      END IF
*     END BVTL
      END
*
      DOUBLE PRECISION FUNCTION PHID(Z)
*
*     Normal distribution probabilities accurate to 1d-15.
*     Reference: J.L. Schonfelder, Math Comp 32(1978), pp 1232-1240.
*
      INTEGER I, IM
      DOUBLE PRECISION A(0:43), BM, B, BP, P, RTWO, T, XA, Z
      PARAMETER( RTWO = 1.414213562373095048801688724209D0, IM = 24 )
      SAVE A
      DATA ( A(I), I = 0, 43 )/
     &    6.10143081923200417926465815756D-1,
     &   -4.34841272712577471828182820888D-1,
     &    1.76351193643605501125840298123D-1,
     &   -6.0710795609249414860051215825D-2,
     &    1.7712068995694114486147141191D-2,
     &   -4.321119385567293818599864968D-3,
     &    8.54216676887098678819832055D-4,
     &   -1.27155090609162742628893940D-4,
     &    1.1248167243671189468847072D-5, 3.13063885421820972630152D-7,
     &   -2.70988068537762022009086D-7, 3.0737622701407688440959D-8,
     &    2.515620384817622937314D-9, -1.028929921320319127590D-9,
     &    2.9944052119949939363D-11, 2.6051789687266936290D-11,
     &   -2.634839924171969386D-12, -6.43404509890636443D-13,
     &    1.12457401801663447D-13, 1.7281533389986098D-14,
     &   -4.264101694942375D-15, -5.45371977880191D-16,
     &    1.58697607761671D-16, 2.0899837844334D-17,
     &   -5.900526869409D-18, -9.41893387554D-19, 2.14977356470D-19,
     &    4.6660985008D-20, -7.243011862D-21, -2.387966824D-21,
     &    1.91177535D-22, 1.20482568D-22, -6.72377D-25, -5.747997D-24,
     &   -4.28493D-25, 2.44856D-25, 4.3793D-26, -8.151D-27, -3.089D-27,
     &    9.3D-29, 1.74D-28, 1.6D-29, -8.0D-30, -2.0D-30 /
*
      XA = ABS(Z)/RTWO
      IF ( XA .GT. 100 ) THEN
         P = 0
      ELSE
         T = ( 8*XA - 30 ) / ( 4*XA + 15 )
         BM = 0
         B  = 0
         DO I = IM, 0, -1
            BP = B
            B  = BM
            BM = T*B - BP  + A(I)
         END DO
         P = EXP( -XA*XA )*( BM - BP )/4
      END IF
      IF ( Z .GT. 0 ) P = 1 - P
      PHID = P
      END
*
      DOUBLE PRECISION FUNCTION BVND( DH, DK, R )
*
*     A function for computing bivariate normal probabilities.
*
*       Alan Genz
*       Department of Mathematics
*       Washington State University
*       Pullman, WA 99164-3113
*       Email : alangenz@wsu.edu
*
*    This function is based on the method described by
*        Drezner, Z and G.O. Wesolowsky, (1989),
*        On the computation of the bivariate normal integral,
*        Journal of Statist. Comput. Simul. 35, pp. 101-107,
*    with major modifications for double precision, and for |R| close to 1.
*
* BVND calculates the probability that X > DH and Y > DK.
*      Note: Prob( X < DH, Y < DK ) = BVND( -DH, -DK, R ).
*
* Parameters
*
*   DH  DOUBLE PRECISION, integration limit
*   DK  DOUBLE PRECISION, integration limit
*   R   DOUBLE PRECISION, correlation coefficient
*
      DOUBLE PRECISION DH, DK, R, TWOPI
      INTEGER I, IS, LG, NG
      PARAMETER ( TWOPI = 6.283185307179586D0 )
      DOUBLE PRECISION X(10,3), W(10,3), AS, A, B, C, D, RS, XS, BVN
      DOUBLE PRECISION PHID, SN, ASR, H, K, BS, HS, HK
*     Gauss Legendre Points and Weights, N =  6
      DATA ( W(I,1), X(I,1), I = 1,3) /
     &  0.1713244923791705D+00,-0.9324695142031522D+00,
     &  0.3607615730481384D+00,-0.6612093864662647D+00,
     &  0.4679139345726904D+00,-0.2386191860831970D+00/
*     Gauss Legendre Points and Weights, N = 12
      DATA ( W(I,2), X(I,2), I = 1,6) /
     &  0.4717533638651177D-01,-0.9815606342467191D+00,
     &  0.1069393259953183D+00,-0.9041172563704750D+00,
     &  0.1600783285433464D+00,-0.7699026741943050D+00,
     &  0.2031674267230659D+00,-0.5873179542866171D+00,
     &  0.2334925365383547D+00,-0.3678314989981802D+00,
     &  0.2491470458134029D+00,-0.1252334085114692D+00/
*     Gauss Legendre Points and Weights, N = 20
      DATA ( W(I,3), X(I,3), I = 1, 10 ) /
     &  0.1761400713915212D-01,-0.9931285991850949D+00,
     &  0.4060142980038694D-01,-0.9639719272779138D+00,
     &  0.6267204833410906D-01,-0.9122344282513259D+00,
     &  0.8327674157670475D-01,-0.8391169718222188D+00,
     &  0.1019301198172404D+00,-0.7463319064601508D+00,
     &  0.1181945319615184D+00,-0.6360536807265150D+00,
     &  0.1316886384491766D+00,-0.5108670019508271D+00,
     &  0.1420961093183821D+00,-0.3737060887154196D+00,
     &  0.1491729864726037D+00,-0.2277858511416451D+00,
     &  0.1527533871307259D+00,-0.7652652113349733D-01/
      SAVE X, W
      IF ( ABS(R) .LT. 0.3 ) THEN
         NG = 1
         LG = 3
      ELSE IF ( ABS(R) .LT. 0.75 ) THEN
         NG = 2
         LG = 6
      ELSE
         NG = 3
         LG = 10
      ENDIF
      H = DH
      K = DK
      HK = H*K
      BVN = 0
      IF ( ABS(R) .LT. 0.925 ) THEN
         IF ( ABS(R) .GT. 0 ) THEN
            HS = ( H*H + K*K )/2
            ASR = ASIN(R)
            DO I = 1, LG
               DO IS = -1, 1, 2
                  SN = SIN( ASR*(  IS*X(I,NG) + 1 )/2 )
                  BVN = BVN + W(I,NG)*EXP( ( SN*HK-HS )/( 1-SN*SN ) )
               END DO
            END DO
            BVN = BVN*ASR/( 2*TWOPI )
         ENDIF
         BVN = BVN + PHID(-H)*PHID(-K)
      ELSE
         IF ( R .LT. 0 ) THEN
            K = -K
            HK = -HK
         ENDIF
         IF ( ABS(R) .LT. 1 ) THEN
            AS = ( 1 - R )*( 1 + R )
            A = SQRT(AS)
            BS = ( H - K )**2
            C = ( 4 - HK )/8
            D = ( 12 - HK )/16
            ASR = -( BS/AS + HK )/2
            IF ( ASR .GT. -100 ) BVN = A*EXP(ASR)
     &             *( 1 - C*( BS - AS )*( 1 - D*BS/5 )/3 + C*D*AS*AS/5 )
            IF ( -HK .LT. 100 ) THEN
               B = SQRT(BS)
               BVN = BVN - EXP( -HK/2 )*SQRT(TWOPI)*PHID(-B/A)*B
     &                    *( 1 - C*BS*( 1 - D*BS/5 )/3 )
            ENDIF
            A = A/2
            DO I = 1, LG
               DO IS = -1, 1, 2
                  XS = ( A*(  IS*X(I,NG) + 1 ) )**2
                  RS = SQRT( 1 - XS )
                  ASR = -( BS/XS + HK )/2
                  IF ( ASR .GT. -100 ) THEN
                     BVN = BVN + A*W(I,NG)*EXP( ASR )
     &                    *( EXP( -HK*XS/( 2*( 1 + RS )**2 ) )/RS
     &                    - ( 1 + C*XS*( 1 + D*XS ) ) )
                  END IF
               END DO
            END DO
            BVN = -BVN/TWOPI
         ENDIF
         IF ( R .GT. 0 ) THEN
            BVN =  BVN + PHID( -MAX( H, K ) )
         ELSE
            BVN = -BVN
            IF ( K .GT. H ) THEN
               IF ( H .LT. 0 ) THEN
                  BVN = BVN + PHID(K)  - PHID(H)
               ELSE
                  BVN = BVN + PHID(-H) - PHID(-K)
               ENDIF
            ENDIF
         ENDIF
      ENDIF
      BVND = BVN
      END

max_stars_repo_name: patel-zeel/lab

max_issues_repo_name: patel-zeel/lab

max_forks_repo_name: patel-zeel/lab


text:
module Cauchylowrank
  implicit none

contains

!!!!!!
    subroutine rrluCauchy(D,F,A,V,tol,U,W,PL,PU,M,N,Rk)
!
! .. Scalar parameter ..
    DOUBLE PRECISION  TOL
    INTEGER           M,N,Rk
! .. Array parameters ..
    DOUBLE PRECISION D(*),F(*),A(*),V(*),U(*),W(*)
    INTEGER          PL(*),PU(*)
!
! Purpose
! =========
! This routine computes a Rank Revealing Schur Complement, RRSC, for a Cauchy-like matrix
! via structured matrices. This matrix C has dimensions M-by-N, with generators D,F,A and V.
! C = (A(i)*V(j) / (D(i)-F(j)) )_{i,j}.
!
! ============
! Written by Sheng-Guo Li, On Sept. 18th, 2012
! ============

! .. Local Scalars
    double precision rho, Amax, Ukk, zero, one, Nmax
    integer          nswap,j,k,mn,prd,flgl,ii,jj
    parameter        (zero = 0.0D0, one=1.0D0)

!  .. Intrinsic Functions ..
    intrinsic    max, abs, maxval,mod

    U(1:M)  = A(1:M)
    Rk = min(M,N)
    mn = Rk
    PL(1:M) = (/ (j, j=1,M) /)
    PU(1:N) = (/ (j, j=1,N) /)
    nswap = 0
    Amax = zero
    rho = 1.12D0
    prd = 10

    do k = 1, mn
       call CauchyPivtG(D,F,U,V,PL,PU,k,A,M,N)

       Ukk = u(k)*V(k)/( D(k)-F(k))
       Amax = max(Amax,abs(Ukk) )
       if (abs(Ukk) .lt. tol*Amax ) then  ! first step converged
          call CauchyPivtG_CP(D,F,U,V,PL,PU,k,A,M,N)
          Ukk = u(k)*V(k)/( D(k)-F(k))
          Amax = max(Amax,abs(Ukk) )
          if (abs(Ukk) .lt. tol*Amax ) then  ! final converged
             Rk = k -1
             exit
          end if
       end if

       U(k+1:M) = U(k+1:M)* (D(k+1:M)-D(k)) / (D(k+1:M)-F(k))
       V(k+1:N) = V(k+1:N)* (F(k+1:N)-F(k)) / (F(k+1:N)-D(k))
       do j = 1, k-1
          W(j) = W(j) * ( (F(k)-D(j))/ (D(k)-D(j)) )
       end do
       W(k) = (D(k)-F(k))/A(k)
       do j = 1, k-1
          W(k) = W(k) * (F(j)-D(k)) / ( D(j)-D(k) )
       end do

       ! swap
       flgl = mod(k,prd)
       do while(flgl .lt. 1)
          flgl = 1
          call searchMax2(U(k+1),W,D(k+1),D,M-k,k,ii,jj,Nmax)

          if(Nmax .gt. rho) then
             nswap = nswap + 1
             flgl = 0
             jj = jj + k
             V(k+1:N)    = V(k+1:N) * ( (F(k+1:N)-D(ii)) / (F(k+1:N)-D(jj)) )
             U(k+1:jj-1) = U(k+1:jj-1) * ( (D(k+1:jj-1)-D(jj)) / (D(k+1:jj-1)-D(ii)) )
             U(jj+1:M)   = U(jj+1:M) * ( (D(jj+1:M)-D(jj)) / (D(jj+1:M)-D(ii)) )
             W(1:ii-1)   = W(1:ii-1) * ( (D(1:ii-1)-D(ii)) / (D(1:ii-1)-D(jj)) )
             W(ii+1:k)   = W(ii+1:k) * ( (D(ii+1:k)-D(ii)) / (D(ii+1:k)-D(jj)) )
             U(jj)       = A(ii) * ( (D(ii)-D(jj)) / (D(ii)-F(ii)) )
             W(ii)       = (D(jj)-F(ii)) / A(jj)
             do j = 1, ii-1
                U(jj) = U(jj) * ( (D(ii)-D(j)) / (D(ii)-F(j)) )
                W(ii) = W(ii) * ( (F(j)-D(jj)) / (D(j)-D(jj)) )
             end do
             do j = ii+1, k
                U(jj) = U(jj) * ( (D(ii)-D(j)) / (D(ii)-F(j)) )
                W(ii) = W(ii) * ( (F(j)-D(jj)) / (D(j)-D(jj)) )
             end do
             call iswap(PL,ii,jj)
             call dswap(D,ii,jj)
             call dswap(A,ii,jj)
!             write(*,*) 'swap once '
          end if ! Nmax

       end do ! while

    end do ! main loop

  end subroutine rrluCauchy

!!!!!!
    subroutine searchMax2(U,W,D2,D1,LDU,LDW,ii,jj,Nmax)

    double precision U(*),W(*),D2(*),D1(*)
    integer LDU,LDW,ii,jj

    double precision zero,junk,Nmax
    integer j,jjL
    parameter(zero = 0.0D0)

    call CauchyMax(D2,D1(1),u,LDU,junk,jjL)
    junk = junk*abs(w(1))
    Nmax = junk
    ii = 1
    jj = jjL

    do j = 2, LDW   ! ii: col, jj: row
       call CauchyMax(D2,D1(j),u,LDU,junk,jjL)
       junk = junk*abs(w(j))
       if(junk .gt. Nmax) then
          Nmax = junk
          ii = j
          jj = jjL
       end if
    end do

  end subroutine searchMax2

!!!!!!
  subroutine CauchyMax(D,F,u,N,junk,jj)

    double precision D(*),u(*),junk,F
    integer N,jj
    integer temp(1)

!  .. Intrinsic Functions ..
    intrinsic    Maxloc, ABS,MAXVAL

    double precision, allocatable :: LK(:) ! change it to LK(N)

    allocate( LK(N) )

    Lk = u(1:N)/ ( D(1:N)-F )
    Lk = abs(Lk)
    junk = maxval(Lk(1:N))
    temp = maxloc(Lk(1:N))
    jj = temp(1)

    deallocate(LK)

  end subroutine CauchyMax

!!!!!!
  subroutine CauchyPivtG( D,F,u,v,PL,PU,k,A,M,N )
!
    integer k,M,N
    double precision D(*),F(*),u(*),v(*),A(*)
    integer PL(*), PU(*)

!  .. Intrinsic Functions ..
    intrinsic    ABS

!   local parameters
    integer jjL,jjU,flg,jj
    double precision junkL,junkU, Piv, pivot,zero,one
    parameter    (zero = 0.0D0, one=1.0D0)

    call CauchyMax(D(k),F(k),u(k),M-k+1,junkL,jjL)
    call CauchyMax(F(k),D(k),v(k),N-k+1,junkU,jjU)
    junkL = junkL * abs(v(k))
    junkU = junkU * abs(u(k))
    Piv   = abs(u(k)*v(k)/ (D(k)-F(k)) )

    if (junkL .le. Piv .and. junkU .le. Piv ) then
       return
    end if

    pivot = zero
    flg = 0
    if(junkL > junkU) flg = 1

    do while (1 < 2)
       pivot = pivot +1
       if (flg == 1) then
          jj = jjL
          call dswap(D,k,jj+k-1)
          call dswap(u,k,jj+k-1)
          call dswap(A,k,jj+k-1)
          call iswap(PL,k,jj+k-1)
          call CauchyMax( F(k),D(k),v(k),N-k+1,junkU,jjU ) ! N-k+1
          if(jjU == 1) return

          flg = 0
          continue
       end if
       jj = jjU
       call dswap(F,k,jj+k-1)
       call dswap(v,k,jj+k-1)
       call iswap(PU,k,jj+k-1)
       call CauchyMax( D(k),F(k),u(k),M-k+1,junkL,jjL )
       if(jjL == 1) return

       flg = 1
    end do

  end subroutine CauchyPivtG

!!!!!!!!
    subroutine CauchyPivtG_CP( D,F,u,v,PL,PU,k,A,M,N )
!
!   This routine chooses the largest entry for the k-th Schur Complement.
!   It uses complete pivoting. If this routine could be faster than Matrix-Matrix
!   multiplication, we can try to modify the codes for bidiagonal SVD too.
!
    integer k,M,N
    double precision D(*),F(*),u(*),v(*),A(*)
    integer PL(*), PU(*)

!  .. Intrinsic Functions ..
    intrinsic    ABS
! ..
! .. local parameters
    integer ii,jj,j
    double precision junk,zero,one
    parameter    (zero = 0.0D0, one=1.0D0)
!   ii : the row index; jj : the column index;
! ..
! .. local array ..
    integer iiU(N-k+1)    ! the row index for each column
    double precision junkU(N-k+1)  ! the largest entry for each column
!
!$OMP PARALLEL DO
    DO j =1, N-k+1
       call CauchyMax(D(k),F(k+j-1),u(k),M-k+1,junkU(j),iiU(j) )
       junkU(j) = abs( junkU(j)*V(k+j-1) )
    END DO
!$OMP END PARALLEL DO

    junk = junkU(1)
    ii = iiU(1)
    jj = 1
    DO j = 2, N-k+1
       IF(junkU(j) > junk ) THEN
          junk = junkU(j)
          ii = iiU(j)
          jj = j
       END IF
    END DO

    call dswap(D,k,ii+k-1)
    call dswap(u,k,ii+k-1)
    call dswap(A,k,ii+k-1)
    call iswap(PL,k,ii+k-1)

    call dswap(F,k,jj+k-1)
    call dswap(v,k,jj+k-1)
    call iswap(PU,k,jj+k-1)

  end subroutine CauchyPivtG_CP

!!!!!!
  subroutine CauchyPivtG_VP( D,F,u,v,PL,PU,k,A,M,N )
!
    integer k,M,N
    double precision D(*),F(*),u(*),v(*),A(*)
    integer PL(*), PU(*)

!  .. Intrinsic Functions ..
    intrinsic    ABS, MIN

!   local parameters
    integer flg,jj,MN,ii,j,ierr
    double precision junk, zero,one
    parameter    (zero = 0.0D0, one=1.0D0)
!
!   local array
    integer, allocatable :: iiU(:)
    double precision, allocatable :: junkU(:)

    flg = 0
    MN = MIN(10, N-K+1)
    allocate( iiU(MN), junkU(MN), stat=ierr )
    if(ierr .ne. 0 ) then
       write(*,*) 'Allocate failed in PivtG_VP'
    end if

    do while ( flg .ne. 1 )
       call CauchyPivtG(D,F,u,v,PL,PU,k,A,M,N)

       ! maximun of first 10 columns
       DO j =1, MN
          call CauchyMax( D(k),F(k+j-1),u(k),M-k+1,junkU(j),iiU(j) )
          junkU(j) = abs( junkU(j)*V(k+j-1) )
       END DO

       junk = junkU(1)
       ii = iiU(1)
       jj = 1
       DO j = 2, MN
          IF(junkU(j) > junk ) THEN
             junk = junkU(j)
             ii = iiU(j)
             jj = j
          END IF
       END DO

       if( jj .eq. 1 ) then
          flg = 1
          return
       else
          call dswap(D,k,ii+k-1)
          call dswap(u,k,ii+k-1)
          call dswap(A,k,ii+k-1)
          call iswap(PL,k,ii+k-1)

          call dswap(F,k,jj+k-1)
          call dswap(v,k,jj+k-1)
          call iswap(PU,k,jj+k-1)
       end if

    end do

    deallocate(iiU, junkU)
  end subroutine CauchyPivtG_VP

!!!!!!
  subroutine dswap(D,mk,nk)
! double 1D array

    double precision D(*), tt
    integer mk,nk
    tt = D(mk)
    D(mk) = D(nk)
    D(nk) = tt

  end subroutine dswap

  subroutine iswap(D,mk,nk)
! integer 1D array

    integer D(*), tt, mk, nk
    tt = D(mk)
    D(mk) = D(nk)
    D(nk) = tt

  end subroutine iswap

!!!!!!
  subroutine compress_cauchy(rowcol,D,F,U,V,tol,M,N,H,Rk)
!
! Scalar parameters
    integer M,N,Rk
    double precision tol
    character(len=1) rowcol
! Array parameters
    double precision D(*),F(*),U(*),V(*),H(*)
!
! Purpose
! ========
! This routine computes a low rank approximation of some off-diagonal blocks of a
! Cauchy-like matrix.
!
! =========
! Written by S.G. Li, on Sept. 1st, 2019
!
! Note that this version only works for 'Row' compression case.
! =========

! local scalars
    integer mn,i,pnh
    logical  CR

! local arrays
    double precision zero, one
    parameter (zero = 0.0D0, one = 1.0D0)
    integer, allocatable :: PL(:), PU(:)
    double precision, allocatable :: Q(:,:),Z(:),W(:),LD(:),LF(:),LU(:),LV(:)

!  .. Intrinsic Functions ..
    intrinsic    Min

!  .. External Functions
      LOGICAL     LSAME
      EXTERNAL    LSAME

      pnh = 1
      mn = min(M,N)
      allocate(Z(M),W(mn),LD(M),LF(N),LU(M),LV(N),PL(M),PU(N) )
      LD(1:M) = D(1:M)
      LU(1:M) = U(1:M)
      LF(1:N) = F(1:N)
      LV(1:N) = V(1:N)

      Z = zero
      W = zero
      CR = lsame(rowcol,'r')
      if( CR )  then
         call rrluCauchy(LD,LF,LU,LV,tol,Z,W,PL,PU,M,N,Rk)
         call invp(PL,M)       ! PL --> InvPL
         allocate( Q(M,Rk) )
         if (Rk .lt. M) then
            call Cauchylike2(Q,LD(Rk+1),LD,Z(Rk+1),W,M,Rk)
            Q(1:M,1:Rk) = Q(PL,1:Rk)
            call dlacpy('A',M,Rk,Q,M,H(pnh),M)        ! copy tall matrix WH to H
            pnh = pnh+M*Rk
            call Cauchylike( H(pnh),LD,F,LU,V,Rk,N )    ! copy fat matrix QT to H
         else
            ! copy identity matrix to generators
            WRITE(*,*) "Warning: This off-diagonal block is full rank."
         end if

      else ! block col

         call rrluCauchy(F,D,V,U,tol,Z,W,PL,PU,M,N,Rk)
         call invp(PL,M)  ! PL --> InvPL
         allocate( Q(Rk,M) )
         if (Rk .lt. M) then
            call Cauchylike2(Q,F(Rk+1),F,Z(Rk+1),W,Rk,M)
            Q(1:Rk,1:M) = Q(1:Rk,PL)
            call dlacpy('A',Rk,M,Q,Rk,H(pnh),Rk)      ! copy V to H
            !call hssexpmm2('V',PH,pnh,Rk,M,nodi,info) ! copy Q to generators
         else
            ! copy identity matrix to generators
            Q(1:M,1:M) = zero
            do i = 1,M
               Q(i,i) = one
            end do
            Q(1:M,1:M) = Q(1:M,PL)
            call dlacpy('A',M,M,Q,M,H(pnh),M)         ! copy V to H
            !call hssexpmm2('V',PH,pnh,Rk,M,nodi,info) ! copy Q to generators
         end if

      end if ! compr type

      deallocate(Z,W,Q,LD,LF,LU,LV,PL,PU)

    end subroutine compress_cauchy

!!!!!!
    subroutine invp(P,N)
      integer N,P(*)
      integer i
      integer :: IP(N)

      IP = 0

      do i = 1, N
         IP( P(i) ) = i
      end do
      P(1:N) = IP(1:N)

    end subroutine invp

subroutine Cauchylike(A,D,F,U,V,M,N)
!   A(i,j) = U(i)V(j) / (D(i)-F(j))

    integer M,N,i,j
    double precision A(*), D(*),F(*),u(*),v(*)

    do j = 1,N
       do i = 1,M
          A((j-1)*M+i) = U(i)*V(j) / (D(i)-F(j))
       end do
    end do

  end subroutine Cauchylike

!!!!!!
  subroutine Cauchylike2(A,D,F,u,v,M,N)
! A(i,j) = u(i)v(j) / (D(i)-F(j))

    integer M,N,i,j,Rk
    double precision A(M,*), D(*),F(*),u(*),v(*)
    double precision zero, one
    parameter(zero = 0.0D0, one = 1.0D0)

    Rk = M-N
    A(1:M,1:N) = zero

    if(Rk .eq. 0) then
       do j = 1, M
          A(j,j) = one
       end do
    else
       if(Rk .gt. 0) then
          do j = 1,N
             A(j,j) = one
             do i = N+1,M
                A(i,j) = u(i-N)*v(j) / (D(i-N)-F(j))
             end do
          end do
       else
          do j = 1, M
             A(j,j) = one
          end do
          do j = M+1, N
             do i = 1, M
                A(i,j) = v(i)*u(j-M) / (D(j-M)-F(i))
             end do
          end do
       end if
    end if

  end subroutine Cauchylike2
!
!!!!!!
       function comp_time()

         real(8) :: comp_time
         integer :: time0(8)
         call date_and_time(values=time0)
         comp_time = time0(5) * 3600 + time0(6)*60 + time0(7) +0.001 * time0(8)

       end function comp_time
!

end module Cauchylowrank

max_stars_repo_name: shengguolsg/PSMMA

max_issues_repo_name: shengguolsg/PSMMA

max_forks_repo_name: shengguolsg/PSMMA