text:
subroutine hbuf_conditionals_init(count,trcount)
  use vars, only: ncondavg, condavgname, condavglongname, &
       dowtg_blossey_etal_JAMES2009
  use rad, only: do_output_clearsky_heating_profiles
  implicit none

  ! Initialize the list of UW statistics variables written in statistics.f90
  integer count,trcount, n

  if(do_output_clearsky_heating_profiles) then
    call add_to_namelist(count,trcount,'RADQRCLW', &
         'Clearsky longwave heating rate','K/d',0)
    call add_to_namelist(count,trcount,'RADQRCSW', &
         'Clearsky shortwave heating rate','K/d',0)
  end if

  if(dowtg_blossey_etal_JAMES2009) then
    call add_to_namelist(count,trcount,'WWTG', &
         'Large-scale W induced by weak temperature gradient approx','m/s',0)
  end if

  if(dowtg_blossey_etal_JAMES2009) then
    call add_to_namelist(count,trcount,'WOBSREF', &
         'Reference Large-scale W Before Modifications by WTG/Scaling','m/s',0)
  end if

  !bloss: setup to add an arbitrary number of conditional statistics
  do n = 1,ncondavg

     !bloss: add all of the conditional statistics here, so that they don't
     !  have to be added to the lst file
     call add_to_namelist(count,trcount,TRIM(condavgname(n)), &
          TRIM(condavglongname(n))//' Fraction',' ',0)
     call add_to_namelist(count,trcount,'W'//TRIM(condavgname(n)), &
          'Mean W in '//TRIM(condavglongname(n)),'m/s',n)
     call add_to_namelist(count,trcount,'U'//TRIM(condavgname(n)), &
          'Mean U in '//TRIM(condavglongname(n)),'m/s',n)
     call add_to_namelist(count,trcount,'V'//TRIM(condavgname(n)), &
          'Mean V in '//TRIM(condavglongname(n)),'m/s',n)
     call add_to_namelist(count,trcount,'MSE'//TRIM(condavgname(n)), &
          'Mean moist static energy in '//TRIM(condavglongname(n)),'K',n)
     call add_to_namelist(count,trcount,'DSE'//TRIM(condavgname(n)), &
          'Mean dry static energy in '//TRIM(condavglongname(n)),'K',n)
     call add_to_namelist(count,trcount,'TL'//TRIM(condavgname(n)), &
          'Mean liquid-ice static energy in '//TRIM(condavglongname(n)),'K',n)
     call add_to_namelist(count,trcount,'TA'//TRIM(condavgname(n)), &
          'Mean TABS in '//TRIM(condavglongname(n)),'K',n)
     call add_to_namelist(count,trcount,'TV'//TRIM(condavgname(n)), &
          'Mean THETAV in '//TRIM(condavglongname(n)),'K',n)
     call add_to_namelist(count,trcount,'TV'//TRIM(condavgname(n))//'A', &
          'Mean THETAV anomaly in '//TRIM(condavglongname(n)),'K',n)
     call add_to_namelist(count,trcount,'QT'//TRIM(condavgname(n)), &
          'Mean QT in '//TRIM(condavglongname(n)),'g/kg',n)
     call add_to_namelist(count,trcount,'QN'//TRIM(condavgname(n)), &
          'Mean QN in '//TRIM(condavglongname(n)),'g/kg',n)
     !bloss: these conditional averages are now computed inside the microphysics
     !         routines.
     !bloss        call add_to_namelist(count,trcount,'QC'//TRIM(condavgname(n)), &
     !bloss             'Mean QC in '//TRIM(condavglongname(n)),'g/kg',n)
     !bloss        call add_to_namelist(count,trcount,'QI'//TRIM(condavgname(n)), &
     !bloss             'Mean QI in '//TRIM(condavglongname(n)),'g/kg',n)
     call add_to_namelist(count,trcount,'QP'//TRIM(condavgname(n)), &
          'Mean QP in '//TRIM(condavglongname(n)),'g/kg',n)
     call add_to_namelist(count,trcount,'S'//TRIM(condavgname(n)), &
          'Mean scalar in '//TRIM(condavglongname(n)),'K',n)
     call add_to_namelist(count,trcount,'W'//TRIM(condavgname(n))//'A', &
          'W in '//TRIM(condavglongname(n))//' averaged over the whole domain','m/s',0)
     call add_to_namelist(count,trcount,'TLW'//TRIM(condavgname(n)), &
          'TLW in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'Km/s',0)
     call add_to_namelist(count,trcount,'TVW'//TRIM(condavgname(n)), &
          'TVW in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'Km/s',0)
     call add_to_namelist(count,trcount,'SW'//TRIM(condavgname(n)), &
          'SW in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'Km/s',0)
     call add_to_namelist(count,trcount,'QTW'//TRIM(condavgname(n)), &
          'QTW in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'g/kg m/s',0)
     call add_to_namelist(count,trcount,'QCW'//TRIM(condavgname(n)), &
          'QCW in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'g/kg m/s',0)
     call add_to_namelist(count,trcount,'QIW'//TRIM(condavgname(n)), &
          'QIW in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'g/kg m/s',0)

     !bloss: frozen moist static energy statistics
     call add_to_namelist(count,trcount,'HF'//TRIM(condavgname(n)), &
          'Mean Frozen MSE in '//TRIM(condavglongname(n)),'K',n)
     call add_to_namelist(count,trcount,'HF'//TRIM(condavgname(n))//'A', &
          'Mean Frozen MSE anomaly in '//TRIM(condavglongname(n)),'K',n)

     !bloss: velocity anomalies
     call add_to_namelist(count,trcount,'U'//TRIM(condavgname(n))//'A', &
          'Mean U anomaly in '//TRIM(condavglongname(n)),'m/s',n)
     call add_to_namelist(count,trcount,'V'//TRIM(condavgname(n))//'A', &
          'Mean V anomaly in '//TRIM(condavglongname(n)),'m/s',n)

     !bloss: pressure gradients
     call add_to_namelist(count,trcount,'UPGF'//TRIM(condavgname(n)), &
          'Zonal pressure gradient in '//TRIM(condavglongname(n)),'m/s2',n)
     call add_to_namelist(count,trcount,'VPGF'//TRIM(condavgname(n)), &
          'Meridional pressure gradient in '//TRIM(condavglongname(n)),'m/s2',n)
     call add_to_namelist(count,trcount,'WPGF'//TRIM(condavgname(n)), &
          'Vertical pressure gradient in '//TRIM(condavglongname(n)),'m/s2',n)

     !bloss: momentum statistics
     call add_to_namelist(count,trcount,'UW'//TRIM(condavgname(n)), &
          'UW in '//TRIM(condavglongname(n)),'m2/s2',n)
     call add_to_namelist(count,trcount,'VW'//TRIM(condavgname(n)), &
          'VW in '//TRIM(condavglongname(n)),'m2/s2',n)
     call add_to_namelist(count,trcount,'UWSB'//TRIM(condavgname(n)), &
          'Subgrid UW in '//TRIM(condavglongname(n)),'m2/s2',n)
     call add_to_namelist(count,trcount,'VWSB'//TRIM(condavgname(n)), &
          'Subgrid VW in '//TRIM(condavglongname(n)),'m2/s2',n)

     !bloss: UW-added mass flux weighted statistics
     call add_to_namelist(count,trcount,'MF'//TRIM(condavgname(n)), &
          'Mass flux in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'kg/m2/s',0)
     call add_to_namelist(count,trcount,'MFH'//TRIM(condavgname(n)), &
          'RHO*W*HF in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'K kg/m2/s',0)
     call add_to_namelist(count,trcount,'MFH'//TRIM(condavgname(n))//'A', &
          'RHO*W*HF anomaly in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'K kg/m2/s',0)
     call add_to_namelist(count,trcount,'MFTL'//TRIM(condavgname(n)), &
          'RHO*W*TL in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'K kg/m2/s',0)
     call add_to_namelist(count,trcount,'MFTL'//TRIM(condavgname(n))//'A', &
          'RHO*W*TL anomaly in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'K kg/m2/s',0)
     call add_to_namelist(count,trcount,'MFTV'//TRIM(condavgname(n)), &
          'RHO*W*TV in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'K kg/m2/s',0)
     call add_to_namelist(count,trcount,'MFTV'//TRIM(condavgname(n))//'A', &
          'RHO*W*TV anomaly in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'K kg/m2/s',0)
     call add_to_namelist(count,trcount,'MFQT'//TRIM(condavgname(n)), &
          'RHO*W*QT in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'g/m2/s',0)
     call add_to_namelist(count,trcount,'MFQT'//TRIM(condavgname(n))//'A', &
          'RHO*W*QT anomaly in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'g/m2/s',0)
     call add_to_namelist(count,trcount,'RUW'//TRIM(condavgname(n)), &
          'RHOUW in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'kg/m/s2',0)
     call add_to_namelist(count,trcount,'RVW'//TRIM(condavgname(n)), &
          'RHOVW in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'kg/m/s2',0)
     call add_to_namelist(count,trcount,'RWW'//TRIM(condavgname(n)), &
          'RHOWW in '//TRIM(condavglongname(n))//' averaged over the whole domain', &
          'kg/m/s2',0)
  end do ! n = 1,ncondavg

end

subroutine add_to_namelist(count,trcount,varname,varlongname,varunits,varavg)
  use hbuffer, only: namelist,deflist,unitlist,status,average_type
  implicit none

  ! add variable to namelist
  integer count, trcount, ntr, n, varstatus, varavg
  character(*) varname
  character(*) varlongname
  character(*) varunits

  count = count + 1
  trcount = trcount + 1
  namelist(count) = trim(varname)
  deflist(count) = trim(varlongname)
  unitlist(count) = trim(varunits)
  status(count) = 1
  average_type(count) = varavg

end subroutine add_to_namelist

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


program backwords
use M_process,  only : process_readall
implicit none
integer                      :: ierr
integer                      :: i
character(len=:),allocatable :: string
character(len=4096)          :: command
   write(*,'(a)',advance='no')'Enter command to read output from:'
   read(*,'(a)')command
   ! NOTE: the file is being read into memory. Doing this with a
   !       large file would use an excessive amount of memory.
   string=process_readall(command,delim=NEW_LINE("A"),ierr=ierr)
   ! parse on newline character and process the data
   ! using routines like split(3f) from M_strings(3fm) or in
   ! this case, just echo it back reversed
   write(*,'(a)')'BACKWARDS'
   write(*,'(*(a))',advance='no')(string(i:i),i=len(string),1,-1)
   write(*,'(a)')
   write(*,*)'FORWARDS'
   write(*,'(a)',advance='no')string
end program backwords


text:
      MODULE kab_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  10:47:22  03/09/06
      SUBROUTINE kab (IA, JA, PK, W, F)
      USE vast_kind_param,ONLY: DOUBLE
      INTEGER, INTENT(IN) :: IA
      INTEGER, INTENT(IN) :: JA
      REAL(DOUBLE), DIMENSION(*), INTENT(IN) :: PK
      REAL(DOUBLE), DIMENSION(*), INTENT(IN) :: W
      REAL(DOUBLE), DIMENSION(*), INTENT(INOUT) :: F
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
c-----------------------------------------------------------------------
c  activate making inicon table by the follwing commands in optns file:
c-----------------------------------------------------------------------
c        set nevent 1
c        set bminim 3      (for example)
c        set bmaxim 5      (for example)
c        set ninicon 1000  (for example)
c        set icocore 1  !or 2
c        make icotable
c        switch  splitting on
c        switch  fusion on
c        switch clusterdecay off
c-----------------------------------------------------------------------
c  inicon table means storage of IcoE, IcoV, IcoF
c     IcoE ... energy density in the comoving frame
c     IcoV ... flow velocity in hyperbolic coordinates
c               v=(vx,vy,veta=vz/tau)
c                with the velocity in the frame moving with rapidity eta
c     IcoV ... net flavor density in the comoving frame
c  the indices of these fields ix,iy,iz also refer to hyperbolic coordinates
c     (x,y,eta) at given tau
c     the corresponding grid is defined in epos.incico
c-----------------------------------------------------------------------

c-----------------------------------------------------------------------
      subroutine IniCon(nin)
c-----------------------------------------------------------------------
      include 'epos.inc'
      include 'epos.incico'
      double precision seedf
      character*80 fn
      logical table

      ii=index(fnhi(1:nfnhi),".histo")-1
      fn=fnhi(1:ii)//".ico"
      inquire(file=fn(1:ii+4),exist=table)

      if(icotabr.eq.1)then  !read from file

        if(.not.table)then
          write(ifmt,'(//10x,3a//)')'STOP: file '
     *             ,fn(1:ii+4),' not found !!!'
          stop
        endif
        write(ifmt,'(3a)')'read from ',fn(1:ii+4),' ...'
        open(97,file=fn,status='old')
        read(97,*) iversn
        read(97,*) laprojx,maprojx,latargx,matargx
        read(97,*) engyx
        read(97,*) bminimx,bmaximx
        read(97,*) tauicox
        read(97,*) iabs_ninicon_dum
        read(97,*) nxicox,nyicox,nzicox
        read(97,*)xminicox,xmaxicox,yminicox,ymaxicox,zminicox,zmaxicox
        if(laprojx.ne.laproj.or.maprojx.ne.maproj
     * .or.latargx.ne.latarg.or.matargx.ne.matarg)stop'(1)'
        if(engyx.ne.engy)    stop'variables dont match (2)    '
        if(bminimx.ne.bminim.or.bmaximx.ne.bmaxim)
     *                       stop'variables dont match (3)    '
        if(tauicox.ne.tauico)stop'variables dont match (4)    '
        if(nxicox.ne.nxico.or.nyicox.ne.nyico
     * .or.nzicox.ne.nzico)stop'(5)'
        if( xminicox.ne.xminico.or.xmaxicox.ne.xmaxico
     * .or.yminicox.ne.yminico.or.ymaxicox.ne.ymaxico
     * .or.zminicox.ne.zminico.or.zmaxicox.ne.zmaxico)stop'(6)'
        read(97,*) IcoE,IcoV,IcoF
        close(97)

      elseif(icotabr.eq.0)then  ! calculate

        !if(icotabm.eq.1.and.table)then
          !write(ifmt,'(//10x,3a//)')'STOP: file '
          !*             ,fn(1:ii+4),' already exists !!!'
          !stop
        !endif

        if(nin.eq.1)then
          write(ifmt,'(2a,i7,a)')' calculate inicon table ',
     &    'based on',iabs(ninicon),' ico-events ...'
          call IcoStr(1)
        endif

        call ranfgt(seedf)
        if(nin.le.iabs(ninicon).and.mod(nin,modsho).eq.0)
     &        write(ifmt,'(a,i7,a,f6.2,a,d27.16)')' ico-event ',nin
     &                   ,'   bimevt:',bimevt,'   seedf:',seedf
        if(nin.le.iabs(ninicon).and.jwseed.eq.1)then
          open(unit=1,file=fnch(1:nfnch-5)//'see',status='unknown')
          write(1,'(a,i7,a,f6.2,a,d27.16)')' ico-event ',nin
     &                   ,'   bimevt:',bimevt,'   seedf:',seedf
         close(1)
        endif
        !if(nin.le.iabs(ninicon).and.mod(nin,modsho).eq.0)
        !&         print*,'+++++ time: ',timefin-timeini
        seedc=seedf
        call IcoStr(2)

        if(nin.eq.iabs(ninicon))then
          write(ifmt,*)'normalize and diagonalize engy-mom tensor'
           call IcoStr(3)

           if(icotabm.eq.1)then  ! make table
            write(ifmt,'(2a)')' write inicon table into ',fn(1:ii+4)
            open(97,file=fn,status='unknown')
            write(97,*) iversn
            write(97,*) laproj,maproj,latarg,matarg
            write(97,*) engy
            write(97,*) bminim,bmaxim
            write(97,*) tauico
            write(97,*) iabs(ninicon)
            write(97,*) nxico,nyico,nzico
            write(97,*) xminico,xmaxico,yminico,ymaxico,zminico,zmaxico
            write(97,*) IcoE,IcoV,IcoF
            close(97)
          endif
         endif

      else

        stop'IniCon: wrong choice'

      endif

      end

c-----------------------------------------------------------------------
      subroutine IcoStr(iflag)
c-----------------------------------------------------------------------
c     energy momentum tensor and flavor current
c     and corresponding contractions from string method
c
c     iflag = 1 initialization
c           = 2 sum up density
c           = 3 normalization
c
c     output: common blocks /Ico3/  /Ico4/  /Ico5/
c-----------------------------------------------------------------------
      include 'epos.inc'
      include 'epos.incico'
      double precision IcoT(4,4,nxico,nyico,nzico)
     *                ,IcoC(3,4,nxico,nyico,nzico)
      common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
     *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
      common /Ico10/ntot,ish0,irs
      double precision p(4),v(4),u(4),eps,T(4,4),gamma,rap,rapx
      integer ic(2),jc(nflav,2)
      common/cranphi/ranphi
      logical ok
      goto (1,10,100),iflag
      return

c......................................................................
c                             initialization
c......................................................................

 1    do ix=1,nxico
        do iy=1,nyico
          do iz=1,nzico
            IcoE(ix,iy,iz)=0.
            do i=1,3
              IcoV(i,ix,iy,iz)=0.
            enddo
            do i=1,3
              IcoF(i,ix,iy,iz)=0.
            enddo
            do i=1,4
            do j=1,4
              IcoT(i,j,ix,iy,iz)=0.
            enddo
            enddo
            do i=1,3
            do j=1,4
              IcoC(i,j,ix,iy,iz)=0.
            enddo
            enddo
          enddo
        enddo
      enddo
      ntot=0
      return

c......................................................................
c                        fill arrays
c......................................................................

 10   ntot=ntot+1
      dxico=(xmaxico-xminico)/nxico
      dyico=(ymaxico-yminico)/nyico
      dzico=(zmaxico-zminico )/nzico
      phinll=phievt+ranphi

      do j=1,nptl
       ok=.true.
       if(icocore.eq.2
     . .and.(ityptl(j).lt.20.or.ityptl(j).ge.40))ok=.false.
       if(ok
     . .and.dezptl(j).lt.1e3.and.j.le.maxfra.and.istptl(j).eq.2)then
        aa=cos(phinll)
        bb=sin(phinll)
        x=xptl(j)*aa+yptl(j)*bb
        cc=-sin(phinll)
        dd=cos(phinll)
        y=xptl(j)*cc+yptl(j)*dd
        z=dezptl(j)
        ix=1+(x-xminico)/dxico
        iy=1+(y-yminico)/dyico
        iz=1+(z-zminico)/dzico
        if(ix.ge.1.and.ix.le.nxico.and.
     .  iy.ge.1.and.iy.le.nyico.and.
     .  iz.ge.1.and.iz.le.nzico)then
          !~~~~~~~~~~~~
          ! T^i^k = \int d3p p^i p^k /E f(\vec x,\vec p)
          ! N^k =   \int d3p p^k /E f_net(\vec x,\vec p)
          !~~~~~~~~~~~~
          rapx=zminico+(iz-0.5)*dzico
          amt2=pptl(5,j)**2+pptl(1,j)**2+pptl(2,j)**2
          pp=pptl(4,j)+abs(pptl(3,j))
          if(amt2.gt.0..and.pp.gt.0.)then
           amt=sqrt(amt2)
           rap=sign(1.,pptl(3,j))*alog(pp/amt)
           p(1)=pptl(1,j)
           p(2)=pptl(2,j)
           p(3)=amt*sinh(rap-rapx)
           p(4)=amt*cosh(rap-rapx)
           do i=1,4
            do k=1,4
              IcoT(i,k,ix,iy,iz)=IcoT(i,k,ix,iy,iz)+p(i)*p(k)/p(4)
            enddo
           enddo
           id=idptl(j)
           ida=iabs(id/10)
           ids=id/iabs(id)
           if(ida.ne.111.and.ida.ne.222.and.ida.ne.333)id=id/10*10
           if(ida.eq.111.or. ida.eq.222.or. ida.eq.333)id=id/10*10+ids
           if(ida.eq.213)id=1230*ids
           ic(1)=idtrai(1,id,1)
           ic(2)=idtrai(2,id,1)
           call iddeco(ic,jc)
           do i=1,3
            fi=jc(i,1)-jc(i,2)
            do k=1,4
             IcoC(i,k,ix,iy,iz)=IcoC(i,k,ix,iy,iz)+p(k)/p(4)*fi
            enddo
           enddo
          endif
          !~~~~~~~~~~~~
        endif
       endif
      enddo

      return

c............................................................................
c                 normalization and diagonalization
c............................................................................

 100  continue

      !~~~normalize

      vol= (xmaxico-xminico)/nxico
     .    *(ymaxico-yminico)/nyico
     .    *(zmaxico-zminico)/nzico  *tauico
      fac=1./ntot/vol

      do ix=1,nxico
        do iy=1,nyico
          do iz=1,nzico
            do i=1,4
             do k=1,4
              IcoT(i,k,ix,iy,iz)=IcoT(i,k,ix,iy,iz)*fac
             enddo
            enddo
            do i=1,3
             do k=1,4
              IcoC(i,k,ix,iy,iz)=IcoC(i,k,ix,iy,iz)*fac
             enddo
            enddo
          enddo
        enddo
      enddo

      !~~~diagonalize T

      do ix=1,nxico
        do iy=1,nyico
          do iz=1,nzico
            do i=1,4
            do k=1,4
             T(i,k)=IcoT(i,k,ix,iy,iz)
            enddo
            enddo
            call DiagTmunu(T,u,v,gamma,eps,nonz,ix,iy,iz)
            if(nonz.eq.1)then
              IcoE(ix,iy,iz)=eps
              IcoV(1,ix,iy,iz)=v(1)
              IcoV(2,ix,iy,iz)=v(2)
               !rapx=zminico+(iz-0.5)*dzico
               !rap=rapx+log(gamma*(1+v(3)))
               !IcoV(3,ix,iy,iz)=tanh(rap)
              IcoV(3,ix,iy,iz)=v(3) / tauico
              do i=1,3
               IcoF(i,ix,iy,iz)=IcoC(i,4,ix,iy,iz)*u(4)
               do  k=1,3
                IcoF(i,ix,iy,iz)=IcoF(i,ix,iy,iz)
     .                  -IcoC(i,k,ix,iy,iz)*u(k)
               enddo
              enddo
            endif
          enddo
        enddo
      enddo

      return
      end

c------------------------------------------------------------------------------------
      subroutine DiagTmunu(Tmunu,u,v,gamma,eps,nonz,ix,iy,iz)
c------------------------------------------------------------------------------------
      ! solve lambda * T * lambda = diag(eps,P,P,P), for symmetric T
      !
      !  lambda =  g      g*vx        g*vy        g*vz          =symmetric!
      !            g*vx   a*vx*vx+1   a*vy*vx     a*vz*vx
      !            g*vy   a*vx*vy     a*vy*vy+1   a*vz*vy
      !            g*vz   a*vx*vz     a*vy*vz     a*vz*vz+1
      ! with g=gamma, a=g*g/(g+1)
      !
      ! so T*lambda(v)=lambda(-v)*diag(eps,P,P,P)
      ! first column: four equations
      !                    eps=T00+T0k*vk
      !                -eps*vi=Ti0+Tik*vk
      ! solved iteratively to get eps, vi
      ! returns u,v,eps if nonz=1 (otherwise zero tensor)
      !  (by T.Kodama)
      !  (K. Werner: modifs)

      include 'epos.incico'
      double precision Tmunu(4,4), u(4),eps,gamma,g,a, Lor(4,4),w(4),sum
      double precision v(4),vx(3),tt(4,4),err,sg(4)

      sg(4)=1d0
      v(4)=1d0
      do i=1,3
       sg(i)=-1d0
      enddo
      sum=0d0
      do i=1,4
       do k=1,4
       sum=sum+dabs(Tmunu(i,k))
       end do
      end do
      nonz=0
      if(sum.eq.0.0d0)return
      nonz=1

      do k=1,3
       v(k)=0.
      end do
      eps=0

      do lrep=1,100
       epsx=eps
       do i=1,3
        vx(i)=v(i)
       end do
       eps=Tmunu(4,4)
       do k=1,3
        eps=eps-Tmunu(4,k)*vx(k)
       end do
       if(eps.le.0d0)then
         print*,'DiagTmunu: sum(abs(Tmunu))=',sum,'   eps=',eps
         print*,'Tmunu(4,nu):',(Tmunu(4,nu),nu=1,4)
         if(eps.gt.-1e-5)then
           nonz=0
           return
         else
           print*,'STOP in DiagTmunu: negative energy density.  '
           stop'(3003200808)'
         endif
       endif
       do i=1,3
        Tv=0d0
        do k=1,3
         Tv=Tv+Tmunu(i,k)*vx(k)
        end do
        v(i)=(Tmunu(i,4)-Tv)/eps
       end do
       if(lrep.gt.60)then
        do i=1,3
         v(i)=0.5d0*(vx(i)+v(i))
        enddo
       endif
       !print*,'Tmunu: ',lrep,abs(eps-epsx),(abs(v(i)-vx(i)),i=1,3)
       err=1d-6
       if(lrep.gt.50)err=1d-5
       if(lrep.gt.89)err=1d-4
       if(abs(eps-epsx).lt.err.and.abs(v(1)-vx(1)).lt.err
     . .and.abs(v(2)-vx(2)).lt.err.and.abs(v(3)-vx(3)).lt.err)goto1
        do i=1,4
          w(i)=0
          do k=1,4
          w(i)=w(i)+Tmunu(i,k)*v(k)*sg(k)
          enddo
          w(i)=w(i)-eps*v(i)
        enddo
        if(lrep.gt.95
     ..and.w(1)*w(1)+w(2)*w(2)+w(3)*w(3)+w(4)*w(4).lt.err)goto1
      end do

  1   v2=0.d0
      do i=1,3
       v2=v2+v(i)**2
      end do
      if(v2.ge.1.)stop'DiagTmunu: v2 ge 1.          '
      gamma=1./sqrt(abs(1.-v2))
      u(4)=gamma
      u(1)=v(1)*gamma
      u(2)=v(2)*gamma
      u(3)=v(3)*gamma

      !~~~check
      g=gamma
      a=g*g/(g+1d0)
      Lor(4,4)=g
      do k=1,3
       Lor(4,k)=-g*v(k)
       Lor(k,4)=Lor(4,k)
      enddo
      do i=1,3
      do k=1,3
       Lor(i,k)=a*v(i)*v(k)
      enddo
      enddo
      do k=1,3
       Lor(k,k)=Lor(k,k)+1
      enddo
      do i=1,4
      do k=1,4
        tt(i,k)=0d0
        do m=1,4
        do n=1,4
          tt(i,k)=tt(i,k)+Lor(i,m)*Tmunu(m,n)*Lor(n,k)
        enddo
        enddo
      enddo
      enddo
      err=err*1d2
      if(tt(1,4).gt.err.or.tt(2,4).gt.err.or.tt(3,4).gt.err)then
        print*,'************** DiagTmunu: nonzero T14 or T24 or T34'
     .  ,'  after ',lrep,' iterations'
        print*,'d_eps or d_v(i): ',abs(eps-epsx),(abs(v(i)-vx(i)),i=1,3)
        print*,'Tmunu ( ix=',ix-nxico/2-1,'  iy=',iy-nyico/2-1,' ) :'
     .  ,'  iz=',iz-nzico/2-1
        do i=1,4
        write(*,'(4f10.5,2x,4f10.5)')
     .  (sngl(tmunu(i,k)),k=1,4),(sngl(tt(i,k)),k=1,4)
        enddo
      endif

      end

ext: f

max_stars_repo_name: ckamtsikis/cmssw

max_issues_repo_name: ckamtsikis/cmssw

max_forks_repo_name: ckamtsikis/cmssw


text:
      SUBROUTINE GETWIN(NODE,LUN,IWIN,JWIN)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    GETWIN (docblock incomplete)
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 1994-01-06
C
C ABSTRACT: THIS SUBROUTINE ....
C
C PROGRAM HISTORY LOG:
C 1994-01-06  J. WOOLLEN -- ORIGINAL AUTHOR
C 1998-07-08  J. WOOLLEN -- REPLACED CALL TO CRAY LIBRARY ROUTINE
C                           "ABORT" WITH CALL TO NEW INTERNAL BUFRLIB
C                           ROUTINE "BORT"
C 1999-11-18  J. WOOLLEN -- THE NUMBER OF BUFR FILES WHICH CAN BE
C                           OPENED AT ONE TIME INCREASED FROM 10 TO 32
C                           (NECESSARY IN ORDER TO PROCESS MULTIPLE
C                           BUFR FILES UNDER THE MPI)
C 2002-05-14  J. WOOLLEN -- REMOVED OLD CRAY COMPILER DIRECTIVES
C 2003-11-04  S. BENDER  -- ADDED REMARKS/BUFRLIB ROUTINE
C                           INTERDEPENDENCIES
C 2003-11-04  D. KEYSER  -- MAXJL (MAXIMUM NUMBER OF JUMP/LINK ENTRIES)
C                           INCREASED FROM 15000 TO 16000 (WAS IN
C                           VERIFICATION VERSION); UNIFIED/PORTABLE FOR
C                           WRF; ADDED DOCUMENTATION (INCLUDING
C                           HISTORY) (INCOMPLETE); OUTPUTS MORE
C                           COMPLETE DIAGNOSTIC INFO WHEN ROUTINE
C                           TERMINATES ABNORMALLY
C DART $Id$
C
C USAGE:    CALL GETWIN (NODE, LUN, IWIN, JWIN)
C   INPUT ARGUMENT LIST:
C     NODE     - INTEGER: ....
C     LUN      - INTEGER: I/O STREAM INDEX INTO INTERNAL MEMORY ARRAYS
C     IWIN     - INTEGER: ....
C     JWIN     - INTEGER: ....
C
C   OUTPUT ARGUMENT LIST:
C     IWIN     - INTEGER: ....
C     JWIN     - INTEGER: ....
C
C REMARKS:
C    THIS ROUTINE CALLS:        BORT     INVWIN   LSTRPC
C    THIS ROUTINE IS CALLED BY: CONWIN   UFBEVN   UFBIN3   UFBRW
C                               Normally not called by any application
C                               programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      INCLUDE 'bufrlib.prm'

      COMMON /USRINT/ NVAL(NFILES),INV(MAXJL,NFILES),VAL(MAXJL,NFILES)

      CHARACTER*128 BORT_STR
      REAL*8        VAL

C----------------------------------------------------------------------
C----------------------------------------------------------------------

      IRPC = LSTRPC(NODE,LUN)

      IF(IRPC.EQ.0) THEN
         IWIN = INVWIN(NODE,LUN,JWIN,NVAL(LUN))
         IF(IWIN.EQ.0 .and. jwin.gt.1) GOTO 100
         IWIN = 1
         JWIN = NVAL(LUN)
         GOTO 100
      ELSE
         IWIN = INVWIN(IRPC,LUN,JWIN,NVAL(LUN))
         IF(IWIN.EQ.0) THEN
            GOTO 100
         ELSEIF(VAL(IWIN,LUN).EQ.0.) THEN
            IWIN = 0
            GOTO 100
         ENDIF
      ENDIF

      JWIN = INVWIN(IRPC,LUN,IWIN+1,NVAL(LUN))
      IF(JWIN.EQ.0) GOTO 900

C  EXITS
C  -----

100   RETURN
900   WRITE(BORT_STR,'("BUFRLIB: GETWIN - SEARCHED BETWEEN",I5," AND"'//
     . ',I5,", MISSING BRACKET")') IWIN+1,NVAL(LUN)
      CALL BORT(BORT_STR)
      END

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                          Futility Development Group                          !
!                             All rights reserved.                             !
!                                                                              !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !
! can be found in LICENSE.txt in the head directory of this repository.        !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!> @brief Utility module for defining matrix types.
!>
!> The types of matrices defined in this module include a sparse matrix type,
!> square and rectangle dense matrix types, and a square tridiagonal matrix
!> type. Presently, the dense matrices are all assumed to be stored in a general
!> format as opposed to packed, banded, or upper/lower triangular formats.
!> Interfaces to BLAS routines for matrix-vector multiplication and
!> matrix-matrix multiplication are also added to the global generic interfaces
!> for @ref BLAS2::BLAS_matvec "BLAS_matvec" and BLAS_matmat interfaces.
!>
!> The objects are initialized with a parameter list. For valid reference lists
!> see @ref MatrixTypes::MatrixTypes_Declare_ValidParams
!> "MatrixTypes_Declare_ValidParams".
!>
!> NOTE: Sparse formats are CSR formats that match the 3-Array Variation Format
!> of the Intel MKL. However, filling these requires traversing the matrix
!> in row-major format (i.e., increment i in inner loop). This is more how a C
!> programmer would do it, not us Fortran folk.  Something to be aware of.
!>
!> TIBWSFB: This documentation doesnt match the current interface
!> @par EXAMPLES
!> @code
!> PROGRAM ExampleMatrix
!>   TYPE(SparseMatrixType) :: sparse
!>   CLASS(MatrixType),POINTER :: matrix_p => NULL()
!>   TYPE(ParamType) :: params
!>
!>   CALL params%add("MatrixType->n",3_SIK)
!>   CALL params%add("MatrixType->nnz",6_SIK)
!>
!>   CALL sparse%init(params)
!>   CALL params%clear()
!>   CALL sparse%setShape(1,1,1._SRK)
!>   CALL sparse%setShape(1,3,2._SRK)
!>   CALL sparse%setShape(2,3,3._SRK)
!>   CALL sparse%setShape(3,1,4._SRK)
!>   CALL sparse%setShape(3,2,5._SRK)
!>   CALL sparse%setShape(3,3,6._SRK)
!>   CALL sparse%set(1,1,10._SRK)
!>   CALL sparse%set(1,3,20._SRK)
!>   CALL sparse%set(2,3,30._SRK)
!>   CALL sparse%set(3,1,40._SRK)
!>   CALL sparse%set(3,2,50._SRK)
!>   CALL sparse%set(3,3,60._SRK)
!>   CALL sparse%clear()
!>
!>   ! Create a PETSc matrix using the factory
!>   CALL params%add("MatrixType->n",3_SIK)
!>   CALL params%add("MatrixType->nnz",6_SIK)
!>   CALL params%add("MatrixType->matType",SPARSE)
!>   CALL params%add("MatrixType->engine",VM_PETSC)
!>   CALL MatrixFactory(matrix_p,params)
!>
!>   ! Clean up
!>   CALL matrix_p%clear()
!>   DEALLOCATE(matrix_p)
!>
!> ENDPROGRAM ExampleMatrix
!> @endcode
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
MODULE MatrixTypes
USE IntrType
USE Allocs
USE BLAS2,           ONLY: BLAS2_matvec => BLAS_matvec
USE BLAS3,           ONLY: BLAS3_matmult => BLAS_matmat
USE ParameterLists
USE VectorTypes
USE MatrixTypes_Base
USE MatrixTypes_Native
USE MatrixTypes_Trilinos
USE MatrixTypes_PETSc
USE trilinos_interfaces
IMPLICIT NONE

#ifdef HAVE_MPI
#include <mpif.h>
#endif

PRIVATE
!
! List of public members
PUBLIC :: MatrixType
PUBLIC :: MatrixFactory
PUBLIC :: DistributedMatrixFactory
PUBLIC :: MatrixResemble
PUBLIC :: SquareMatrixType
PUBLIC :: RectMatrixType
PUBLIC :: DistributedMatrixType
! Matrix structure enumerations
PUBLIC :: SPARSE,DENSESQUARE,DENSERECT,TRIDIAG, &
    BANDED,DISTRIBUTED_BANDED,DISTR_BLOCKBANDED
! Matrix-Vector engine enumerations
PUBLIC :: VM_PETSC,VM_TRILINOS,VM_NATIVE
! Parameter list setup/teardown
PUBLIC :: MatrixTypes_Declare_ValidParams
PUBLIC :: MatrixTypes_Clear_ValidParams
! Native implementations
PUBLIC :: DenseSquareMatrixType
PUBLIC :: DenseRectMatrixType
PUBLIC :: TriDiagMatrixType
PUBLIC :: BandedMatrixType
PUBLIC :: DistributedBandedMatrixType
PUBLIC :: DistributedBlockBandedMatrixType
PUBLIC :: SparseMatrixType
INTEGER(SIK),PARAMETER,PUBLIC :: MATVEC_SLOTS=10
! PETSc implementations
#ifdef FUTILITY_HAVE_PETSC
PUBLIC :: PETScMatrixType
#endif
#ifdef FUTILITY_HAVE_Trilinos
PUBLIC :: TrilinosMatrixType
#endif

PUBLIC :: eMatrixType
PUBLIC :: BLAS_matvec
PUBLIC :: BLAS_matmult

!> @brief Adds to the @ref BLAS2::BLAS_matvec "BLAS_matvec" interface so that
!> the matrix types defined in this module are also supported.
INTERFACE BLAS_matvec
  !> @copybrief MatrixTypes::matvec_MatrixType
  !> @copydetails MatrixTypes::matvec_MatrixType
  MODULE PROCEDURE matvec_MatrixType
  !> @copybrief MatrixTypes::matvec_MatrixTypeVectorType
  !> @copydetails MatrixTypes::matvec_MatrixTypeVectorType
  MODULE PROCEDURE matvec_MatrixTypeVectorType
ENDINTERFACE BLAS_matvec

!> @brief Adds to the @ref BLAS3::BLAS_matmult "BLAS_matmat" interface so that
!> the matrix types defined in this module are also supported.
INTERFACE BLAS_matmult
  !> @copybrief MatrixTypes::matmult_MatrixType
  !> @copydetails MatrixTypes::matmult_MatrixType
  MODULE PROCEDURE matmult_MatrixType
ENDINTERFACE BLAS_matmult

!> @brief Adds strsv_all and dtrsv_all routines for sparse matrices to the
!> @ref BLAS2::BLAS_matvec "BLAS_matvec" interface
!  INTERFACE trsv_sparse
  !> @copybrief MatrixTypes::strsv_all_sparse
  !> @copydetails MatrixTypes::strsv_all_sparse
!    MODULE PROCEDURE strsv_all_sparse
  !> @copybrief MatrixTypes::dtrsv_all_sparse
  !> @copydetails MatrixTypes::dtrsv_all_sparse
!    MODULE PROCEDURE dtrsv_all_sparse
!  ENDINTERFACE trsv_sparse

!> Name of module
CHARACTER(LEN=*),PARAMETER :: modName='MATRIXTYPES'

!
!===============================================================================
CONTAINS
!
!-------------------------------------------------------------------------------
!> @brief Abstract factory for all enabled MatrixTypes
!>
!> @param matrix a pointer to the base MatrixType that will store the
!> constructed matrix. Should be NULL
!> @param params a parameter list to use to determine the type of and construct
!> the matrix
!>
!> This is an abstract factory routine for the base MatrixTypes. It uses the
!> matType and engine parameters on the passed parameter list to determine which
!> type of MatrixType to which the passed pointer should be allocated. It then
!> uses the same parameter list to initialize the matrix.
!>
!> The valid matrix engines (if enabled at configure time) are MV_NATIVE,
!> MV_TRILINOS and MV_PETSC.
!> The valid matTypes are those enumerated in the MatrixTypes_Base module.
!>
!> The purpose of such a routine is to hide from the user which specific
!> MatrixTypes are available and reduce the need for SELECT TYPEing the matrices
!> later. This also hides the #ifdef guards needed to safely interact with each
!> of the optional engines.
SUBROUTINE MatrixFactory(matrix, params)
  CHARACTER(LEN=*),PARAMETER :: myName="MatrixFactory"
  CLASS(MatrixType),POINTER,INTENT(INOUT) :: matrix
  CLASS(ParamType),INTENT(IN) :: params
  !
  INTEGER(SIK) :: engine,matType
  CLASS(DistributedMatrixType),POINTER :: dist_p

  IF(ASSOCIATED(matrix)) THEN
    CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
        "Matrix pointer already allocated")
    RETURN
  ENDIF

  dist_p => NULL()

  ! Default to dative, dense matrix
  engine=VM_NATIVE
  matType=DENSESQUARE

  IF(params%has("MatrixType->engine")) THEN
    CALL params%get("MatrixType->engine", engine)
  ENDIF

  IF(params%has("MatrixType->matType")) THEN
    CALL params%get("MatrixType->matType", matType)
  ELSE
    CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
        "no matrix TYPE specified")
  ENDIF

  SELECTCASE(engine)
  CASE(VM_NATIVE)
    SELECTCASE(matType)
    CASE(DENSESQUARE)
      ALLOCATE(DenseSquareMatrixType :: matrix)
    CASE(DENSERECT)
      ALLOCATE(DenseRectMatrixType :: matrix)
    CASE(SPARSE)
      ALLOCATE(SparseMatrixType :: matrix)
    CASE(TRIDIAG)
      ALLOCATE(TriDiagMatrixType :: matrix)
    CASE(BANDED)
      ALLOCATE(BandedMatrixType :: matrix)
    CASE(DISTRIBUTED_BANDED)
      ALLOCATE(DistributedBandedMatrixType :: matrix)
    CASE(DISTR_BLOCKBANDED)
      ALLOCATE(DistributedBlockBandedMatrixType :: matrix)
    CASE DEFAULT
      CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
          "Unrecognized matrix structure requested")
    ENDSELECT
  CASE(VM_TRILINOS)
    CALL DistributedMatrixFactory(dist_p, params)
    matrix => dist_p
    RETURN
  CASE(VM_PETSC)
    CALL DistributedMatrixFactory(dist_p, params)
    matrix => dist_p
    RETURN
  CASE DEFAULT
    CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
        "Unsupported matrix engine requested.")
  ENDSELECT

  CALL matrix%init(params)
ENDSUBROUTINE MatrixFactory
!
!-------------------------------------------------------------------------------
!> @brief Abstract factory for all enabled DistributedMatrixTypes
!>
!> @param matrix a pointer to the DistributedMatrixType that will store the
!> constructed matrix. Should be NULL
!> @param params a parameter list to use to determine the type of and construct
!> the matrix
!>
!> This is similar to the MatrixType factory, but specifically interacts with
!> DistributedMatrixType pointers. This distinction in necessary for client code
!> that knows that it wants one of the DistributedMatrixTypes and therefore uses
!> a pointer to such.
SUBROUTINE DistributedMatrixFactory(matrix, params)
  CHARACTER(LEN=*),PARAMETER :: myName="DistributedMatrixFactory"
  CLASS(DistributedMatrixType),POINTER,INTENT(INOUT) :: matrix
  CLASS(ParamType),INTENT(IN) :: params
  !
  INTEGER(SIK) :: engine,matType

  IF(ASSOCIATED(matrix)) THEN
    CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
        "Matrix pointer already allocated")
    RETURN
  ENDIF

  matType=SPARSE

  IF(params%has("MatrixType->engine")) THEN
    CALL params%get("MatrixType->engine", engine)
  ENDIF

  IF(params%has("MatrixType->matType")) THEN
    CALL params%get("MatrixType->matType", matType)
  ENDIF

  SELECTCASE(engine)
  CASE(VM_TRILINOS)
#ifdef FUTILITY_HAVE_Trilinos
    IF(matType == SPARSE) THEN
      ALLOCATE(TrilinosMatrixType :: matrix)
    ELSE
      CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
          "Trilinos matrix should be sparse")
    ENDIF
#else
    CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
        "Futility was not compiled with Trilinos support. Recompile "// &
      "with Trilinos to create PETSc matrices.")
#endif
  CASE(VM_PETSC)
#ifdef FUTILITY_HAVE_PETSC
    IF(matType == SPARSE .OR. matType == DENSESQUARE) THEN
      ALLOCATE(PETScMatrixType :: matrix)
    ELSE
      CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
          "PETSc matrix should be sparse")
    ENDIF
#else
    CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
        "Futility was not compiled with PETSc support. Recompile with "// &
      "PETSc to create PETSc matrices.")
#endif
  CASE DEFAULT
    CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
        "Unsupported distributed matrix engine requested.")
  ENDSELECT

  CALL matrix%init(params)
ENDSUBROUTINE
!
!-------------------------------------------------------------------------------
!> @brief Abstractly reproduce a new Matrix of the same type passed in
!>
!> @param dest the matrix pointer to allocate
!> @param source the matrix to use when determining the proper type for dest
!> @param params the parameter list to use when constructing the dest matrix
!>
!> This routine uses the matrix passed in as source to allocate a new matrix,
!> dest, of the same type. The passed parameter list is then used to initialize
!> the dest matrix.
!> Such a routine exists to make it easy for client code to create a new matrix
!> that will be compatible with an existing matrix, without needed to know which
!> type the existing matrix is. This cuts down on #ifdef guards and SELECT TYPEs
SUBROUTINE MatrixResemble(dest, source, params)
  CHARACTER(LEN=*),PARAMETER :: myName="MatrixResemble"
  CLASS(MatrixType),POINTER,INTENT(INOUT) :: dest
  CLASS(MatrixType),POINTER,INTENT(IN) :: source
  CLASS(ParamType),INTENT(INOUT) :: params

  IF(.NOT. ASSOCIATED(source)) THEN
    CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
        "Source matrix is not associated")
    RETURN
  ENDIF

  IF(ASSOCIATED(dest)) THEN
    CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
        "Destination pointer is already associated")
    RETURN
  ENDIF

  SELECTTYPE(source)
  TYPE IS(DenseSquareMatrixType)
    ALLOCATE(DenseSquareMatrixType :: dest)
  TYPE IS(DenseRectMatrixType)
    ALLOCATE(DenseRectMatrixType :: dest)
  TYPE IS(TriDiagMatrixType)
    ALLOCATE(TriDiagMatrixType :: dest)
  TYPE IS(BandedMatrixType)
    ALLOCATE(BandedMatrixType :: dest)
  TYPE IS(DistributedBandedMatrixType)
    ALLOCATE(DistributedBandedMatrixType :: dest)
  TYPE IS(DistributedBlockBandedMatrixType)
    ALLOCATE(DistributedBlockBandedMatrixType :: dest)
  TYPE IS(SparseMatrixType)
    ALLOCATE(SparseMatrixType :: dest)
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScMatrixType)
    ALLOCATE(PETScMatrixType :: dest)
#endif
#ifdef FUTILITY_HAVE_Trilinos
  TYPE IS(TrilinosMatrixType)
    ALLOCATE(TrilinosMatrixType :: dest)
#endif
  CLASS DEFAULT
    CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
        "Unsupported source matrix type")
  ENDSELECT
  CALL dest%init(params)
ENDSUBROUTINE MatrixResemble
!
!-------------------------------------------------------------------------------
!> @brief Subroutine provides an interface to matrix vector multiplication for
!> the MatrixType.
!> @param trans single character input indicating whether or not to use the
!>        transpose of @c A
!> @param thisMatrix derived matrix type.
!> @param alpha the scalar used to scale @c x
!> @param x the vector to multiply with @c A
!> @param beta the scalar used to scale @c y
!> @param y the vector to add to the product of @c A and @c x
!> @param uplo character indicating if @c thisMatrix is upper or lower triangular
!> @param diag character indicating if diagonal of @c thisMatrix should be treated
!> @param incx_in integer containing distance between elements in @c x
!>
SUBROUTINE matvec_MatrixType(thisMatrix,trans,alpha,x,beta,y,uplo,diag,incx_in)
  CHARACTER(LEN=*),PARAMETER :: myName='matvec_MatrixType'
  CLASS(MatrixType),INTENT(INOUT) :: thisMatrix
  CHARACTER(LEN=1),OPTIONAL,INTENT(IN) :: trans
  REAL(SRK),INTENT(IN),OPTIONAL :: alpha
  REAL(SRK),INTENT(IN) :: x(:)
  REAL(SRK),INTENT(IN),OPTIONAL :: beta
  REAL(SRK),INTENT(INOUT) :: y(:)
  CHARACTER(LEN=1),INTENT(IN),OPTIONAL :: uplo
  CHARACTER(LEN=1),INTENT(IN),OPTIONAL :: diag
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx_in
  !
  CHARACTER(LEN=1) :: t,ul,d
  REAL(SRK) :: a,b
  INTEGER(SIK),ALLOCATABLE :: idxMult(:)
  INTEGER(SIK) :: bIdx,incx

  SELECTTYPE(mat => thisMatrix)
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScMatrixType)
    CALL matvec_PETSc(mat,trans,alpha,x,beta,y,uplo,diag,incx_in)
    RETURN
#endif
  ENDSELECT

  IF(thisMatrix%isInit) THEN
    t='n'
    ul='n'
    d='n'
    incx=1_SIK
    IF(PRESENT(trans)) t=trans
    IF(PRESENT(uplo)) ul=uplo
    IF(PRESENT(diag)) d=diag
    IF(PRESENT(incx_in)) incx=incx_in

    SELECT TYPE(thisMatrix)
    TYPE IS(DenseSquareMatrixType)
      IF(ul /= 'n') THEN
        y=x
        CALL BLAS2_matvec(ul,t,d,thisMatrix%a,y,incx)
      ELSEIF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
        CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%n, &
            alpha,thisMatrix%a,thisMatrix%n,x,1,beta,y,1)
      ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
        CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%n, &
            alpha,thisMatrix%a,thisMatrix%n,x,1,y,1)
      ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
        CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%n, &
            thisMatrix%a,thisMatrix%n,x,1,beta,y,1)
      ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
        CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%n, &
            thisMatrix%a,thisMatrix%n,x,1,y,1)
      ENDIF
    TYPE IS(DenseRectMatrixType)
      IF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
        CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%m, &
            alpha,thisMatrix%a,thisMatrix%n,x,1,beta,y,1)
      ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
        CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%m, &
            alpha,thisMatrix%a,thisMatrix%n,x,1,y,1)
      ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
        CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%m, &
            thisMatrix%a,thisMatrix%n,x,1,beta,y,1)
      ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
        CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%m, &
            thisMatrix%a,thisMatrix%n,x,1,y,1)
      ENDIF
    TYPE IS(SparseMatrixType)
      IF(ul /= 'n') THEN
        y=x
        CALL trsv_sparse(ul,t,d,thisMatrix%a,thisMatrix%ia,thisMatrix%ja,y,incx)
      ELSEIF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
        CALL BLAS2_matvec(thisMatrix%n,thisMatrix%nnz,thisMatrix%ia, &
            thisMatrix%ja,thisMatrix%a,alpha,x,beta,y)
      ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
        CALL BLAS2_matvec(thisMatrix%n,thisMatrix%nnz,thisMatrix%ia, &
            thisMatrix%ja,thisMatrix%a,alpha,x,y)
      ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
        CALL BLAS2_matvec(thisMatrix%n,thisMatrix%nnz,thisMatrix%ia, &
            thisMatrix%ja,thisMatrix%a,x,beta,y)
      ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
        CALL BLAS2_matvec(thisMatrix%n,thisMatrix%nnz,thisMatrix%ia, &
            thisMatrix%ja,thisMatrix%a,x,y)
      ENDIF
    TYPE IS(BandedMatrixType)
      IF(ul /= 'n' .OR. d /= 'n' .OR. t /= 'n' .OR. incx /= 1) THEN
        CALL eMatrixType%raiseError('Incorrect call to '// &
            modName//'::'//myName//' - This interface is not implemented.')
      ENDIF

      a=1.0_SRK
      b=1.0_SRK
      IF(PRESENT(alpha)) a=alpha
      IF(PRESENT(beta)) b=beta

      IF(b == 0.0_SRK) THEN
        y=0.0_SRK
      ELSEIF(b /= 1.0_SRK) THEN
        y=y*b
      ENDIF

      ! This can probably be optimized for the a /= 1 case
      IF(a == 1.0_SRK) THEN
        DO bIdx=1,SIZE(thisMatrix%bandIdx)
          idxMult=thisMatrix%bands(bIdx)%jIdx-thisMatrix%bandIdx(bIdx)
          y(idxMult)=y(idxMult)+thisMatrix%bands(bIdx)%elem*x(thisMatrix%bands(bIdx)%jIdx)
        ENDDO
      ELSEIF(a == -1.0_SRK) THEN
        DO bIdx=1,SIZE(thisMatrix%bandIdx)
          idxMult=thisMatrix%bands(bIdx)%jIdx-thisMatrix%bandIdx(bIdx)
          y(idxMult)=y(idxMult)-thisMatrix%bands(bIdx)%elem*x(thisMatrix%bands(bIdx)%jIdx)
        ENDDO
      ELSEIF(a /= 0.0_SRK) THEN
        DO bIdx=1,SIZE(thisMatrix%bandIdx)
          idxMult=thisMatrix%bands(bIdx)%jIdx-thisMatrix%bandIdx(bIdx)
          y(idxMult)=y(idxMult)+a*thisMatrix%bands(bIdx)%elem*x(thisMatrix%bands(bIdx)%jIdx)
        ENDDO
      ENDIF

    CLASS IS(DistributedBandedMatrixType)
      CALL eMatrixType%raiseError('Incorrect call to '// &
          modName//'::'//myName//' - This interface is not available.')
    CLASS DEFAULT
      CALL eMatrixType%raiseError('Incorrect call to '// &
          modName//'::'//myName//' - This interface is not available.')
    ENDSELECT
  ENDIF
ENDSUBROUTINE matvec_MatrixType
!
!-------------------------------------------------------------------------------
!> @brief Subroutine provides an interface to matrix vector multiplication for
!> the MatrixType.
!> @param trans single character input indicating whether or not to use the
!>        transpose of @c A
!> @param thisMatrix derived matrix type.
!> @param alpha the scalar used to scale @c x
!> @param x the vector to multiply with @c A
!> @param beta the scalar used to scale @c y
!> @param y the vector to add to the product of @c A and @c x
!> @param uplo character indicating if @c thisMatrix is upper or lower triangular
!> @param diag character indicating if diagonal of @c thisMatrix should be treated
!> @param incx_in integer containing distance between elements in @c x
!>
!> TODO: This is more of a mess than it needs to be, and should be resturctured
!> to match the matvec implementation. Split up the functionality and store
!> closer to their respective implementations.
SUBROUTINE matvec_MatrixTypeVectorType(thisMatrix,trans,alpha,x,beta,y,uplo,diag,incx_in)
  CHARACTER(LEN=*),PARAMETER :: myName='matvec_MatrixTypeVectorType'
  CLASS(MatrixType),INTENT(INOUT) :: thisMatrix
  CLASS(VectorType),INTENT(INOUT) :: x
  CHARACTER(LEN=1),OPTIONAL,INTENT(IN) :: trans
  REAL(SRK),INTENT(IN),OPTIONAL :: alpha
  REAL(SRK),INTENT(IN),OPTIONAL :: beta
  CLASS(VectorType),INTENT(INOUT) :: y
  CHARACTER(LEN=1),INTENT(IN),OPTIONAL :: uplo
  CHARACTER(LEN=1),INTENT(IN),OPTIONAL :: diag
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx_in
  CHARACTER(LEN=1) :: t,ul,d
  INTEGER(SIK) :: incx
  REAL(SRK) :: a,b

#ifdef FUTILITY_HAVE_PETSC
  SELECTTYPE(mat => thisMatrix); TYPE IS(PETScMatrixType)
    CALL matvec_PETScVector(mat,trans,alpha,x,beta,y,uplo,diag,incx_in)
    RETURN
  ENDSELECT
#endif

#ifdef FUTILITY_HAVE_Trilinos
  SELECTTYPE(mat => thisMatrix); TYPE IS(TrilinosMatrixType)
    CALL matvec_TrilinosVector(mat,trans,alpha,x,beta,y,uplo,diag,incx_in)
    RETURN
  ENDSELECT
#endif

  IF(thisMatrix%isInit) THEN
    t='n'
    ul='n'
    d='n'
    incx=1_SIK
    a=1
    b=1
    IF(PRESENT(trans)) t=trans
    IF(PRESENT(uplo)) ul=uplo
    IF(PRESENT(diag)) d=diag
    IF(PRESENT(incx_in)) incx=incx_in
    IF(PRESENT(alpha)) a=alpha
    IF(PRESENT(beta))  b=beta

    SELECT TYPE(x)
    TYPE IS(RealVectorType)
      SELECT TYPE(y)
      TYPE IS(RealVectorType)
        SELECT TYPE(thisMatrix)
        TYPE IS(DenseSquareMatrixType)
          IF(ul /= 'n') THEN
            y%b=x%b
            CALL BLAS2_matvec(ul,t,d,thisMatrix%a,y%b,incx)
          ELSEIF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%n, &
                alpha,thisMatrix%a,thisMatrix%n,x%b,1,beta,y%b,1)
          ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%n, &
                alpha,thisMatrix%a,thisMatrix%n,x%b,1,y%b,1)
          ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%n, &
                thisMatrix%a,thisMatrix%n,x%b,1,beta,y%b,1)
          ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%n, &
                thisMatrix%a,thisMatrix%n,x%b,1,y%b,1)
          ENDIF
        TYPE IS(DenseRectMatrixType)
          IF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%m, &
                alpha,thisMatrix%a,thisMatrix%n,x%b,1,beta,y%b,1)
          ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%m, &
                alpha,thisMatrix%a,thisMatrix%n,x%b,1,y%b,1)
          ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%m, &
                thisMatrix%a,thisMatrix%n,x%b,1,beta,y%b,1)
          ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS2_matvec(t,thisMatrix%n,thisMatrix%m, &
                thisMatrix%a,thisMatrix%n,x%b,1,y%b,1)
          ENDIF
        TYPE IS(SparseMatrixType)
          IF(ul /= 'n') THEN
            y%b=x%b
            CALL trsv_sparse(ul,t,d,thisMatrix%a,thisMatrix%ia,thisMatrix%ja,y%b,incx)
          ELSEIF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS2_matvec(thisMatrix%n,thisMatrix%nnz,thisMatrix%ia, &
                thisMatrix%ja,thisMatrix%a,alpha,x%b,beta,y%b)
          ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS2_matvec(thisMatrix%n,thisMatrix%nnz,thisMatrix%ia, &
                thisMatrix%ja,thisMatrix%a,alpha,x%b,y%b)
          ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS2_matvec(thisMatrix%n,thisMatrix%nnz,thisMatrix%ia, &
                thisMatrix%ja,thisMatrix%a,x%b,beta,y%b)
          ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS2_matvec(thisMatrix%n,thisMatrix%nnz,thisMatrix%ia, &
                thisMatrix%ja,thisMatrix%a,x%b,y%b)
          ENDIF
        TYPE IS(BandedMatrixType)
          CALL matvec_MatrixType(thisMatrix,trans=t,alpha=a,X=x%b,beta=b, &
              Y=y%b,uplo=ul,diag=d,incx_in=incx)
        CLASS IS(DistributedBandedMatrixType)
          CALL matvec_MatrixType(thisMatrix,trans=t,alpha=a,X=x%b,beta=b, &
              Y=y%b,uplo=ul,diag=d,incx_in=incx)
        CLASS DEFAULT
          CALL eMatrixType%raiseError('Incorrect call to '// &
              modName//'::'//myName//' - This interface is not available.')
        ENDSELECT
      CLASS DEFAULT
        CALL eMatrixType%raiseError('Incorrect call to '// &
            modName//'::'//myName//' - This interface is not available.')
      ENDSELECT
#ifdef HAVE_MPI
    TYPE IS(NativeDistributedVectorType)
      SELECT TYPE(y)
      TYPE IS(NativeDistributedVectorType)
        SELECT TYPE(thisMatrix)
        CLASS IS(DistributedBandedMatrixType)
          CALL matvec_DistrBandedMatrixType(thisMatrix,x%b,y%b,t,ul,d,incx,a,b)
        CLASS DEFAULT
          CALL eMatrixType%raiseError('Incorrect call to '//modName//'::'// &
              myName//' - This interface is not available.')
        ENDSELECT
      CLASS DEFAULT
        CALL eMatrixType%raiseError('Incorrect call to '// &
            modName//'::'//myName//' - This interface is not available.')
      ENDSELECT
#endif
    CLASS DEFAULT
      CALL eMatrixType%raiseError('Incorrect call to '// &
          modName//'::'//myName//' - This interface is not available.')
    ENDSELECT
  ENDIF
ENDSUBROUTINE matvec_MatrixTypeVectorType

!
!-------------------------------------------------------------------------------
!> @brief Parallel algorithm to perform a matvec operation with banded matrices.
!> @param thisMatrix derived matrix type.
!> @param x the vector to multiply with @c A
!> @param y the vector to add to the product of @c A and @c x
!> @param t single character input indicating whether or not to use the
!>        transpose of @c A
!> @param ul character indicating if @c thisMatrix is upper or lower triangular
!> @param d character indicating if diagonal of @c thisMatrix should be treated
!> @param incx integer containing distance between elements in @c x
!> @param a the scalar used to scale @c x
!> @param b the scalar used to scale @c y
!>
SUBROUTINE matvec_DistrBandedMatrixType(thisMatrix,x,y,t,ul,d,incx,a,b)
  CLASS(DistributedBandedMatrixType),INTENT(INOUT) :: thisMatrix
  REAL(SRK),INTENT(INOUT) :: x(:)
  CHARACTER(LEN=1),INTENT(IN) :: t
  REAL(SRK),INTENT(IN) :: a
  REAL(SRK),INTENT(IN) :: b
  REAL(SRK),INTENT(INOUT) :: y(:)
  CHARACTER(LEN=1),INTENT(IN) :: ul
  CHARACTER(LEN=1),INTENT(IN) :: d
  INTEGER(SIK),INTENT(IN) :: incx
  INTEGER(SIK) :: rank,k,lowIdx,highIdx
  REAL(SRK),ALLOCATABLE :: recvResult(:,:),sendResult(:,:),tmpProduct(:)
#ifdef HAVE_MPI
  ! Get rank
  INTEGER(SIK) :: sendRequests(MATVEC_SLOTS),recvRequests(MATVEC_SLOTS)
  INTEGER(SIK) :: ctRecv(MATVEC_SLOTS),srcRank,destRank
  INTEGER(SIK) :: i,idxTmp,ctDefault,mpierr,nproc
  CALL MPI_Comm_rank(thisMatrix%comm,rank,mpierr)
  CALL MPI_Comm_Size(thisMatrix%comm,nProc,mpierr)
  sendRequests=MPI_REQUEST_NULL
  recvRequests=MPI_REQUEST_NULL
#else
  rank=0
#endif
  ! NOTE: incx,ul,d,t are NOT USED

  ! We will store our send/recv arrays in an array of length MATVEC_SLOTS
  ! This allows one to be used for computation and the rest for
  ! communication.
  ! Which one to write to will be determined by *Counter MOD MATVEC_SLOTS
  ! The recv/sendResult array will sometimes hold it's full length,
  ! and sometimes less.
  ALLOCATE(recvResult(thisMatrix%iOffsets(2),MATVEC_SLOTS))
  ALLOCATE(sendResult(thisMatrix%iOffsets(2),MATVEC_SLOTS))
  ALLOCATE(tmpProduct(thisMatrix%iOffsets(rank+2)- &
      thisMatrix%iOffsets(rank+1)))
  tmpProduct=0.0_SRK

#ifdef HAVE_MPI
  ! On each rank, loop over the chunks held (on diagonal moving down)
  DO i=1,SIZE(thisMatrix%iOffsets)-2
    ! Destination rank:
    destRank = MODULO(rank+i,nProc)
    ctDefault=thisMatrix%iOffsets(destRank+2)-thisMatrix%iOffsets(destRank+1)

    ! First do local computation and send
    IF (thismatrix%chunks(destRank+1)%isInit) THEN
      ! Decide whether to send whole vector or multiple sparse
      IF (thisMatrix%chunks(destRank+1)%nnz >= thisMatrix%chunks(destRank+1)%n) THEN
        ! Check if we can safely write to sendRequests
        CALL pop_send(sendResult,idxTmp,sendRequests)
        CALL BLAS_matvec(THISMATRIX=thisMatrix%chunks(destRank+1),X=x,&
            y=sendResult(1:ctDefault,idxTmp),ALPHA=1.0_SRK,BETA=0.0_SRK)
        CALL MPI_ISend(sendResult(1:ctDefault,idxTmp),ctDefault, &
            MPI_DOUBLE_PRECISION,destRank,0,thisMatrix%comm,sendRequests(idxTmp),mpierr)
      ELSE
        ! Gather and send several sparse vectors
        CALL pop_send(sendResult,idxTmp,sendRequests)
        lowIdx=1
        DO k=1,SIZE(thisMatrix%chunks(destRank+1)%bandIdx)
          highidx=lowIdx+SIZE(thisMatrix%chunks(destRank+1)%bands(k)%jIdx)-1
          sendResult(lowIdx:highIdx,idxTmp)=thisMatrix%chunks(destRank+1)%bands(k)%elem &
              *x(thisMatrix%chunks(destRank+1)%bands(k)%jIdx)
          lowIdx=highIdx+1
        ENDDO
        CALL MPI_ISend(sendResult(1:highIdx,idxTmp),highIdx,MPI_DOUBLE_PRECISION, &
            destRank,0,thisMatrix%comm,sendRequests(idxTmp), mpierr)
      ENDIF
    ENDIF
    ! We might receive data from rank MOD(rank-i,nproc)
    srcRank = MODULO(rank-i,nProc)
    IF (thisMatrix%incIdxStt(srcRank+1) == -1) THEN
      ! We are receiving a whole vector at once
      ! If we've filled up the available storage, we need
      ! to wait for communication to finish up
      CALL pop_recv(tmpProduct,recvResult,thisMatrix, &
          ctRecv,idxTmp,recvRequests)
      ctRecv(idxTmp)=SIZE(tmpProduct)
      CALL MPI_IRecv(recvResult(1:ctRecv(idxTmp),idxTmp), ctRecv(idxTmp), &
          MPI_DOUBLE_PRECISION,srcRank,0,thisMatrix%comm, &
          recvRequests(idxTmp),mpierr)
    ELSEIF (thisMatrix%incIdxStt(srcRank+1) > 0) THEN
      ! We are receiving multiple sparse chunks gathered together
      CALL pop_recv(tmpProduct,recvResult,thisMatrix, &
          ctRecv,idxTmp,recvRequests)
      ctRecv(idxTmp)= -srcRank-1
      lowIdx=thisMatrix%incIdxStt(srcRank+1)
      highIdx=thisMatrix%incIdxStp(srcRank+1)
      CALL MPI_IRecv(recvResult(1:(highIdx-lowIdx+1),idxTmp),highIdx-lowIdx+1, &
          MPI_DOUBLE_PRECISION,srcRank,0,thisMatrix%comm,recvRequests(idxTmp), mpierr)
    ENDIF
  ENDDO
#endif
  ! Now, take care of locally held data.
  SELECT TYPE(thisMatrix)
  TYPE IS(DistributedBlockBandedMatrixType)
    IF(thisMatrix%chunks(rank+1)%isInit) THEN
      CALL BLAS_matvec(THISMATRIX=thisMatrix%chunks(rank+1),X=x, &
          y=tmpProduct,ALPHA=1.0_SRK,BETA=1.0_SRK)
    ENDIF
    IF(.NOT. thisMatrix%blockMask) THEN
      DO k=1,thisMatrix%nlocalBlocks
        lowIdx=(k-1)*thisMatrix%blockSize+1
        highIdx=lowIdx-1+thisMatrix%blockSize
        CALL matvec_MatrixType(THISMATRIX=thisMatrix%blocks(k),X=x(lowIdx:highIdx), &
            Y=tmpProduct(lowIdx:highIdx),ALPHA=1.0_SRK,BETA=1.0_SRK)
      ENDDO
    ENDIF
  TYPE IS(DistributedBandedMatrixType)
    IF(thisMatrix%chunks(rank+1)%isInit) THEN
      CALL BLAS_matvec(THISMATRIX=thisMatrix%chunks(rank+1),X=x, &
          y=tmpProduct,ALPHA=1.0_SRK,BETA=1.0_SRK)
    ENDIF
  ENDSELECT
#if HAVE_MPI
  ! Wait for remaining requests to finish:
  DO k=1,MATVEC_SLOTS
    idxTmp=k
    CALL pop_recv(tmpProduct,recvResult,thisMatrix,ctRecv,idxTmp,recvRequests,.TRUE.)
  ENDDO
  DO k=1,MATVEC_SLOTS
    idxTmp=k
    CALL pop_send(sendResult,idxTmp,sendRequests,.TRUE.)
  ENDDO
#endif

  ! do y=alpha*reduce+beta*y
  ! We take many special cases in order to effectively use FMA instructions
  ! or reduce the required flops.
  IF(a == 1.0_SRK) THEN
    IF(b == 1.0_SRK) THEN
      y=tmpProduct+y
    ELSEIF(b == 0.0_SRK) THEN
      y=tmpProduct
    ELSEIF(b == -1.0_SRK) THEN
      y=tmpProduct-y
    ELSE
      y=tmpProduct+b*y
    ENDIF
  ELSEIF(a == 0.0_SRK) THEN
    IF(b == 0.0_SRK) THEN
      y=0.0_SRK
    ELSEIF(b == -1.0_SRK) THEN
      y=-y
    ELSEIF(b /= 1.0_SRK) THEN
      y=b*y
    ENDIF
  ELSEIF(a == -1.0_SRK) THEN
    IF(b == 1.0_SRK) THEN
      y=y-tmpProduct
    ELSEIF(b == 0.0_SRK) THEN
      y=-tmpProduct
    ELSEIF(b == -1.0_SRK) THEN
      y=-tmpProduct-y
    ELSE
      y=b*y-tmpProduct
    ENDIF
  ELSE
    IF(b == 1.0_SRK) THEN
      y=a*tmpProduct+y
    ELSEIF(b == 0.0_SRK) THEN
      y=a*tmpProduct
    ELSEIF(b == -1.0_SRK) THEN
      y=a*tmpProduct-y
    ELSE
      y=a*tmpProduct+b*y
    ENDIF
  ENDIF
ENDSUBROUTINE matvec_DistrBandedMatrixType

#ifdef HAVE_MPI
!
!-------------------------------------------------------------------------------
!> @brief Helper routine that handles the movement of data out of the
!>        recieving buffer. This buffer has multiple "slots" that are either
!>        being worked on directly or receiving data from MPI in background.
!>        Empty slots are denoted with MPI_REQUEST_NULL
!> @param acc The accumulating vector
!> @param valBuf The recv buffer for values
!> @param thisMat The matrixtype used in multiplication
!> @param ctBuf Count of elements in each buffer slot
!> @param idx The empty buffer slot that was popped
!> @param req The array of active requests for data
!> @param f Optional bool to flush all requests
!>
SUBROUTINE pop_recv(acc,valBuf,thisMat,ctBuf,idx,req,f)
  !> Local vector data
  REAL(SRK), INTENT(INOUT) :: acc(:)
  !> List of buffers
  REAL(SRK), INTENT(INOUT) :: valBuf(:,:)
  !> Matrix used in SpMV
  CLASS(DistributedBandedMatrixType),INTENT(INOUT) :: thisMat
  !> Array of buffer sizes
  INTEGER(SIK), INTENT(INOUT) :: ctBuf(MATVEC_SLOTS)
  !> The buffer slot popped
  INTEGER(SIK), INTENT(OUT) :: idx
  !> List of MPI Requests
  INTEGER(SIK), INTENT(INOUT) :: req(MATVEC_SLOTS)
  !> Optional argument to flush all buffers
  LOGICAL(SBK),OPTIONAL,INTENT(IN) :: f
  INTEGER(SIK) :: stt,stp,mpierr,i,rank
  LOGICAL(SBK) :: force
  CALL MPI_Comm_Rank(thisMat%comm,rank,mpierr)

  force=.FALSE.
  IF (PRESENT(f)) force=f
  IF (force) THEN
    IF (req(idx) == MPI_REQUEST_NULL) RETURN
    CALL MPI_Wait(req(idx),MPI_STATUS_IGNORE,mpierr)
  ELSE
    IF (ANY(req == MPI_REQUEST_NULL)) THEN
      DO i=1,SIZE(req)
        IF (req(i) == MPI_REQUEST_NULL) THEN
          idx=i
          RETURN
        ENDIF
      ENDDO
    ELSE
      CALL MPI_WaitAny(SIZE(req),req,idx,MPI_STATUS_IGNORE,mpierr)
    ENDIF
  ENDIF

  IF(ctBuf(idx) > 0) THEN
    acc(1:ctBuf(idx))=acc(1:ctBuf(idx))+valBuf(1:ctBuf(idx),idx)
  ELSE
    stt=thisMat%incIdxStt(-ctBuf(idx))
    stp=thisMat%incIdxStp(-ctBuf(idx))
    DO i=stt,stp
      acc(thisMat%incIdxMap(i))=acc(thisMat%incIdxMap(i))+valBuf(i-stt+1,idx)
    ENDDO
  ENDIF
ENDSUBROUTINE pop_recv

!
!-------------------------------------------------------------------------------
!> @brief Helper routine that keeps track of data in the send buffer while
!>        MPI is still working. This buffer has multiple "slots" that are
!>        either being worked on directly or contain data being sent by
!>        MPI in background. Empty slots are denoted with MPI_REQUEST_NULL
!> @param valBuf The send buffer for values
!> @param idx The buffer slot popped
!> @param req The array of requests for data
!> @param f Optional bool to flush all requests
!>
SUBROUTINE pop_send(valBuf,idx,req,f)
  !> Set of data buffers
  REAL(SRK), INTENT(INOUT) :: valBuf(:,:)
  !> The buffer slot popped
  INTEGER(SIK), INTENT(OUT) :: idx
  !> The array of MPI requests
  INTEGER(SIK), INTENT(INOUT) :: req(MATVEC_SLOTS)
  !> Optional argument to flush all buffers
  LOGICAL(SBK),OPTIONAL,INTENT(IN) :: f
  INTEGER(SIK) :: mpierr,i
  LOGICAL(SBK) :: force

  force=.FALSE.
  IF (PRESENT(f)) force=f
  IF (force) THEN
    IF (req(idx) == MPI_REQUEST_NULL) RETURN
    CALL MPI_Wait(req(idx),MPI_STATUS_IGNORE,mpierr)
  ELSE
    IF (ANY(req == MPI_REQUEST_NULL)) THEN
      DO i=1,SIZE(req)
        IF (req(i) == MPI_REQUEST_NULL) THEN
          idx=i
          RETURN
        ENDIF
      ENDDO
    ELSE
      CALL MPI_WaitAny(SIZE(req),req,idx,MPI_STATUS_IGNORE,mpierr)
    ENDIF
  ENDIF
ENDSUBROUTINE pop_send
#endif
!
!-------------------------------------------------------------------------------
!> @brief Subroutine solves a sparse triangular matrix linear system.
!> @param uplo single character input indicating if an upper (U) or lower (L)
!>        maxtrix is stored in @c A
!> @param trans single character input indicating whether or not to use the
!>        transpose of @c A
!> @param diag single character input indicating whether or not a unity
!>        diagonal is used
!> @param a the double-precision matrix to multiply with @c x, stored in CSR format
!> @param ia the indices of the @c a for the first element in each row
!> @param ja the column indices for each element in @c a
!> @param x the double-precision vector to multiply with @c a
!> @param incx_in the increment to use when looping over elements in @c x
PURE SUBROUTINE trsv_sparse(uplo,trans,diag,a,ia,ja,x,incx_in)
  CHARACTER(LEN=1),INTENT(IN) :: uplo
  CHARACTER(LEN=1),INTENT(IN) :: trans
  CHARACTER(LEN=1),INTENT(IN) :: diag
  REAL(SRK),INTENT(IN) :: a(:)
  INTEGER(SIK),INTENT(IN) :: ia(:)
  INTEGER(SIK),INTENT(IN) :: ja(:)
  REAL(SRK),INTENT(INOUT) :: x(:)
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx_in

  LOGICAL(SBK) :: nounit
  INTEGER(SIK) :: i,ix,j,jx,kx,n,incx
  REAL(SRK),PARAMETER :: ZERO=0.0_SRK
  INTRINSIC MAX
  n=SIZE(ia)-1
  IF(PRESENT(incx_in)) THEN
    incx=incx_in
  ELSE
    incx=1_SIK ! This won't ever be used in the current interfaces
  ENDIF
  ! Check inputs
  IF((trans == 't' .OR. trans == 'T' .OR. trans == 'c' .OR. trans == 'C' .OR. &
      trans == 'n' .OR. trans == 'N') .AND. &
      (uplo == 'u' .OR. uplo == 'U' .OR. uplo == 'l' .OR. uplo == 'L') .AND. &
      (diag == 'u' .OR. diag == 'U' .OR. diag == 'n' .OR. diag == 'N')) THEN

    ! Check whether diagonal is treated as unity or not.
    nounit=.FALSE.
    IF(diag == 'n' .OR. diag == 'N') nounit=.TRUE.
    ! Determine how elements of x are stored, and find "first" element accordingly
    IF(incx <= 0_SIK) THEN
      kx=1-(n-1)*incx ! Elements stored in reverse order (highest index to lowest)
    ELSEIF(incx /= 1) THEN
      kx=1 ! Elements stored from lowest index to highest
    ENDIF

    ! Don't use transpose of matrix
    IF(trans == 'n' .OR. trans == 'N') THEN  ! Form  x := inv( A )*x.
      ! Multiply by upper triangular part of matrix
      IF(uplo == 'u' .OR. uplo == 'U') THEN
        ! Elements of x are stored in increments of 1
        IF(incx == 1) THEN
          DO i=n,1,-1
            DO j=ia(i+1)-1,ia(i),-1
              IF(ja(j) == i) EXIT
              x(i)=x(i)-a(j)*x(ja(j))
            ENDDO
            IF(nounit) x(i)=x(i)/a(j)
          ENDDO
        ! Elements are stored in reverse order and/or by increments other than 1
        ELSE
          ix=ABS(incx)*n-kx
          DO i=n,1,-1
            jx=ABS(incx)*n-kx
            DO j=ia(i+1)-1,ia(i),-1
              IF(ja(j) <= i) EXIT
              x(ix)=x(ix)-a(j)*x(jx)
              IF(j == ia(i+1)-1) THEN
                jx=jx-incx*(n-ja(j)+1)
              ELSE
                jx=jx-incx*(ja(j+1)-ja(j))
              ENDIF
            ENDDO
            IF(nounit) x(ix)=x(ix)/a(j)
            ix=ix-incx
          ENDDO
        ENDIF
      ! Multiply by lower triangular part of matrix
      ELSE
        ! Elements of x are stored in increments of 1
        IF(incx == 1) THEN
          DO i=1,n
            DO j=ia(i),ia(i+1)-1
              IF(ja(j) >= i) EXIT
              x(i)=x(i)-a(j)*x(ja(j))
              IF(j == SIZE(ja)) EXIT
            ENDDO
            IF((ja(j) == i) .AND. nounit) x(i)=x(i)/a(j)
          ENDDO
        ! Elements are stored in reverse order and/or by increments other than 1
        ELSE
          ix=kx
          DO i=1,n
            jx=kx
            DO j=ia(i),ia(i+1)-1
              IF(ja(j) >= i) EXIT
              IF(ja(j) > 1) THEN
                jx=jx+incx*(ja(j)-ja(j-1))
              ELSE
                jx=jx+incx*(ja(j)-1)
              ENDIF
              x(ix)=x(ix)-a(j)*x(jx)
              IF(j == SIZE(ja)) EXIT
            ENDDO
            IF(nounit) x(ix)=x(ix)/a(j)
            ix=ix+incx
          ENDDO
        ENDIF
      ENDIF
    ! Use transpose of matrix
    ELSE
      ! Multiply by upper trianbular part of matrix
      IF(uplo == 'u' .OR. uplo == 'U') THEN
        IF(incx == 1) THEN
          DO j=1,SIZE(ia)-1
            IF(.NOT.(x(j) .APPROXEQA. ZERO)) THEN
              IF(nounit) x(j)=x(j)/a(ia(j))
              DO i=ia(j)+1,ia(j+1)-1
                IF(ja(i) <= j) CYCLE
                x(ja(i))=x(ja(i))-x(j)*a(i)
              ENDDO
            ENDIF
          ENDDO
        ! Elements are stored in reverse order and/or by increments other than 1
        ELSE
          ix=kx
          DO j=1,SIZE(ia)-1
            IF(.NOT.(x(ix) .APPROXEQA. ZERO)) THEN
              IF(nounit) x(ix)=x(ix)/a(ia(j))
              jx=ix
              DO i=ia(j)+1,ia(j+1)-1
                IF(ja(i) <= j) CYCLE
                jx=jx+incx*(ja(i)-ja(i-1))
                x(jx)=x(jx)-x(ix)*a(i)
              ENDDO
            ENDIF
            ix=ix+incx
          ENDDO
        ENDIF
      ! Multiply by lower triangular part of matrix
      ELSE
        ! Elements of x are stored in increments of 1
        IF(incx == 1) THEN
          DO i=n,1,-1
            IF(.NOT.(x(i) .APPROXEQA. ZERO)) THEN
              IF(nounit) x(i)=x(i)/a(ia(i+1)-1)
              DO j=ia(i+1)-2,ia(i),-1
                IF(ja(j) > i) CYCLE
                x(ja(j))=x(ja(j))-a(j)*x(i)
              ENDDO
            ENDIF
          ENDDO
        ! Elements are stored in reverse order and/or by increments other than 1
        ELSE
          ix=ABS(incx)*n-kx
          DO i=n,1,-1
            IF(.NOT.(x(ix) .APPROXEQA. ZERO)) THEN
              DO j=ia(i+1)-1,ia(i),-1
                IF(ja(j) == i) EXIT
              ENDDO
              IF(nounit) x(ix)=x(ix)/a(j)
              jx=ABS(incx)*n-kx
              DO j=ia(i+1)-1,ia(i),-1
                IF(j == ia(i+1)-1) THEN
                  jx=jx-incx*(n-ja(j))
                ELSE
                  jx=jx-incx*(ja(j+1)-ja(j))
                ENDIF
                IF(ja(j) >= i) CYCLE
                x(jx)=x(jx)-a(j)*x(ix)
              ENDDO
            ENDIF
            ix=ix-incx
          ENDDO
        ENDIF
      ENDIF
    ENDIF
  ENDIF
ENDSUBROUTINE trsv_sparse
!
!-------------------------------------------------------------------------------
!> @brief Subroutine provides an interface to matrix matrix multiplication for
!> the MatrixType. alpha*A*B+beta*C
!> @param A a derived matrix type to multiply by @c B
!> @param B a derived matrix type to multiply by @c A
!> @param C a derived matrix type to add to the product of @c A and @c B
!>          If C is not initialized when the subroutine is called, it will be
!>          initialized and set to all zeros
!> @param transA single character input indicating whether or not to use the
!>        transpose of @c A
!> @param alpha the scalar used to scale the product of @c A and @cB
!> @param transB single character input indicating whether or not to use the
!>        transpose of @c B
!> @param beta the scalar used to scale @c C
!>
SUBROUTINE matmult_MatrixType(A,B,C,alpha,beta,transA,transB)
  CLASS(MatrixType),INTENT(INOUT) :: A
  CLASS(MatrixType),INTENT(INOUT) :: B
  CLASS(MatrixType),INTENT(INOUT) :: C
  REAL(SRK),INTENT(IN),OPTIONAL :: alpha
  REAL(SRK),INTENT(IN),OPTIONAL :: beta
  CHARACTER(LEN=1),OPTIONAL,INTENT(IN) :: transA
  CHARACTER(LEN=1),OPTIONAL,INTENT(IN) :: transB
  CHARACTER(LEN=1) :: tA
  CHARACTER(LEN=1) :: tB
#ifdef FUTILITY_HAVE_PETSC
  CHARACTER(LEN=*),PARAMETER :: myName='matmult_MatrixType'
  SELECTTYPE(A)
  TYPE IS(PETScMatrixType)
    SELECTTYPE(B)
    TYPE IS(PETScMatrixType)
      SELECTTYPE(C)
      TYPE IS(PETScMatrixType)
        CALL matmult_PETSc(A,B,C,alpha,beta,transA,transB)
      CLASS DEFAULT
        CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
            "Interface not implemented")
      ENDSELECT
    CLASS DEFAULT
      CALL eMatrixType%raiseError(modName//"::"//myName//" - "// &
          "Interface not implemented")
    ENDSELECT
    RETURN
  ENDSELECT
#endif

  IF(A%isInit) THEN
    tA='n'
    IF(PRESENT(transA)) tA=transA
  ENDIF
  IF(B%isInit) THEN
    tB='n'
    IF(PRESENT(transB)) tB=transB
  ENDIF

  IF(A%isInit .AND. B%isInit .AND. C%isInit) THEN
    SELECTTYPE(A)
    TYPE IS(DenseSquareMatrixType)
      SELECTTYPE(B)
      TYPE IS(DenseSquareMatrixType)
        SELECTTYPE(C)
        TYPE IS(DenseSquareMatrixType)
          ! A: Square   B: Square  C: Square
          IF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,alpha,A%a,B%a,beta,C%a)
          ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,alpha,A%a,B%a,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,A%a,B%a,beta,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,A%a,B%a,C%a)
          ENDIF
        TYPE IS(DenseRectMatrixType)
          ! A: Square   B: Square  C: Rect
          IF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,alpha,A%a,B%a,beta,C%a)
          ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,alpha,A%a,B%a,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,A%a,B%a,beta,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,A%a,B%a,C%a)
          ENDIF
        TYPE IS(SparseMatrixType)
          ! NOT SUPPORTED
        ENDSELECT
      TYPE IS(DenseRectMatrixType)
        SELECTTYPE(C)
        TYPE IS(DenseSquareMatrixType)
          ! A: Square   B: Rect  C: Square
          IF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,alpha,A%a,B%a,beta,C%a)
          ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,alpha,A%a,B%a,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,A%a,B%a,beta,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,A%a,B%a,C%a)
          ENDIF
        TYPE IS(DenseRectMatrixType)
          ! A: Square   B: Rect  C: Rect
          IF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,alpha,A%a,B%a,beta,C%a)
          ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,alpha,A%a,B%a,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,A%a,B%a,beta,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,A%a,B%a,C%a)
          ENDIF
        TYPE IS(SparseMatrixType)
          ! NOT SUPPORTED
        ENDSELECT
      TYPE IS(SparseMatrixType)
        ! NOT SUPPORTED
      ENDSELECT
    TYPE IS(DenseRectMatrixType)
      SELECTTYPE(B)
      TYPE IS(DenseSquareMatrixType)
        SELECTTYPE(C)
        TYPE IS(DenseSquareMatrixType)
          ! A: Rect   B: Square  C: Square
          IF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,alpha,A%a,B%a,beta,C%a)
          ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,alpha,A%a,B%a,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,A%a,B%a,beta,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,A%a,B%a,C%a)
          ENDIF
        TYPE IS(DenseRectMatrixType)
          ! A: Rect   B: Square  C: Rect
          IF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,alpha,A%a,B%a,beta,C%a)
          ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,alpha,A%a,B%a,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,A%a,B%a,beta,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,A%a,B%a,C%a)
          ENDIF
        TYPE IS(SparseMatrixType)
          ! NOT SUPPORTED
        ENDSELECT
      TYPE IS(DenseRectMatrixType)
        SELECTTYPE(C)
        TYPE IS(DenseSquareMatrixType)
          ! A: Rect   B: Rect  C: Square
          IF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,alpha,A%a,B%a,beta,C%a)
          ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,alpha,A%a,B%a,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,A%a,B%a,beta,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%n,B%n,A%a,B%a,C%a)
          ENDIF
        TYPE IS(DenseRectMatrixType)
          ! A: Rect   B: Rect  C: Rect
          IF(PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,alpha,A%a,B%a,beta,C%a)
          ELSEIF(PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,alpha,A%a,B%a,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,A%a,B%a,beta,C%a)
          ELSEIF(.NOT.PRESENT(alpha) .AND. .NOT.PRESENT(beta)) THEN
            CALL BLAS3_matmult(tA,tB,C%n,C%m,B%n,A%a,B%a,C%a)
          ENDIF
        TYPE IS(SparseMatrixType)
          ! NOT SUPPORTED
        ENDSELECT
      TYPE IS(SparseMatrixType)
        ! NOT SUPPORTED
      ENDSELECT
    TYPE IS(SparseMatrixType)
      ! NOT SUPPORTED
    ENDSELECT
  ENDIF
ENDSUBROUTINE matmult_MatrixType
!
ENDMODULE MatrixTypes

max_stars_repo_name: picmc/Futility

max_issues_repo_name: picmc/Futility

max_forks_repo_name: picmc/Futility


text:
Describe Users/Samantha here.

20070719 18:46:33 nbsp Welcome to the Wiki. Please consider looking at both Welcome to the Wiki/Business Owner and Reviews. Users/JasonAller

20070719 22:57:30 nbsp Hi! I wanted to add to something to what Users/Davidlm mentioned in his last message to you. As a loyal client, you would be one of the best people to expand on the services Spa Central offers (currently there are only a few pictures and not much more). Also, yes, there may be bitingif not vindictivecomments on many pages, but if I am like some other dwiki users, a few things about Users/ElisaWellers comment will stand out: first, its very angry (and not just the Oh, my soup was cold and the waiter forgot to refill my water kind of angry). Second, shes a new user (blank User page). This combo would prompt me to check her other http://daviswiki.org/Users/ElisaWeller?actionuserinfo contributions where we find. . . well, there are no others. Just two attacks against one business. Granted, not everyone will check every dubious comment (though it only takes a few clicks to do), but some may. For the rest, a nice positive Reviews review would likely suffice. Users/KevinChin

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*DIM, analysisFilename, STRING, 200
analysisFilename(1) = '202_HDPE_2Div_fundering_multiple'

! GeometryFile(1) = 'KLP_Rechthoek'
! SleeperType = '202'

Materiaal_S = Kunststof_HDPE
Materiaal_W = Staal
Materiaal_RH = Staal
Materiaal_Rail = Staal
AantalSleepers = 10

LoadLocation_X = 0
W_Divisions = 2

MakeUseOfSymmetry = 'false'
*DEL, Sleepers, NOPR
*DIM, Sleepers, Array, 6, AantalSleepers

! Tabel met data voor elke individuele ligger. 1) type-model [1 = KLP_Rechthoek, 2 = KLP_Rechthoek_Simple] 2) type-ligger, 3) divisions, 4) foundation
Sleepers(1,1) = 1, 202, W_Divisions, K_foundation, fundering_OndersteundPercentage, fundering_MultiplierOndersteund

!Octagonal reinforcement from here on:
Sleepers(1,2) = 2, 202, 2, K_foundation, fundering_OndersteundPercentage, fundering_MultiplierOndersteund
Sleepers(1,3) = 2, 202, 2, K_foundation, fundering_OndersteundPercentage, fundering_MultiplierOndersteund
Sleepers(1,4) = 2, 202, 1, K_foundation, fundering_OndersteundPercentage, fundering_MultiplierOndersteund
Sleepers(1,5) = 2, 202, 1, K_foundation, fundering_OndersteundPercentage, fundering_MultiplierOndersteund
Sleepers(1,6) = 2, 202, 1, K_foundation, fundering_OndersteundPercentage, fundering_MultiplierOndersteund
Sleepers(1,7) = 2, 202, 1, K_foundation, fundering_OndersteundPercentage, fundering_MultiplierOndersteund
Sleepers(1,8) = 2, 202, 1, K_foundation, fundering_OndersteundPercentage, fundering_MultiplierOndersteund
Sleepers(1,9) = 2, 202, 1, K_foundation, fundering_OndersteundPercentage, fundering_MultiplierOndersteund
Sleepers(1,10) = 2, 202, 1, K_foundation, fundering_OndersteundPercentage, fundering_MultiplierOndersteund

max_stars_repo_name: dryjognos/AnsysSleeperModel

max_issues_repo_name: dryjognos/AnsysSleeperModel

max_forks_repo_name: dryjognos/AnsysSleeperModel