text:
!=========================================================================================
!Copyright (c) 2009-2019, The Regents of the University of Massachusetts, Amherst.
!E. Polizzi research lab
!All rights reserved.
!
!Redistribution and use in source and binary forms, with or without modification,
!are permitted provided that 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. Neither the name of the University nor the names of its contributors may be used to endorse or promote
!    products derived from this software without specific prior written permission.
!
!THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 AUTHOR 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 THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!==========================================================================================

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!! SPARSE PRIMITIVES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!! SINGLE PRECISION VERSION

!!!!! List of routines

!saddcsr, caddcsr           :: add two square csr matrices and return a new one (C=alpha*A+beta*B)
!csaddcsr                   :: same as zaddcsr but A, B are real;  alpha, beta, C are complex
!cinc_addcsr                :: Increment matrix B by adding matrix A times alpha (sparsity pattern of A included in B)

!scsr2csr_up,  ccsr2csr_up  :: extract upper triangular part of a full csr matrix
!scsr2csr_low, ccsr2csr_low :: extract lower triangular part of a full csr matrix

!scsrmm, ccsrmm             :: mat-vec multiplication =>  b=alpha*A*x+beta*b
!chcsrmm                    :: mat-vec multiplication =>  b=alpha*A*x+beta*b (A is complex Hermitian if UPLO=L,U)
!sccsrmm                    :: same as zcsrmm but A is real;  alpha, beta, b are complex

!scsr_transpose, ccsr_transpose :: transpose a csr matrix
!ccsr_htranspose                :: transpose conjugate a csr matrix

!scsr_uplo_to_csr, ccsr_uplo_to_csr  :: convert a upper or lower csr matrix into a full csr one
!chcsr_uplo_to_csr                   :: convert a upper or lower (Hermitian) csr matrix into a full csr one

!scsr_convert_upper                 :: convert a upper,lower or full csr matrix into an upper csr matrix
!chcsr_convert_full                 :: convert a upper,lower or full csr Hermitian matrix into a full csr matrix

!scoo2csr :: transform coordinate format into csr format
!ccoo2csr :: transform coordinate format into csr format

!ssubmat                            :: extract csr submatrix from csr matrix
!csubmat                            :: extract csr submatrix from csr matrix

subroutine saddcsr(N,M,opt,alpha,sa,isa,jsa,beta,sb,isb,jsb,sc,isc,jsc)
  !  Purpose
  !  =======
  !
  !  Addition of two (NxM any M) sparse CSR matrices A and B and return a new (NxM) CSR matrix C
  !  C=alpha*A+beta*B
  !
  !  A,B,C are real
  !
  !  Arguments
  !  =========
  !
  !  N,M    (input) INTEGER
  !         The number of rows/columns of the matrices A,B,C.  N,M >= 0.
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs addition, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine computes only values of the array isc of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value ic(n+1) - 1 is the actual
  !         number of the elements in the arrays sc and jsc.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array jsc
  !
  !         if opt=3, the routine has been called previously with the parameter opt=2. On exit it will
  !         return the array sc
  !
  !
  !  alpha  (input) REAL
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa REAL)
  !
  !  beta  (input) REAL
  !
  !  sb,isb,jsb (input) CSR format for the matrix B (sb REAL)
  !
  !  sc,isc,jsc (input/output) CSR format for the matrix C (sc REAL).
  !             see definition of 'opt' to determine input or output.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  real,dimension(*) :: sa,sb,sc
  integer,dimension(*) :: isa,isb,isc,jsa,jsb,jsc
  real :: alpha,beta
  integer :: opt,N,M

  integer :: i,ka,kb,kc,ja,jb
  logical :: next,row,col,val

  row=.false.
  col=.false.
  val=.false.

  if (opt==0) then
     row=.true.
     col=.true.
     val=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     col=.true.
  elseif (opt==3) then
     val=.true.
  endif

  kc=1
  if (row) isc(1)=kc
  do i=1,N
     ka=isa(i)
     kb=isb(i)

     next=.false.
     do while (.not.next)

        if (ka<=isa(i+1)-1) then
           ja=jsa(ka)
        else
           ja=m+1
        endif

        if (kb<=isb(i+1)-1) then
           jb=jsb(kb)
        else
           jb=m+1
        endif

        if ((ja==jb).and.(ja/=M+1)) then
           if (val) sc(kc)=alpha*sa(ka)+beta*sb(kb)
           if (col) jsc(kc)=ja
           ka=ka+1
           kb=kb+1
           kc=kc+1
        elseif (ja<jb)  then
           if (val) sc(kc)=alpha*sa(ka)
           if (col) jsc(kc)=ja
           ka=ka+1
           kc=kc+1
        elseif (jb<ja)  then
           if (val) sc(kc)=beta*sb(kb)
           if (col) jsc(kc)=jb
           kb=kb+1
           kc=kc+1
        end if

        if ((ka>isa(i+1)-1).and.(kb>isb(i+1)-1)) next=.true.

     end do

     if (row) isc(i+1)=kc

  end do

end subroutine saddcsr

subroutine caddcsr(N,M,opt,alpha,sa,isa,jsa,beta,sb,isb,jsb,sc,isc,jsc)
  !  Purpose
  !  =======
  !
  !  Addition of two (NxM any M) sparse CSR matrices A and B and return a new (NxM) CSR matrix C
  !  C=alpha*A+beta*B
  !
  !  A,B,C are complex
  !
  !  Arguments
  !  =========
  !
  !  N,M    (input) INTEGER
  !         The number of rows/columns of the matrices A,B,C.  N,M >= 0.
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs addition, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine computes only values of the array isc of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value ic(n+1) - 1 is the actual
  !         number of the elements in the arrays sc and jsc.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array jsc
  !
  !         if opt=3, the routine has been called previously with the parameter opt=2. On exit it will
  !         return the array sc
  !
  !
  !  alpha  (input) COMPLEX
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa COMPLEX)
  !
  !  beta  (input) COMPLEX
  !
  !  sb,isb,jsb (input) CSR format for the matrix B (sb COMPLEX)
  !
  !  sc,isc,jsc (input/output) CSR format for the matrix C (sc COMPLEX).
  !             see definition of 'opt' to determine input or output.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  complex ,dimension(*) :: sa,sb,sc
  integer,dimension(*) :: isa,isb,isc,jsa,jsb,jsc
  complex  :: alpha,beta
  integer :: opt,N,M
  integer :: i,ka,kb,kc,ja,jb
  logical :: next,row,col,val

  row=.false.
  col=.false.
  val=.false.

  if (opt==0) then
     row=.true.
     col=.true.
     val=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     col=.true.
  elseif (opt==3) then
     val=.true.
  endif

  kc=1
  if (row) isc(1)=kc

  do i=1,N
     ka=isa(i)
     kb=isb(i)

     next=.false.
     do while (.not.next)

        if (ka<=isa(i+1)-1) then
           ja=jsa(ka)
        else
           ja=m+1
        endif

        if (kb<=isb(i+1)-1) then
           jb=jsb(kb)
        else
           jb=m+1
        endif

        if ((ja==jb).and.(ja/=M+1)) then
           if (val) sc(kc)=alpha*sa(ka)+beta*sb(kb)
           if (col) jsc(kc)=ja
           ka=ka+1
           kb=kb+1
           kc=kc+1
        elseif (ja<jb) then
           if (val) sc(kc)=alpha*sa(ka)
           if (col) jsc(kc)=ja
           ka=ka+1
           kc=kc+1
        elseif (jb<ja) then
           if (val) sc(kc)=beta*sb(kb)
           if (col) jsc(kc)=jb
           kb=kb+1
           kc=kc+1
        end if

        if ((ka>isa(i+1)-1).and.(kb>isb(i+1)-1)) next=.true.

     end do

     if (row) isc(i+1)=kc

  end do

end subroutine caddcsr

subroutine csaddcsr(N,M,opt,alpha,sa,isa,jsa,beta,sb,isb,jsb,sc,isc,jsc)
  !  Purpose
  !  =======
  !
  !  Addition of two (NxM any M) sparse CSR matrices A and B and return a new (NxM) CSR matrix C
  !  C=alpha*A+beta*B
  !
  !  A is complex, B and C are real, (alpha,beta are COMPLEX)
  !
  !  Arguments
  !  =========
  !
  !  N,M    (input) INTEGER
  !         The number of rows/columns of the matrices A,B,C.  N,M >= 0.
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs addition, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine computes only values of the array isc of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value ic(n+1) - 1 is the actual
  !         number of the elements in the arrays sc and jsc.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array jsc
  !
  !         if opt=3, the routine has been called previously with the parameter opt=2. On exit it will
  !         return the array sc
  !
  !
  !  alpha  (input) COMPLEX
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa REAL)
  !
  !  beta  (input) COMPLEX
  !
  !  sb,isb,jsb (input) CSR format for the matrix B (sb REAL)
  !
  !  sc,isc,jsc (input/output) CSR format for the matrix C (sc COMPLEX).
  !             see definition of 'opt' to determine input or output.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  real,dimension(*) :: sa,sb
  complex ,dimension(*) :: sc
  integer,dimension(*) :: isa,isb,isc,jsa,jsb,jsc
  complex  :: alpha,beta
  integer :: opt,N,M
  integer :: i,ka,kb,kc,ja,jb
  logical :: next,row,col,val

  row=.false.
  col=.false.
  val=.false.

  if (opt==0) then
     row=.true.
     col=.true.
     val=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     col=.true.
  elseif (opt==3) then
     val=.true.
  endif

  kc=1
  if (row) isc(1)=kc

  do i=1,N
     ka=isa(i)
     kb=isb(i)

     next=.false.
     do while (.not.next)

        if (ka<=isa(i+1)-1) then
           ja=jsa(ka)
        else
           ja=m+1
        endif

        if (kb<=isb(i+1)-1) then
           jb=jsb(kb)
        else
           jb=m+1
        endif

        if ((ja==jb).and.(ja/=M+1)) then
           if (val) sc(kc)=alpha*sa(ka)+beta*sb(kb)
           if (col) jsc(kc)=ja
           ka=ka+1
           kb=kb+1
           kc=kc+1
        elseif (ja<jb)  then
           if (val) sc(kc)=alpha*sa(ka)
           if (col) jsc(kc)=ja
           ka=ka+1
           kc=kc+1
        elseif (jb<ja)  then
           if (val) sc(kc)=beta*sb(kb)
           if (col) jsc(kc)=jb
           kb=kb+1
           kc=kc+1
        end if

        if ((ka>isa(i+1)-1).and.(kb>isb(i+1)-1)) next=.true.

     end do

     if (row) isc(i+1)=kc

  end do

end subroutine csaddcsr

subroutine cinc_addcsr(N,M,alpha,sa,isa,jsa,sb,isb,jsb)
  !  Purpose
  !  =======
  !
  !  Increment matrix B by adding matrix A times alpha
  !  B=B+alpha*A
  !
  !  A,B are complex single precision
  !  We consider (NxM) CSR matrices and that the csr format of B "includes" the csr format of A (matrices should be compatible, i.e. B contains *at least* all the same non-zeros elements than A)
  !
  !  Arguments
  !  =========
  !
  !  N,M      (input) INTEGER
  !         The number of rows/columns of the matrices A,B,C.  N,M >= 0.
  !
  !  alpha  (input) COMPLEX SINGLE PRECISION
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa COMPLEX SINGLE PRECISION)
  !
  !
  !  sb,isb,jsb (input/output) CSR format for the matrix B (sb COMPLEX SINGLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  complex,dimension(*) :: sa,sb
  integer,dimension(*) :: isa,isb,jsa,jsb
  complex :: alpha
  integer :: N,M

  integer :: i,ka,kb

  do i=1,N
     do ka=isa(i),isa(i+1)-1 ! element in row A
        do kb=isb(i),isb(i+1)-1 ! find where it belongs in row B
           if (jsa(ka)==jsb(kb)) then
              sb(kb)=sb(kb)+alpha*sa(ka)
              cycle
           end if
        end do
    enddo
 end do

end subroutine cinc_addcsr

subroutine scsr2csr_up(opt,N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine extract the upper triangular part (including diagonal) of a real CSR matrix A
  !      into a new CSR matrix sA
  !      the size of sa, and sja for csr upper format may be overestimated
  !
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs extraction, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine returns only values of the array sia of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value sia(n+1) - 1 is the actual
  !         number of the elements in the arrays sa and sja.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array sa and sja
  !
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a REAL)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa REAL)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  real,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  real,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k
  logical :: row,colval
!!!!
  row=.false.
  colval=.false.
  if (opt==0) then
     row=.true.
     colval=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     colval=.true.
  end if

  k=0
  if (row) sia(1)=1
  do i=1,n
     if (row) sia(i+1)=sia(i)
     do j=ia(i),ia(i+1)-1
        if (ja(j)>=i) then
           k=k+1
           if (row) sia(i+1)=sia(i+1)+1
           if (colval) then
              sja(k)=ja(j)
              sa(k)=a(j)
           end if
        endif
     end do
  end do
end subroutine scsr2csr_up

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine ccsr2csr_up(opt,N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine extract the upper triangular part (including diagonal) of a real CSR matrix A
  !      into a new CSR matrix sA
  !      the size of sa, and sja for csr upper format may be overestimated
  !
  !
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs extraction, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine returns only values of the array sia of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value sia(n+1) - 1 is the actual
  !         number of the elements in the arrays sa and sja.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array sa and sja
  !
  !
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  complex ,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex ,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k
  logical :: row,colval
!!!!
  row=.false.
  colval=.false.
  if (opt==0) then
     row=.true.
     colval=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     colval=.true.
  end if

  k=0
  if (row) sia(1)=1
  do i=1,n
     if (row) sia(i+1)=sia(i)
     do j=ia(i),ia(i+1)-1
        if (ja(j)>=i) then
           k=k+1
           if (row) sia(i+1)=sia(i+1)+1
           if (colval) then
              sja(k)=ja(j)
              sa(k)=a(j)
           end if
        endif
     end do
  end do
end subroutine ccsr2csr_up

subroutine scsr2csr_low(opt,N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine extract the lower triangular part (including diagonal) of a real CSR matrix A
  !      into a new CSR matrix sA
  !      the size of sa, and sja for csr lower format may be overestimated
  !
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs extraction, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine returns only values of the array sia of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value sia(n+1) - 1 is the actual
  !         number of the elements in the arrays sa and sja.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array sa and sja
  !
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a REAL)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa REAL)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  real,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  real,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k
  logical :: row,colval
!!!!
  row=.false.
  colval=.false.
  if (opt==0) then
     row=.true.
     colval=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     colval=.true.
  end if

  k=0
  if (row) sia(1)=1
  do i=1,n
     if (row) sia(i+1)=sia(i)
     do j=ia(i),ia(i+1)-1
        if (ja(j)<=i) then
           k=k+1
           if (row) sia(i+1)=sia(i+1)+1
           if (colval) then
              sja(k)=ja(j)
              sa(k)=a(j)
           end if
        endif
     end do
  end do
end subroutine scsr2csr_low
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine ccsr2csr_low(opt,N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine extract the lower triangular part (including diagonal) of a COMPLEX CSR matrix A
  !      into a new CSR matrix sA
  !      the size of sa, and sja for csr lower format may be overestimated
  !
  !
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs extraction, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine returns only values of the array sia of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value sia(n+1) - 1 is the actual
  !         number of the elements in the arrays sa and sja.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array sa and sja
  !
  !
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  complex ,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex ,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k
  logical :: row,colval
!!!!
  row=.false.
  colval=.false.
  if (opt==0) then
     row=.true.
     colval=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     colval=.true.
  end if

  k=0
  if (row) sia(1)=1
  do i=1,n
     if (row) sia(i+1)=sia(i)
     do j=ia(i),ia(i+1)-1
        if (ja(j)<=i) then
           k=k+1
           if (row) sia(i+1)=sia(i+1)+1
           if (colval) then
              sja(k)=ja(j)
              sa(k)=a(j)
           end if
        endif
     end do
  end do
end subroutine ccsr2csr_low

SUBROUTINE scsrmm(UPLO,TRANS,N,M,rhs,alpha,a,ia,ja,X,beta,B)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine peforms
  !
  !      b=alpha*o(A)*x+beta*b
  !
  !      where the NxM matrix A is stored in CSR format (a,ia,ja),
  !      and X is Mxrhs and B is Nxrhs
  !
  !      REAL SINGLE PRECISION version
  !
  !      Arguments:
  !
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !       if UPLO='U' only the upper part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !
  !  TRANS (input)  CHARACTER*1
  !       if TRANS='N' o(A)=A -- normal mode
  !       if TRANS='T' o(A)=A^T -- transpose mode
  !
  !  N    (input) INTEGER
  !        The number of row of the matrix A and row of matrix B.  N >= 0.
  !  M    (input) INTEGER
  !        The number of column of matrix A and row of matrix X. M>=0; M=N (square matrix) for UPLO=L,U
  !  rhs  (input) INTEGER
  !        The number of column of matrices B and C
  !
  !  alpha (input) SINGLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a SINGLE PRECISION)
  !
  !  X     (input) DOUBLE PRECISION
  !        matrix of size (Mxrhs) or (Nxrhs) if TRANS='T'
  !
  !  beta  (input) SINGLE PRECISION
  !
  !  B     (input/output) SINGLE PRECISION
  !        matrix of size (Nxrhs) or (Mxrhs) if TRANS='T'. On exit contains the solution.
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  character(len=1) :: UPLO,TRANS
  integer :: N,M,rhs
  real :: alpha,beta
  real,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  real,dimension(*):: x
  real,dimension(*) ::b
!!!!!!!!

  integer ::i,j,k,s,n1
  real,parameter :: SZERO=0.0e0

  if (TRANS=='N') then
     n1=N
  elseif (TRANS=='T') then
     n1=M
  endif

!!!!!!!! Initialization

  if (beta/=SZERO) then
     call SSCAL(n1*rhs,beta,b,1)
  else
     b(1:n1*rhs)=SZERO
  end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (UPLO=='F') then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if (TRANS=='N') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(i+(s-1)*N)=b(i+(s-1)*N)+alpha*a(k)*x(j+(s-1)*M)
              enddo
           end do
        end do

     elseif (TRANS=='T') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(j+(s-1)*M)=b(j+(s-1)*M)+alpha*a(k)*x(i+(s-1)*N)
              enddo
           end do
        end do

     end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  elseif ((UPLO=='L').or.(UPLO=='U')) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     do i=1,N
        do k=ia(i),ia(i+1)-1
           j=ja(k)
           do s=1,rhs
              b(i+(s-1)*N)=b(i+(s-1)*N)+alpha*a(k)*x(j+(s-1)*M)
              if (j/=i) b(j+(s-1)*M)=b(j+(s-1)*M)+alpha*a(k)*x(i+(s-1)*N)
           enddo
        end do
     end do

  end if

end SUBROUTINE scsrmm

SUBROUTINE ccsrmm(UPLO,TRANS,N,M,rhs,alpha,a,ia,ja,X,beta,B)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine peforms
  !
  !      b=alpha*o(A)*x+beta*b
  !
  !      where the NxM matrix A is stored in CSR format (a,ia,ja),
  !      and X is Mxrhs and B is Nxrhs
  !
  !      COMPLEX SINGLE PRECISION version
  !
  !
  !      Arguments:
  !
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !       if UPLO='U' only the upper part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !
  !  TRANS (input)  CHARACTER*1
  !       if TRANS='N' o(A)=A   -- normal mode
  !       if TRANS='T' o(A)=A^T -- transpose mode
  !       if TRANS='C' o(A)=A^H -- transpose conjugate mode
  !
  !  N    (input) INTEGER
  !        The number of row of the matrix A and row of matrix B.  N >= 0.
  !  M    (input) INTEGER
  !        The number of column of matrix A and row of matrix X. M>=0; M=N (square matrix) for UPLO=L,U
  !  rhs  (input) INTEGER
  !        The number of column of matrices B and C
  !
  !  alpha (input) SINGLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a SINGLE PRECISION)
  !
  !  X     (input) SINGLE PRECISION
  !        matrix of size (Mxrhs) or (Nxrhs) if TRANS='T','C'.
  !
  !  beta  (input) SINGLE PRECISION
  !
  !  B     (input/output) SINGLE PRECISION
  !        matrix of size (Nxrhs) or (Mxrhs) if TRANS='T','C'. On exit contains the solution.
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  character(len=1) :: UPLO,TRANS
  integer :: N,M,rhs
  complex  :: alpha,beta
  complex ,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex ,dimension(*):: x
  complex ,dimension(*) ::b
!!!!!!!!

  integer ::i,j,k,s,n1
  complex,parameter :: ZEROC=(0.0e0,0.0e0)

  if (TRANS=='N') then
     n1=N
  else
     n1=M
  endif

!!!!!!!! Initialization

  if (beta/=ZEROC) then
     call CSCAL(n1*rhs,beta,b,1)
  else
     b(1:n1*rhs)=ZEROC
  end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (UPLO=='F') then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if (TRANS=='N') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(i+(s-1)*N)=b(i+(s-1)*N)+alpha*a(k)*x(j+(s-1)*M)
              end do
           end do
        end do

     elseif (TRANS=='T') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(j+(s-1)*M)=b(j+(s-1)*M)+alpha*a(k)*x(i+(s-1)*N)
              enddo
           end do
        end do

     elseif (TRANS=='C') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(j+(s-1)*M)=b(j+(s-1)*M)+alpha*conjg(a(k))*x(i+(s-1)*N)
              end do
           end do
        end do

     end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  elseif ((UPLO=='L').or.(UPLO=='U')) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if ((TRANS=='N').or.(TRANS=='T')) then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(i+(s-1)*N)=b(i+(s-1)*N)+alpha*a(k)*x(j+(s-1)*M)
                 if (j/=i) b(j+(s-1)*M)=b(j+(s-1)*M)+alpha*a(k)*x(i+(s-1)*N)
              end do
           end do
        end do

     elseif (TRANS=='C') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(j+(s-1)*M)=b(j+(s-1)*M)+alpha*conjg(a(k))*x(i+(s-1)*N)
                 if (j/=i) b(i+(s-1)*N)=b(i+(s-1)*N)+alpha*conjg(a(k))*x(j+(s-1)*M)
              end do
           end do
        end do

     end if

  end if

END SUBROUTINE ccsrmm

SUBROUTINE chcsrmm(UPLO,TRANS,N,M,rhs,alpha,a,ia,ja,X,beta,B)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine peforms
  !
  !      b=alpha*o(A)*x+beta*b
  !
  !      where the NxM matrix A is stored in CSR format (a,ia,ja),
  !      and X is Mxrhs and B is Nxrhs
  !
  !      Note: A is complex Hermitian if UPLO='L or UPLO='U'
  !
  !      Arguments:
  !
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided- symmetric/hermitian CSR (N must be equal to M)
  !       if UPLO='U' only the upper part of the matrix A is provided- symmetric/hermitian CSR (N must be equal to M)
  !
  !  TRANS (input)  CHARACTER*1
  !       if TRANS='N' o(A)=A   -- normal mode
  !       if TRANS='T' o(A)=A^T -- transpose mode
  !       if TRANS='C' o(A)=A^H -- transpose conjugate mode
  !
  !  N    (input) INTEGER
  !        The number of row of the matrix A and row of matrix B.  N >= 0.
  !  M    (input) INTEGER
  !        The number of column of matrix A and row of matrix X. M>=0; M=N (square matrix) for UPLO=L,U
  !  rhs  (input) INTEGER
  !        The number of column of matrices B and C
  !
  !  alpha (input) COMPLEX SINGLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX SINGLE PRECISION)
  !
  !  X     (input) COMPLEX SINGLE PRECISION
  !        matrix of size (Mxrhs) or (Nxrhs) if TRANS='T.'
  !
  !  beta  (input) COMPLEX SINGLE PRECISION
  !
  !  B     (input/output) COMPLEX SINGLE PRECISION
  !        matrix of size (Nxrhs) or (Mxrhs) if TRANS='T','C'. On exit contains the solution.
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  character(len=1) :: UPLO,TRANS
  integer :: N,M,rhs
  complex :: alpha,beta
  complex,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex,dimension(*) :: x
  complex,dimension(*) :: b
!!!!!!!!

  integer ::i,j,k,s,n1
  real,parameter :: SZERO=0.0e0
  complex,parameter :: ZEROC=(SZERO,SZERO)

  if (TRANS=='N') then
     n1=N
  else
     n1=M
  endif

!!!!!!!! Initialization

  if (beta/=ZEROC) then
     call CSCAL(n1*rhs,beta,b,1)
  else
     b(1:n1*rhs)=ZEROC
  end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (UPLO=='F') then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if (TRANS=='N') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(i+(s-1)*N)=b(i+(s-1)*N)+alpha*a(k)*x(j+(s-1)*M)
              end do
           end do
        end do

     elseif (TRANS=='T') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(j+(s-1)*M)=b(j+(s-1)*M)+alpha*a(k)*x(i+(s-1)*N)
              end do
           end do
        end do

     elseif (TRANS=='C') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(j+(s-1)*M)=b(j+(s-1)*M)+alpha*conjg(a(k))*x(i+(s-1)*N)
              end do
           end do
        end do

     end if

  elseif ((UPLO=='L').or.(UPLO=='U')) then

     if ((TRANS=='N').or.(TRANS=='C')) then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(i+(s-1)*N)=b(i+(s-1)*N)+alpha*a(k)*x(j+(s-1)*M)
                 if (j/=i) b(j+(s-1)*M)=b(j+(s-1)*M)+alpha*conjg(a(k))*x(i+(s-1)*N)
              end do
           end do
        end do

     elseif (TRANS=='T') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(j+(s-1)*M)=b(j+(s-1)*M)+alpha*conjg(a(k))*x(i+(s-1)*N)
                 if (j/=i) b(i+(s-1)*N)=b(i+(s-1)*N)+alpha*a(k)*x(j+(s-1)*M)
              end do
           end do
        end do

     end if

  end if

end SUBROUTINE chcsrmm

SUBROUTINE sccsrmm(UPLO,TRANS,N,M,rhs,alpha,a,ia,ja,X,beta,B)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine peforms
  !
  !      b=alpha*o(A)*x+beta*b
  !
  !      where the NxM matrix A is stored in CSR format (a,ia,ja),
  !      and X is Mxrhs and B is Nxrhs
  !
  !      SINGLE PRECISION version
  !
  !      Note: same as ccsrmm but A is real
  !
  !      Arguments:
  !
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !       if UPLO='U' only the upper part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !
  !  TRANS (input)  CHARACTER*1
  !       if TRANS='N' o(A)=A   -- normal mode
  !       if TRANS='T' o(A)=A^T -- transpose mode
  !       if TRANS='C' o(A)=A^H -- transpose conjugate mode
  !
  !  N    (input) INTEGER
  !        The number of row of the matrix A and row of matrix B.  N >= 0.
  !  M    (input) INTEGER
  !        The number of column of matrix A and row of matrix X. M>=0; M=N (square matrix) for UPLO=L,U
  !  rhs  (input) INTEGER
  !        The number of column of matrices B and C
  !
  !  alpha (input) COMPLEX SINGLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a SINGLE PRECISION)
  !
  !  X     (input) COMPLEX SINGLE PRECISION
  !        matrix of size (Mxrhs) or (Nxrhs) if TRANS='T','C'.
  !
  !  beta  (input) COMPLEX SINGLE PRECISION
  !
  !  B     (input/output) COMPLEX SINGLE PRECISION
  !        matrix of size (Nxrhs) or (Mxrhs) if TRANS='T','C'.. On exit contains the solution.
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  character(len=1) :: UPLO,TRANS
  integer :: N,M,rhs
  complex :: alpha,beta
  real,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex,dimension(*):: x
  complex,dimension(*) ::b
!!!!!!!!

  integer ::i,j,k,s,n1
  complex,parameter :: ZEROC=(0.0e0,0.0e0)

  if (TRANS=='N') then
     n1=N
  else
     n1=M
  endif

!!!!!!!! Initialization

  if (beta/=ZEROC) then
     call CSCAL(n1*rhs,beta,b,1)
  else
     b(1:n1*rhs)=ZEROC
  end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (UPLO=='F') then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if (TRANS=='N') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(i+(s-1)*N)=b(i+(s-1)*N)+alpha*a(k)*x(j+(s-1)*M)
              end do
           end do
        end do

     elseif ((TRANS=='T').or.(TRANS=='C')) then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              do s=1,rhs
                 b(j+(s-1)*M)=b(j+(s-1)*M)+alpha*a(k)*x(i+(s-1)*N)
              end do
           end do
        end do

     end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  elseif ((UPLO=='L').or.(UPLO=='U')) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     do i=1,N
        do k=ia(i),ia(i+1)-1
           j=ja(k)
           do s=1,rhs
              b(i+(s-1)*N)=b(i+(s-1)*N)+alpha*a(k)*x(j+(s-1)*M)
              if (j/=i) b(j+(s-1)*M)=b(j+(s-1)*M)+alpha*a(k)*x(i+(s-1)*N)
           end do
        end do
     end do

  end if

END SUBROUTINE sccsrmm

subroutine scsr_transpose(N,M,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine transpose a double real CSR matrix A
  !      into a new CSR matrix sA
  !
  !  N,M    (input) INTEGER
  !          The number of rows/columns of the matrices A.  N,M >= 0.
  !          The number of columns/rows of the matrices sA. N,M >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a SINGLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa SINGLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: N,M
  real,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  real,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
  integer :: i,j,k,istart
  integer,dimension(:),allocatable :: tab
!!!!

  allocate(tab(M))
!!!!
!!!!

  !! #of elements by row
  tab=0
  istart=N
  do i=1,N
     do j=ia(i),ia(i+1)-1
        tab(ja(j))=tab(ja(j))+1
        if (ja(j)<istart) istart=ja(j)
     end do
  end do

  !! find sia
  do i=1,istart-1
     sia(i)=1 !!<<0
  end do
  sia(istart)=1
  do i=istart+1,M+1
     sia(i)=sia(i-1)+tab(i-1)
  end do

  tab=sia(1:M)
  !! find jsa,sa
  do i=1,N
     do k=ia(i),ia(i+1)-1
        sja(tab(ja(k)))=i
        sa(tab(ja(k)))=a(k)
        tab(ja(k))=tab(ja(k))+1
     end do
  end do

  deallocate(tab)

end subroutine scsr_transpose

subroutine ccsr_transpose(N,M,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine performs the transpose  of
  !      a double complex CSR matrix A into a new CSR matrix sA
  !
  !
  !  N,M    (input) INTEGER
  !          The number of rows/columns of the matrices A.  N,M >= 0.
  !          The number of columns/rows of the matrices sA. N,M >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX SINGLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX SINGLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: N,M
  complex,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
  integer :: i,j,k,istart
  integer,dimension(:),allocatable :: tab
!!!!

  allocate(tab(M))
!!!!
!!!!

  !! #of elements by row
  tab=0
  istart=N
  do i=1,N
     do j=ia(i),ia(i+1)-1
        tab(ja(j))=tab(ja(j))+1
        if (ja(j)<istart) istart=ja(j)
     end do
  end do

  !! find sia
  do i=1,istart-1
     sia(i)=1 !!<<0
  end do
  sia(istart)=1
  do i=istart+1,M+1
     sia(i)=sia(i-1)+tab(i-1)
  end do

  tab=sia(1:M)
  !! find jsa,sa
  do i=1,N
     do k=ia(i),ia(i+1)-1
        sja(tab(ja(k)))=i
        sa(tab(ja(k)))=a(k)
        tab(ja(k))=tab(ja(k))+1
     end do
  end do

  deallocate(tab)

end subroutine ccsr_transpose

subroutine ccsr_htranspose(N,M,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine performs the transpose conjugate of
  !      a double complex CSR matrix A into a new CSR matrix sA
  !
  !
  !  N,M    (input) INTEGER
  !          The number of rows/columns of the matrices A.  N,M >= 0.
  !          The number of columns/rows of the matrices sA. N,M >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX SINGLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX SINGLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: N,M
  complex,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
  integer :: i,j,k,istart
  integer,dimension(:),allocatable :: tab
!!!!

  allocate(tab(M))
!!!!
!!!!

  !! #of elements by row
  tab=0
  istart=N
  do i=1,N
     do j=ia(i),ia(i+1)-1
        tab(ja(j))=tab(ja(j))+1
        if (ja(j)<istart) istart=ja(j)
     end do
  end do

  !! find sia
  do i=1,istart-1
     sia(i)=1 !!<< integer,dimension(:),allocatable :: tab
!!!!

     allocate(tab(M))
!!!!
  end do
  sia(istart)=1
  do i=istart+1,M+1
     sia(i)=sia(i-1)+tab(i-1)
  end do

  tab=sia(1:M)
  !! find jsa,sa
  do i=1,N
     do k=ia(i),ia(i+1)-1
        sja(tab(ja(k)))=i
        sa(tab(ja(k)))=conjg(a(k))
        tab(ja(k))=tab(ja(k))+1
     end do
  end do

  deallocate(tab)

end subroutine ccsr_htranspose

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine scsr_uplo_to_csr(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a real upper or lower CSR matrix A
  !      into a new full CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a REAL)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa REAL)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  !  include 'f90_noruntime_interface.fi'
  integer :: N
  real,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  real,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,k,nnza,info,opt
  real,dimension(:),allocatable ::sb
  integer,dimension(:),allocatable :: sib
  integer,dimension(:),allocatable :: sjb
  real,parameter :: SONE=1.0e0,SZERO=0.0e0
!!!!
  nnza=ia(n+1)-1

  allocate(sb(nnza))
  allocate(sib(N+1))
  allocate(sjb(nnza))

!!! transpose
  call scsr_transpose(N,N,a,ia,ja,sb,sib,sjb)
!!! nullify diagonal
  do i=1,N
     do k=sib(i),sib(i+1)-1
        if (sjb(k)==i) sb(k)=SZERO
     end do
  enddo
  !! add the two matrices to form the full matrix
  opt=0
  call saddcsr(N,N,opt,SONE,a,ia,ja,SONE,sb,sib,sjb,sa,sia,sja)

  deallocate(sb)
  deallocate(sib)
  deallocate(sjb)

end subroutine scsr_uplo_to_csr

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine ccsr_uplo_to_csr(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a complex symetric upper or lower CSR matrix A
  !      into a new full CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  !  include 'f90_noruntime_interface.fi'
  integer :: N
  complex,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,k,nnza,info,opt
  complex,dimension(:),allocatable ::sb
  integer,dimension(:),allocatable :: sib
  integer,dimension(:),allocatable :: sjb
  real,parameter :: SONE=1.0e0,SZERO=0.0e0
  complex,parameter :: ONEC=(SONE,SZERO),ZEROC=(SZERO,SZERO)
!!!!
  nnza=ia(n+1)-1

  allocate(sb(nnza))
  allocate(sib(N+1))
  allocate(sjb(nnza))

!!! transpose
  call ccsr_transpose(N,N,a,ia,ja,sb,sib,sjb)
!!! nullify diagonal
  do i=1,N
     do k=sib(i),sib(i+1)-1
        if (sjb(k)==i) sb(k)=ZEROC
     end do
  enddo
  !! add the two matrices to form the full matrix
  opt=0
  call caddcsr(N,N,opt,ONEC,a,ia,ja,ONEC,sb,sib,sjb,sa,sia,sja)

  deallocate(sb)
  deallocate(sib)
  deallocate(sjb)

end subroutine ccsr_uplo_to_csr

subroutine chcsr_uplo_to_csr(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a complex Hermitian upper or lower CSR matrix A
  !      into a new full CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  !  include 'f90_noruntime_interface.fi'
  integer :: N
  complex ,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex ,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,k,nnza,info,opt
  complex ,dimension(:),allocatable ::sb
  integer,dimension(:),allocatable :: sib
  integer,dimension(:),allocatable :: sjb
  real,parameter :: SONE=1.0e0,SZERO=0.0e0
  complex ,parameter :: ONEC=(SONE,SZERO),ZEROC=(SZERO,SZERO)
!!!!
  nnza=ia(n+1)-1

  allocate(sb(nnza))
  allocate(sib(N+1))
  allocate(sjb(nnza))

!!! transpose
  call ccsr_htranspose(N,N,a,ia,ja,sb,sib,sjb)
!!! nullify diagonal
  do i=1,N
     do k=sib(i),sib(i+1)-1
        if (sjb(k)==i) sb(k)=ZEROC
     end do
  enddo
  !! add the two matrices to form the full matrix
  opt=0
  call caddcsr(N,N,opt,ONEC,a,ia,ja,ONEC,sb,sib,sjb,sa,sia,sja)

  deallocate(sb)
  deallocate(sib)
  deallocate(sjb)

end subroutine chcsr_uplo_to_csr

subroutine scsr_convert_upper(n,UPLO,sa,isa,jsa,ssa,sisa,sjsa)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a double real upper, lower or full CSR matrix A
  !      into a new CSR upper matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided
  !       if UPLO='U' only the upper part of the matrix A is provided
  !
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (a SINGLE PRECISION)
  !
  !  ssa,sisa,sjsa (output) CSR format for the matrix sA (sa SINGLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: n
  character(len=1) :: UPLO
  real,dimension(*) :: sa,ssa
  integer,dimension(*) :: isa,jsa,sisa,sjsa
  integer :: opt,nnz

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  if ((UPLO=='F').or.(UPLO=='f')) then !!! full csr to upper-csr

     opt=0
     call scsr2csr_up(opt,N,sa,isa,jsa,ssa,sisa,sjsa)

!!!!!!!!!!!!!!!!!!!!!!!!!
  elseif ((UPLO=='U').or.(UPLO=='u')) then !! upper-csr already
     nnz=isa(n+1)-1
     call SCOPY(nnz,sa,1,ssa,1)
     sisa(1:n+1)= isa(1:n+1)
     sjsa(1:nnz)= jsa(1:nnz)

  elseif ((UPLO=='L').or.(UPLO=='l')) then !!! lower-csr to upper-csr'
     call scsr_transpose(N,N,sa,isa,jsa,ssa,sisa,sjsa)
  endif

end subroutine scsr_convert_upper

subroutine chcsr_convert_full(n,UPLO,sa,isa,jsa,ssa,sisa,sjsa)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a single complex upper, lower or full Hermitian CSR matrix A
  !      into a new full CSR  matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided
  !       if UPLO='U' only the upper part of the matrix A is provided
  !
  !
  !  sa,isa,jsa (input) CSR format for the Hermitian matrix A (a COMPLEX SINGLE PRECISION)
  !
  !  ssa,sisa,sjsa (output) CSR format for the Hermitian matrix sA (sa COMPLEX SINGLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: n
  character(len=1) :: UPLO
  complex,dimension(*) :: sa,ssa
  integer,dimension(*) :: isa,jsa,sisa,sjsa
  integer :: opt,nnz

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if ((UPLO=='F').or.(UPLO=='f')) then !!! direct copy

     nnz=isa(n+1)-1
     call CCOPY(nnz,sa,1,ssa,1)
     sisa(1:n+1)= isa(1:n+1)
     sjsa(1:nnz)= jsa(1:nnz)

!!!!!!!!!!!!!!!!!!!!!!!!!
  else!! upper-csr or lower-csr to full csr

     call chcsr_uplo_to_csr(N,sa,isa,jsa,ssa,sisa,sjsa)
  end if

end subroutine chcsr_convert_full

subroutine scoo2csr(n,nnz,ic,jc,c,isa,jsa,sa)
  implicit none
  integer :: n,nnz
  integer,dimension(*) :: ic,jc,isa,jsa
  real,dimension(*) :: c,sa
!!!
  integer :: k,k1,i,j,idum
  integer,dimension(n) :: iloc
  real :: adum

  isa(1:N+1) = 0
  !find how many elements in row
  do  k=1, nnz
     isa(ic(k)) = isa(ic(k))+1
  end do

  !build isa
  k = 1
  do  i=1,N+1
     k1 = isa(i)
     isa(i) = k
     k = k+k1
  end do

  iloc(1:n)=isa(1:n)
  !Build jsa, sa - increment local row counter
  do  k=1, nnz
     sa(iloc(ic(k))) =  c(k)
     jsa(iloc(ic(k))) = jc(k)
     iloc(ic(k)) = iloc(ic(k))+1
  end do
  ! Reorder by increasing column
  do i=1,n
     do k=isa(i),isa(i+1)-1
        do k1=k,isa(i+1)-1
           if (jsa(k1)<jsa(k)) then
              idum=jsa(k)
              jsa(k)=jsa(k1)
              jsa(k1)=idum
              adum=sa(k)
              sa(k)=sa(k1)
              sa(k1)=adum
           endif
        enddo
     enddo
  enddo

end subroutine scoo2csr

subroutine ccoo2csr(n,nnz,ic,jc,c,isa,jsa,sa)
  implicit none
  integer :: n,nnz
  integer,dimension(*) :: ic,jc,isa,jsa
  complex,dimension(*) :: c,sa
!!!
  integer :: k,k1,i,j,idum
  integer,dimension(n) :: iloc
  complex:: adum

  isa(1:N+1) = 0
  !find how many elements in row
  do  k=1, nnz
     isa(ic(k)) = isa(ic(k))+1
  end do

  !build isa
  k = 1
  do  i=1,N+1
     k1 = isa(i)
     isa(i) = k
     k = k+k1
  end do

  iloc(1:n)=isa(1:n)
  !Build jsa, sa - increment local row counter
  do  k=1, nnz
     sa(iloc(ic(k))) =  c(k)
     jsa(iloc(ic(k))) = jc(k)
     iloc(ic(k)) = iloc(ic(k))+1
  end do
  ! Reorder by increasing column
  do i=1,n
     do k=isa(i),isa(i+1)-1
        do k1=k,isa(i+1)-1
           if (jsa(k1)<jsa(k)) then
              idum=jsa(k)
              jsa(k)=jsa(k1)
              jsa(k1)=idum
              adum=sa(k)
              sa(k)=sa(k1)
              sa(k1)=adum
           endif
        enddo
     enddo
  enddo

end subroutine ccoo2csr

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine ssubmat (n,i1,i2,j1,j2,a,ja,ia,ao,jao,iao)
  !*******************************************************************************
  ! SUBMAT extracts the submatrix A(i1:i2,j1:j2)
  !*******************************************************************************
  !a,ja,ia = original matrix in CSR format
  !ao,jao,iao = extracted matrix in CSR format
  implicit none
  integer :: i1,i2,j1,j2,n
  real,dimension(*):: a
  integer,dimension(*) :: ia,ja
  real,dimension(*) :: ao
  integer,dimension(*) :: iao,jao

  integer :: klen,i,ii,k1,k2,k,j,nr,nc

  nr = i2-i1+1
  nc = j2-j1+1
  klen = 0
  !!
  do  i = 1,nr
     ii = i1+i-1
     k1 = ia(ii)
     k2 = ia(ii+1)-1
     iao(i) = klen+1
     do  k=k1,k2
        j = ja(k)
        if ((j .ge. j1).and. (j .le. j2)) then
           klen = klen+1
           ao(klen) = a(k)
           jao(klen) = j - j1+1
        endif
     end do
  end do
  iao(nr+1) = klen+1
end subroutine ssubmat
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine csubmat (n,i1,i2,j1,j2,a,ja,ia,ao,jao,iao)
  !*******************************************************************************
  ! SUBMAT extracts the submatrix A(i1:i2,j1:j2)
  !*******************************************************************************
  !a,ja,ia = original matrix in CSR format
  !ao,jao,iao = extracted matrix in CSR format
  implicit none
  integer :: i1,i2,j1,j2,n
  complex,dimension(*):: a
  integer,dimension(*) :: ia,ja
  complex,dimension(*) :: ao
  integer,dimension(*) :: iao,jao

  integer :: klen,i,ii,k1,k2,k,j,nr,nc

  nr = i2-i1+1
  nc = j2-j1+1
  klen = 0
  !!
  do  i = 1,nr
     ii = i1+i-1
     k1 = ia(ii)
     k2 = ia(ii+1)-1
     iao(i) = klen+1
     do  k=k1,k2
        j = ja(k)
        if ((j .ge. j1).and. (j .le. j2)) then
           klen = klen+1
           ao(klen) = a(k)
           jao(klen) = j - j1+1
        endif
     end do
  end do
  iao(nr+1) = klen+1
end subroutine csubmat
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

max_stars_repo_name: lkusch/Kratos

max_issues_repo_name: lkusch/Kratos

max_forks_repo_name: lkusch/Kratos


text:
      SUBROUTINE ARCELL(XSTART,YSTART,ANGLE,ECCEN)
C
C          ------------------------------------------------
C          ROUTINE NO. ( 103)   VERSION (A8.6)    20:JUN:86
C          ------------------------------------------------
C
C          THIS DRAWS THE ARC OF AN ELLIPSE WHICH IS
C          CENTRED ON THE CURRENT PLOTTING POSITION.
C
C
C          <XSTART> IS THE STARTING X-COORDINATE,
C          <YSTART> IS THE STARTING Y-COORDINATE,
C          <ANGLE>  IS THE ARC ANGLE ANTI-CLOCKWISE
C                   FROM THE STARTING POSITION,
C          <ECCEN>  IS THE ECCENTRICITY OF THE ELLIPSE.
C
C
      REAL    RDATA(4)
      INTEGER IDATA(1)
C
      COMMON /T0ACON/ ANGCON
      COMMON /T0MAPP/ X1MAP0,X2MAP0,Y1MAP0,Y2MAP0
      COMMON /T0PPOS/ XPLOT0,YPLOT0
      COMMON /T0TRAC/ IPRINT
      COMMON /T0TRAR/ RTRAC1,RTRAC2,RTRAC3,RTRAC4
      COMMON /T3CONS/ PI
      COMMON /T3NBYR/ NBYTR
C
      DATA IDATA /0/
C
C
      RTRAC1= XSTART
      RTRAC2= YSTART
      RTRAC3= ANGLE
      RTRAC4= ECCEN
      IF (IPRINT.EQ.1) CALL G0MESG(45,4)
C
      RDATA(1)= XSTART-XPLOT0
      RDATA(2)= YSTART-YPLOT0
      RDATA(4)= ABS(ECCEN)
C
C          CALCULATE THE TRUE ROTATION ANGLE ALLOWING FOR
C          DIFFERENT SCALING ON THE X AND Y AXES.
C
      YCOORD= RDATA(2)*(X2MAP0-X1MAP0)/(Y2MAP0-Y1MAP0)
      ANG= ANGLE*ANGCON
      SINANG= SIN(ANG)
      COSANG= COS(ANG)
      ANUM= RDATA(4)*SINANG*(RDATA(1)*RDATA(1)+YCOORD*YCOORD)
      DENOM= COSANG*(RDATA(1)*RDATA(1)+
     &       YCOORD*YCOORD*RDATA(4)*RDATA(4))+
     &       SINANG*RDATA(1)*YCOORD*(RDATA(4)*RDATA(4)-1)
      RDATA(3)= ATAN2(ANUM,DENOM)
      IF ((RDATA(3)*ANG).LT.0.0)
     &     RDATA(3)= RDATA(3)+PI*SIGN(1.0,ANG)
    1   IF (ABS(RDATA(3)-ANG).LT.PI*0.5) GO TO 2
        RDATA(3)= RDATA(3)+PI*SIGN(1.0,ANG)
       GO TO 1
C
    2 CALL G3LINK(0,9,4*NBYTR,IDATA,RDATA)
C
      RETURN
      END

ext: f

max_stars_repo_name: ZedThree/GHOST

max_issues_repo_name: ZedThree/GHOST

max_forks_repo_name: ZedThree/GHOST


text:
      subroutine chkinz(ier,imchmx,imech,iopt,jcode,kmax,kxt,nelect,
     $ noptmx,noutpt,no2gaq,nrct,nrctmx,nrk,nstmax,nttyo,rkb,ureac,
     $ uspeca,uzveci,zvclgi)
c
c     This subroutine checks the code input for various kinds of
c     errors and inconsistencies. Here ier accumulates the
c     number of errors caught.
c
c     This subroutine is called by:
c
c       EQ6/eq6.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c
c     Principal output:
c
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer imchmx,kmax,noptmx,nrctmx,nstmax
c
      integer noutpt,nttyo
c
      integer imech(2,nrctmx),iopt(noptmx),jcode(nrctmx),nrk(2,nrctmx)
c
      integer ier,kxt,nelect,no2gaq,nrct
c
      character*48 uspeca(nstmax),uzveci(kmax)
      character*24 ureac(nrctmx)
c
      real*8 rkb(imchmx,2,nrctmx),zvclgi(kmax)
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer i,jlen,j2,kcol,nrc
c
      integer ilnobl
c
      character*56 uspn56
c
c-----------------------------------------------------------------------
c
      ier = 0
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Rate law data.
c
c       nrk(1,nrc)= forward rate law code:
c
c         -1 = Use the backward rate law form (legal only if
c              nrk(2,nrc) = 2)
c          0 = Illegal value
c          1 = Relative rate
c          2 = Transition state theory net rate
c          3 = Linear rate
c
c       For the case nrk(1,nrc) = 2, the reactant must be
c       either a pure mineral (jcode = 0) or a solid solution
c       (jcode = 1).
c
c       nrk(2,nrc)= backward rate law code:
c
c         -1 = Use the forward rate law form (legal only if
c              nrk(1,nrc) = 2)
c          0 = No rate law specified; the reaction may be controlled
c              by partial equilibrium.
c          1 = Relative rate
c          2 = Transition state theory net rate
c          3 = Linear rate
c
c       For the case nrk(2,nrc) = 2, the reactant must be
c       either a pure mineral (jcode = 0) or a solid solution
c       (jcode = 1)
c
      do nrc = 1,nrct
        j2 = ilnobl(ureac(nrc))
c
        if (iopt(2) .le. 0) then
c
c         Check for rate law choices that can't be used without
c         an explicit time frame.
c
          if (nrk(2,nrc) .ge. 2) then
            write (noutpt,1000) ureac(nrc)(1:j2),nrk(2,nrc)
            write (nttyo,1000) ureac(nrc)(1:j2),nrk(2,nrc)
 1000       format(/' * Error - (EQ6/chkinz) The forward rate law code',
     $      /7x,'for ',a,' is ',i2,', which is not valid',
     $      /7x,'unless iopt(2) = 1 (the model has a time frame).')
            ier = ier + 1
          endif
c
          if (nrk(2,nrc) .ge. 2) then
            write (noutpt,1010) ureac(nrc)(1:j2),nrk(2,nrc)
            write (nttyo,1010) ureac(nrc)(1:j2),nrk(2,nrc)
 1010       format(/' * Error - (EQ6/chkinz) The backward rate law',
     $      ' code',/7x,'for ',a,' is ',i2,', which is not valid',
     $      /7x,'unless iopt(2) = 1 (the model has a time frame).')
            ier = ier + 1
          endif
        endif
c
c       Check for invalid forward and backward rate law combinations.
c       The following combinations are invalid:
c
c         nrk(1,nrc) = 1 and nrk(2,nrc) = 1
c         nrk(1,nrc) = 3 and nrk(2,nrc) = 3
c         nrk(1,nrc) = -1 and nrk(2,nrc) is not 2
c         nrk(2,nrc) = 0
c
c
        if (nrk(1,nrc).eq.1 .and. nrk(2,nrc).eq.1) then
          write (noutpt,1020) ureac(nrc)(1:j2)
          write (nttyo,1020) ureac(nrc)(1:j2)
 1020     format(' * Error - (EQ6/chkinz) The forward and backward',
     $    /7x,'rate law codes for ',a,' are both 1, which is not a',
     $    /7x,'valid combination.')
          ier = ier + 1
        endif
c
        if (nrk(1,nrc).eq.3 .and. nrk(2,nrc).eq.3) then
          write (noutpt, 1030) ureac(nrc)(1:j2)
          write (nttyo, 1030) ureac(nrc)(1:j2)
 1030     format(' * Error - (EQ6/chkinz) The forward and backward',
     $    /7x,'rate law codes for ',a,' are both 3, which is not a',
     $    /7x,'valid combination.')
          ier = ier + 1
        endif
c
        if (nrk(1,nrc) .eq. -1) then
          if (nrk(2,nrc) .ne. 2) then
            write (noutpt, 1040) ureac(nrc)(1:j2)
            write (nttyo, 1040) ureac(nrc)(1:j2)
 1040       format(' * Error - (EQ6/chkinz) The forward rate law code',
     $      /7x,'for ',a,' is -1, which is not a valid value  unless',
     $      /7x,'the backward rate law code is 2.')
            ier = ier + 1
          endif
        endif
c
        if (nrk(1,nrc) .eq. 0) then
          write (noutpt,1050) ureac(nrc)(1:j2)
          write (nttyo,1050) ureac(nrc)(1:j2)
 1050     format(' * Error - (EQ6/chkinz) The forward rate law code',
     $    /7x,'for ',a,' is 0, which is not a valid value.')
          ier = ier + 1
        endif
c
c       Check the parameters for the forward rate laws.
c
        if (nrk(1,nrc) .eq. 2) then
c
c         Transition state theory.
c
          do i = 1,imech(1,nrc)
            if (rkb(i,1,nrc) .lt. 0.) then
              write (noutpt,1060) i,rkb(i,1,nrc),ureac(nrc)(1:j2),
     $        nrk(1,nrc)
              write (nttyo,1060) i,rkb(i,1,nrc),ureac(nrc)(1:j2),
     $        nrk(1,nrc)
 1060         format (/' * Error - (EQ6/chkinz) Forward rate constant',
     $        /7x,i2,' has a value of ',1pe11.4,' for reactant'
     $        /7x,a,'. A negative rate constant is not valid for a',
     $        /7x,'forward rate code of ',i2,'.')
              ier = ier + 1
            endif
          enddo
c
          if (jcode(nrc) .gt. 1) then
            write (noutpt,1070) nrk(1,nrc),ureac(nrc)(1:j2)
            write (nttyo,1070) nrk(1,nrc),ureac(nrc)(1:j2)
 1070       format(/' * Error - (EQ6/chkinz) Have specified a',
     $      /7x,'forward rate law code of ',i2,' for reactant',
     $      /7x,a,'. This rate law incorporates an affinity',
     $      /7x,'dependence and can only be used for a reactant which',
     $      /7x,'is a mineral or solid solution.')
            ier = ier + 1
          endif
        elseif (nrk(1,nrc) .eq. 3) then
c
c         Linear rate law
c
          if (rkb(1,1,nrc) .lt. 0.) then
            write (noutpt,1060) i,rkb(1,1,nrc),ureac(nrc)(1:j2),
     $      nrk(1,nrc)
            write (nttyo,1060) i,rkb(1,1,nrc),ureac(nrc)(1:j2),
     $      nrk(1,nrc)
            ier = ier + 1
          endif
        endif
c
c       Check the parameters for the backward rate laws.
c
        if (nrk(2,nrc) .eq. -1) then
c
c         Use the forward rate law for the backward case.
c
          if (nrk(1,nrc) .ne. 2) then
            write (noutpt,1080) ureac(nrc)(1:j2),nrk(1,nrc)
            write (nttyo,1080) ureac(nrc)(1:j2),nrk(1,nrc)
 1080       format(' * Error - (EQ6/chkinz) The backward rate law',
     $      /7x,'code for ',a,' is  -1, which is invalid because',
     $      /7x,'the forward rate law code is ',i2,'.')
            ier = ier + 1
          endif
        elseif (nrk(2,nrc) .eq. 2) then
c
c         Transition state rate law.
c
          do i = 1,imech(2,nrc)
            if (rkb(i,2,nrc) .lt. 0.) then
              write (noutpt,1090) i,rkb(i,2,nrc),ureac(nrc)(1:j2),
     $        nrk(2,nrc)
              write (nttyo,1090) i,rkb(i,2,nrc),ureac(nrc)(1:j2),
     $        nrk(2,nrc)
 1090         format (/' * Error - (EQ6/chkinz) Backward rate constant',
     $        /7x,i2,' has a value of ',1pe11.4,' for reactant'
     $        /7x,a,'. A negative rate constant is not valid for a',
     $        /7x,'forward rate code of ',i2,'.')
              write (noutpt,1060) ureac(nrc)(1:j2)
              write (nttyo,1060) ureac(nrc)(1:j2)
              ier = ier + 1
            endif
          enddo
c
          if (jcode(nrc) .gt. 1) then
            write (noutpt,1100) nrk(2,nrc),ureac(nrc)(1:j2)
            write (nttyo,1100) nrk(2,nrc),ureac(nrc)(1:j2)
 1100       format(/' * Error - (EQ6/chkinz) Have specified a',
     $      /7x,'backward rate law code of ',i2,' for reactant',
     $      /7x,a,'. This rate law incorporates an affinity',
     $      /7x,'dependence and can only be used for a reactant which',
     $      /7x,'is a mineral or solid solution.')
            write (noutpt,1070)
            write (nttyo,1070)
            ier = ier + 1
          endif
c
        elseif (nrk(2,nrc) .eq. 3) then
c
c         Linear  rate law.
c
          if (rkb(1,2,nrc) .lt. 0.) then
            write (noutpt,1090) i,rkb(1,2,nrc),ureac(nrc)(1:j2),
     $      nrk(2,nrc)
            write (nttyo,1090) i,rkb(1,2,nrc),ureac(nrc)(1:j2),
     $      nrk(2,nrc)
            ier = ier + 1
          endif
c
          if (jcode(nrc) .gt. 1) then
            if (nrk(1,nrc).eq.3 .or. nrk(1,nrc).eq.1) then
              write (noutpt,1110) ureac(nrc)(1:j2),nrk(1,nrc)
              write (nttyo,1110) ureac(nrc)(1:j2),nrk(1,nrc)
 1110         format(' * Error - (EQ6/chkinz) The forward rate law',
     $        /7x,'code for reactant ',a,' may not be 3 if the',
     $        /7x,'backward rate law code is',i2,'. This is permitted',
     $        /7x,'only for a reactant which is a pure mineral or',
     $        /7x,'a solid solution..')
              ier = ier + 1
            endif
          endif
c
        endif
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Matrix variable values.
c
c        zvclgi = logarithmic basis variable, starting estimate
c           or pick-up value
c
      do kcol = 1,kxt
        if (zvclgi(kcol) .le. -99999.) then
          if (uzveci(kcol)(1:8).ne.uspeca(no2gaq)(1:8) .and.
     $      uzveci(kcol)(1:8).ne.uspeca(nelect)(1:8)) then
            call fmspnx(jlen,uzveci(kcol),uspn56)
            write (noutpt,1200) uspn56(1:jlen)
            write (nttyo,1200) uspn56(1:jlen)
 1200       format(/' * Error - (EQ6/chkinz) The basis species ',a,
     $      /7x,'has a log number of moles (zvclgi) value less than',
     $      ' or equal to -99999.',/7x,'on the input file.')
            ier = ier + 1
          endif
        endif
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      end


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 performing stochastic sampling.
!>
!> This package contains a Linear Congruential Random Number Generator (LCRNG)
!> that is contained in MCNP. The LCRNG format is common and could easily be
!> extended to have other LCRNGs. The LCRNG is initialized by providing the
!> routine with an initial seed. In addition to the LCRNG, several statistical
!> distribution functions have been created to allow for the sampling of
!> non-uniform distrubitions.
!>
!> @par EXAMPLES
!> @code
!> PROGRAM ExampleRNG
!>   TYPE(StochasticSamplingType) :: myRNG
!>   INTEGER(SIK) :: randint
!>   REAL(SDK) :: randvar
!>
!>   CALL myRNG%init(RNG_MCNP_STD,SEED0=5_SLK**19)
!>   randvar=myRNG%rng()
!>   randvar=myRNG%uniform(-1.0,1.0)
!>   randvar=myRNG%exponential(0.75)
!>   randvar=myRNG%normal(0.0,1.0)
!>   randint=myRNG%histogram( (/ 0.2, 0.5, 0.1, 0.3 /) )
!> ENDPROGRAM ExampleRNG
!> @endcode
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
MODULE StochasticSampling

USE IntrType
USE ExceptionHandler
USE ParallelEnv
IMPLICIT NONE
PRIVATE
!
! List of public members
PUBLIC :: eStochasticSampler
PUBLIC :: StochasticSamplingType
PUBLIC :: RNG_MCNP_STD
PUBLIC :: RNG_LEcuyer1
PUBLIC :: RNG_LEcuyer2
PUBLIC :: RNG_LEcuyer3

! Ennumeration list for RNG types
INTEGER(SIK),PARAMETER :: RNG_MCNP_STD=1
INTEGER(SIK),PARAMETER :: RNG_LEcuyer1=2
INTEGER(SIK),PARAMETER :: RNG_LEcuyer2=3
INTEGER(SIK),PARAMETER :: RNG_LEcuyer3=4

!> Maximum length for the name of a random number generator
INTEGER(SIK),PARAMETER :: MAX_RNG_NAME_LEN=8

!> Pi
REAL(SRK),PARAMETER,PRIVATE :: PI=3.141592653589793_SRK

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

!> Add description
TYPE :: RNGdataType
  !> Random Number Data
  INTEGER(SLK) :: RNmult=-1
  !> Random Number Additive Term
  INTEGER(SLK) :: RNadd=-1
  !> log2 of Modulus, must be < 64
  INTEGER(SIK) :: RNlog2mod=-1
  !> Random Number Stride
  INTEGER(SLK) :: RNstride=-1
  !> Random Number Initial Seed
  INTEGER(SLK) :: RNseed0=-1
  !> Random Number Name
  CHARACTER(LEN=MAX_RNG_NAME_LEN) :: name=''
ENDTYPE

!> Add description
INTEGER(SIK),PRIVATE,PARAMETER :: nRN=4

!> Add description
TYPE(RNGdataType),PRIVATE,PARAMETER :: generators(nRN)=(/ &
    RNGdataType(              5_SLK**19, 0_SLK, 48, 152917_SLK, 5_SLK**19, 'mcnp std'),  &
    RNGdataType(9219741426499971445_SLK, 1_SLK, 63, 152917_SLK, 1_SLK,     'LEcuyer1'),  &
    RNGdataType(2806196910506780709_SLK, 1_SLK, 63, 152917_SLK, 1_SLK,     'LEcuyer2'),  &
    RNGdataType(3249286849523012805_SLK, 1_SLK, 63, 152917_SLK, 1_SLK,     'LEcuyer3') /)
!                           mult              add  log2mod  stride   seed0        name

!> Add description
TYPE :: StochasticSamplingType
  !> Initialization status
  LOGICAL(SBK) :: isInit=.FALSE.
  !> Random Number Seed
  INTEGER(SLK) :: RNseed=-1
  !> Random Number Data
  INTEGER(SLK) :: RNmult=-1
  !> Random Number Additive Term
  INTEGER(SLK) :: RNadd=-1
  !> Random Number Mask
  INTEGER(SLK) :: RNmask=-1
  !> Random Number Mod
  INTEGER(SLK) :: RNmod=-1
  !> Random Number Normilization
  REAL(SDK) :: RNnorm=-1
  !> Random Number Count
  INTEGER(SLK) :: counter=0
!
!List of type bound procedures
  CONTAINS
    !> @copybrief StochasticSampling::init_Sampler
    !> @copydetails StochasticSampling::init_Sampler
    PROCEDURE,PASS :: init => init_Sampler
    !> @copybrief StochasticSampling::clear_Sampler
    !> @copydetails StochasticSampling::clear_Sampler
    PROCEDURE,PASS :: clear => clear_Sampler
    !> @copybrief StochasticSampling::rng_Sampler
    !> @copydetails StochasticSampling::rng_Sampler
    PROCEDURE,PASS :: rng => rng_Sampler
    !> @copybrief StochasticSampling::unif_Sampler
    !> @copydetails StochasticSampling::unif_Sampler
    PROCEDURE,PASS :: uniform => unif_Sampler
    !> @copybrief StochasticSampling::exp_Sampler
    !> @copydetails StochasticSampling::exp_Sampler
    PROCEDURE,PASS :: exponential => exp_Sampler
    !> @copybrief StochasticSampling::norm_Sampler
    !> @copydetails StochasticSampling::norm_Sampler
    PROCEDURE,PASS :: normal => norm_Sampler
    !> @copybrief StochasticSampling::lognorm_Sampler
    !> @copydetails StochasticSampling::lognorm_Sampler
    PROCEDURE,PASS :: lognormal => lognorm_Sampler
    !> @copybrief StochasticSampling::maxw_Sampler
    !> @copydetails StochasticSampling::maxw_Sampler
    PROCEDURE,PASS :: maxwellian => maxw_Sampler
    !> @copybrief StochasticSampling::watt_Sampler
    !> @copydetails StochasticSampling::watt_Sampler
    PROCEDURE,PASS :: watt => watt_Sampler
    !> @copybrief StochasticSampling::evap_Sampler
    !> @copydetails StochasticSampling::evap_Sampler
    PROCEDURE,PASS :: evaporation => evap_Sampler
!  Kalbach-mann not yet implemented
!      !> @copybrief StochasticSampling::kalb_Sampler
!      !> @copydetails StochasticSampling::kalb_Sampler
!      PROCEDURE,PASS :: kalbachmann => kalb_Sampler
!
    !> @copybrief StochasticSampling::nhist_Sampler
    !> @copydetails StochasticSampling::nhist_Sampler
    PROCEDURE,PASS :: histogram => nhist_Sampler
    !> @copybrief StochasticSampling::uhist_Sampler
    !> @copydetails StochasticSampling::uhist_Sampler
    PROCEDURE,PASS :: unnormhistogram => uhist_Sampler
    !> @copybrief StochasticSampling::nchist_Sampler
    !> @copydetails StochasticSampling::nchist_Sampler
    PROCEDURE,PASS :: conthistogram => nchist_Sampler
    !> @copybrief StochasticSampling::uchist_Sampler
    !> @copydetails StochasticSampling::uchist_Sampler
    PROCEDURE,PASS :: unnormconthistogram => uchist_Sampler
    !> @copybrief StochasticSampling::npwl_Sampler
    !> @copydetails StochasticSampling::npwl_Sampler
    PROCEDURE,PASS :: pwlinear => npwl_Sampler
    !> @copybrief StochasticSampling::upwl_Sampler
    !> @copydetails StochasticSampling::upwl_Sampler
    PROCEDURE,PASS :: unnormpwlinear => upwl_Sampler
    !> @copybrief StochasticSampling::reject_Sampler
    !> @copydetails StochasticSampling::reject_Sampler
    PROCEDURE,PASS :: rejection => reject_Sampler
    !> @copybrief StochasticSampling::rejectarg_Sampler
    !> @copydetails StochasticSampling::rejectarg_Sampler
    PROCEDURE,PASS :: rejectionarg => rejectarg_Sampler
    !> @copybrief StochasticSampling::pwlreject_Sampler
    !> @copydetails StochasticSampling::pwlreject_Sampler
!       Implemented but buggy
!      PROCEDURE,PASS :: pwlrejection => pwlreject_Sampler
ENDTYPE StochasticSamplingType

!> Exception Handler for use in MatrixTypes
TYPE(ExceptionHandlerType),SAVE :: eStochasticSampler
!
!===============================================================================
CONTAINS
!
!-------------------------------------------------------------------------------
!> @brief Constructor for a stochastic sampler
!> @param sampler the variable to initialize
!> @param seed0 is the initial seed
!>
!> There are essentially only a few ways to use this routine. If it is not
!> called in this way then the routine will return the sampler uninitialized.
!>
!> @code
!> TYPE(StochasticSampler) :: sampler
!> CALL sampler%initialize(19073486328125_SLK)
!> @endcode
!>
SUBROUTINE init_Sampler(sampler,RNGid,seed0,skip,MPIparallelEnv,OMPparallelEnv)
  CHARACTER(LEN=*),PARAMETER :: myName='init_Sampler'
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  INTEGER(SIK),INTENT(IN) :: RNGid
  INTEGER(SLK),INTENT(IN),OPTIONAL :: seed0
  INTEGER(SLK),INTENT(IN),OPTIONAL :: skip
  TYPE(MPI_EnvType),INTENT(IN),OPTIONAL :: MPIparallelEnv
  TYPE(OMP_EnvType),INTENT(IN),OPTIONAL :: OMPparallelEnv

  INTEGER(SIK) :: mpirank,omprank,nproc,nthread
  INTEGER(SLK) :: myskip,period
  TYPE(RNGdataType) :: RNGdata

  mpirank=0
  omprank=0
  nproc=1
  nthread=1
  myskip=0_SLK

  RNGdata=generators(RNGid)

  IF(PRESENT(MPIparallelEnv)) THEN
    IF(MPIparallelEnv%isInit()) THEN
      mpirank=MPIparallelEnv%rank
      nproc=MPIparallelEnv%nproc
    ELSE
      CALL eStochasticSampler%raiseDebug(modName//'::'//myName// &
          ' - MPI Env is not initialized, and will not be used.')
    ENDIF
  ENDIF
  IF(PRESENT(OMPparallelEnv)) THEN
    IF(OMPparallelEnv%isInit()) THEN
      omprank=OMPparallelEnv%rank
      nthread=OMPparallelEnv%nproc
    ELSE
      CALL eStochasticSampler%raiseDebug(modName//'::'//myName// &
          ' - OMP Env is not initialized, and will not be used.')
    ENDIF
  ENDIF

  IF(PRESENT(skip)) myskip=skip

  ! Reduced the number of bits for period where add is non-zero to prevent
  ! overflow of the period
  IF(RNGdata%RNadd == 0) THEN
    period=ISHFT(1_SLK,RNGdata%RNlog2mod-2)
  ELSE
    period=ISHFT(1_SLK,RNGdata%RNlog2mod-1)
  ENDIF

  myskip=myskip+INT(mpirank,SLK)*INT(period/INT(nproc,SLK),SLK)+ &
      INT(omprank,SLK)*INT(period/INT(nproc*nthread,SLK),SLK)

  sampler%RNseed=RNGdata%RNseed0
  ! Add checks for constraints on seed0
  IF(PRESENT(seed0)) sampler%RNseed=seed0

  IF(myskip /= 0_SLK) sampler%RNseed=RNskip(RNGdata,sampler%RNseed,myskip)

  sampler%RNmult=RNGdata%RNmult
  sampler%RNadd=RNGdata%RNadd
  sampler%RNmask=2_SLK**RNGdata%RNlog2mod-1_SLK
  sampler%RNmod=2_SLK**RNGdata%RNlog2mod
  sampler%RNnorm=1.0_SDK/2.0_SDK**RNGdata%RNlog2mod

  sampler%isInit=.TRUE.
  sampler%counter=0
ENDSUBROUTINE init_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine clears the data in a stochastic sampler type variable
!> @param sampler the type variable to clear
!>
ELEMENTAL SUBROUTINE clear_Sampler(sampler)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  sampler%RNseed=-1
  sampler%RNseed=-1
  sampler%RNmult=-1
  sampler%RNadd=-1
  sampler%RNmask=-1
  sampler%RNmod=-1
  sampler%RNnorm=-1
  sampler%counter=0
  sampler%isInit=.FALSE.
ENDSUBROUTINE clear_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns the next random number for the stochastic sampler
!> @param sampler the type variable to sample from
!>
FUNCTION rng_Sampler(sampler) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK) :: rang
  sampler%RNseed=IAND(IAND(sampler%RNmult*sampler%RNseed,sampler%RNmask)+ &
      sampler%RNadd,sampler%RNmask)
  rang=sampler%RNseed*sampler%RNnorm
  sampler%counter=sampler%counter+1
ENDFUNCTION rng_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from a uniform distribution between a
!> and b
!> @param sampler the type variable to sample from
!> @param xmin the minimum value in the uniform distribution
!> @param xmax the maximum value in the uniform distribution
!>
FUNCTION unif_Sampler(sampler,xmin,xmax) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: xmin
  REAL(SDK),INTENT(IN) :: xmax
  REAL(SDK) :: rang
  rang=0._SDK
  rang=xmin+(xmax-xmin)*sampler%rng()
ENDFUNCTION unif_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from an exponential distribution with
!> the coefficient a
!> @param sampler the type variable to sample from
!> @param a is the coefficient of the exponential f(x)=a*EXP(-a*x)
!>
FUNCTION exp_Sampler(sampler,a) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: a
  REAL(SDK) :: rang
  rang=0._SDK
  rang=-LOG(sampler%rng())/a
ENDFUNCTION exp_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from a normal distribution with
!> the mean value xbar and standard deviation sigma
!> @param sampler the type variable to sample from
!> @param mu is the mean value
!> @param sigma is the standard deviation
!>
FUNCTION norm_Sampler(sampler,mu,sigma) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: mu
  REAL(SDK),INTENT(IN) :: sigma
  REAL(SDK) :: rn1,rn2
  REAL(SDK) :: rang
  rang=0._SDK
  rn1=-LOG(sampler%rng())
  rn2=COS(2.0_SDK*PI*sampler%rng())
  rang=mu+sigma*SQRT(2.0_SDK*rn1)*rn2
ENDFUNCTION norm_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from a log-normal distribution with
!> the mean value xbar and standard deviation sigma
!> @param sampler the type variable to sample from
!> @param mu is the mean value
!> @param sigma is the standard deviation
!>
FUNCTION lognorm_Sampler(sampler,mu,sigma) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: mu
  REAL(SDK),INTENT(IN) :: sigma
  REAL(SDK) :: rang
  rang=0._SDK
  rang=EXP(sampler%normal(mu,sigma))
ENDFUNCTION lognorm_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from a Maxwellian distribution with
!> the temperature T
!> @param sampler the type variable to sample from
!> @param T is the temperature
!>
FUNCTION maxw_Sampler(sampler,T) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: T
  REAL(SDK) :: rn1,rn2,rn3
  REAL(SDK) :: rang
  rang=0._SDK
  rn1=-LOG(sampler%rng())
  rn2=-LOG(sampler%rng())
  rn3=COS(0.5_SDK*PI*sampler%rng())
  rn3=rn3*rn3
  rang=T*(rn1+rn2*rn3)
ENDFUNCTION maxw_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from a Watt fission distribution with
!> the coefficients a and b
!> @param sampler the type variable to sample from
!> @param a is the coefficient of the Watt fission spectrum
!> @param b is the coefficient of the Watt fission spectrum
!>
FUNCTION watt_Sampler(sampler,a,b) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: a
  REAL(SDK),INTENT(IN) :: b
  REAL(SDK) :: w,a2b
  REAL(SDK) :: rang
  a2b=a*a*b
  w=sampler%maxwellian(a)
  rang=w+0.25_SDK*a2b+sampler%uniform(-1.0_SDK,1.0_SDK)*SQRT(a2b*w)
ENDFUNCTION watt_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from a Evaporation distribution with
!> the coefficient theta
!> @param sampler the type variable to sample from
!> @param thata is the coefficient of the Evaporation spectrum
!>
FUNCTION evap_Sampler(sampler,theta) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: theta
  REAL(SDK) :: rn1
  REAL(SDK) :: rang
  rang=0._SDK
  rn1=sampler%rng()
  rang=-theta*LOG(rn1*sampler%rng())
ENDFUNCTION evap_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from a Discrete Histogram with
!> components y
!> @param sampler the type variable to sample from
!> @param y is a vector of the values of the histogram
!>
FUNCTION nhist_Sampler(sampler,y) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: y(:)
  REAL(SDK) :: rn1, sum
  INTEGER(SIK) :: rang
  INTEGER(SIK) :: n

  n=SIZE(y,DIM=1)
  rn1=sampler%rng()
  sum=0.0_SDK
  DO rang=1,n
    sum=sum+y(rang)
    IF(rn1 < sum) EXIT
  ENDDO
ENDFUNCTION nhist_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from an unnormalized discrete
!> histogram with components y
!> @param sampler the type variable to sample from
!> @param y is a vector of the values of the histogram
!>
FUNCTION uhist_Sampler(sampler,y) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: y(:)
  INTEGER(SIK) :: rang
  rang=sampler%histogram(y/SUM(y))
ENDFUNCTION uhist_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from a normalized continuous
!> histogram with components y and bounds x
!> @param sampler the type variable to sample from
!> @param y is a vector of the values of the histogram
!> @param x is a vector of the bounds of y
!>
FUNCTION nchist_Sampler(sampler,y,x) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: y(:)
  REAL(SDK),INTENT(IN) :: x(:)
  REAL(SDK) :: rang
  INTEGER(SIK) :: interval
  interval=0
  interval=sampler%histogram(y)
  rang=sampler%uniform(x(interval),x(interval+1))
ENDFUNCTION nchist_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from an unnormalized continuous
!> histogram with components y and bounds x
!> @param sampler the type variable to sample from
!> @param y is a vector of the values of the histogram
!> @param x is a vector of the bounds of y
!>
FUNCTION uchist_Sampler(sampler,y,x) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: y(:)
  REAL(SDK),INTENT(IN) :: x(:)
  REAL(SDK) :: rang
  rang=0._SDK
  rang=sampler%conthistogram(y/SUM(y),x)
ENDFUNCTION uchist_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from a normalized piece-wise linear
!> distribution with points x and y
!> @param sampler the type variable to sample from
!> @param y is a vector of the y values of the piece-wise linear function
!> @param x is a vector of the x values of the piece-wise linear function
!>
FUNCTION npwl_Sampler(sampler,y,x) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: y(:)
  REAL(SDK),INTENT(IN) :: x(:)
  REAL(SDK) :: rang
  REAL(SDK) :: rn1,rn2,sum
  INTEGER(SIK) :: i,n

  n=SIZE(x,DIM=1)-1

  rn1=sampler%rng()
  rn2=sampler%rng()

  sum=0.0_SDK
  DO i=1,n
    sum=sum+(y(i)+y(i+1))/2*(x(i+1)-x(i))
    IF(rn1 <sum) EXIT
  ENDDO

  IF(rn1 < y(i)/(y(i+1)+y(i))) THEN
    rang=x(i+1)-(x(i+1)-x(i))*SQRT(rn2)
  ELSE
    rang=x(i)+(x(i+1)-x(i))*SQRT(rn2)
  ENDIF
ENDFUNCTION npwl_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from a unnormalized piece-wise linear
!> distribution with points x and y
!> @param sampler the type variable to sample from
!> @param y is a vector of the y values of the piece-wise linear function
!> @param x is a vector of the x values of the piece-wise linear function
!>
FUNCTION upwl_Sampler(sampler,y,x) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: y(:)
  REAL(SDK),INTENT(IN) :: x(:)
  REAL(SDK) :: rang
  REAL(SDK) :: sum
  INTEGER(SIK) :: i,n

  n=SIZE(y,DIM=1)-1

  sum=0.0_SDK
  DO i=1,n
    sum=sum+(y(i)+y(i+1))/2*(x(i+1)-x(i))
  ENDDO
  rang=sampler%pwlinear(y/sum,x)
ENDFUNCTION upwl_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from an arbitrary function func using
!> rejection sampling
!> @param sampler the type variable to sample from
!> @param func is the function which is sampled
!> @param xmin is the minimum value of x
!> @param xmax is the maximum value of x
!> @param ymax is the maximum value of func in the range xmin to xmax
!>
FUNCTION reject_Sampler(sampler,func,xmin,xmax,ymax) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: xmin
  REAL(SDK),INTENT(IN) :: xmax
  REAL(SDK),INTENT(IN) :: ymax
  REAL(SDK) :: rang

  INTERFACE
    FUNCTION func(x)
      IMPORT :: SDK
      REAL(SDK),INTENT(IN) :: x
      REAL(SDK) :: func
    ENDFUNCTION
  ENDINTERFACE

  DO
    rang=sampler%uniform(xmin,xmax)
    IF(sampler%uniform(0.0_SDK,ymax) <= func(rang)) RETURN
  ENDDO
ENDFUNCTION reject_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from an arbitrary function func with
!> extra argument vector arg using rejection sampling
!> @param sampler the type variable to sample from
!> @param func is the function which is sampled
!> @param xmin is the minimum value of x
!> @param xmax is the maximum value of x
!> @param ymax is the maximum value of func in the range xmin to xmax
!> @param arg is a vector of arguments that gets passed to func
!>
FUNCTION rejectarg_Sampler(sampler,func,xmin,xmax,ymax,arg) RESULT(rang)
  CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
  REAL(SDK),INTENT(IN) :: xmin
  REAL(SDK),INTENT(IN) :: xmax
  REAL(SDK),INTENT(IN) :: ymax
  REAL(SDK),INTENT(IN) :: arg(:)
  REAL(SDK) :: rang

  INTERFACE
    FUNCTION func(x,arg)
      IMPORT :: SDK
      REAL(SDK),INTENT(IN) :: x
      REAL(SDK),INTENT(IN) :: arg(:)
      REAL(SDK) :: func
    ENDFUNCTION
  ENDINTERFACE

  DO
    rang=sampler%uniform(xmin,xmax)
    IF(sampler%uniform(0.0_SDK,ymax) <= func(rang,arg)) RETURN
  ENDDO
ENDFUNCTION rejectarg_Sampler
!!
!!-------------------------------------------------------------------------------
!!> @brief Routine returns a random number from an arbitrary function func using
!!> rejection sampling which is bound by a piece-wise linear function
!!> @param sampler the type variable to sample from
!!> @param func is the function which is sampled
!!> @param yval is the y components of a piece-wise linear function bounding func
!!> @param xval is the x components of a piece-wise linear function bounding func
!!> @param c is an optional scalar which scales the piece-wise linear fucntion
!!>        bounding func.  It is important to note that if c is present the pwl
!!>        function is assumed to be normalized, if it is not present the pwl
!!>        function is not scaled and used as is
!!>
!*** this function is commented out because it is buggy and isn't a needed feature right now
!    FUNCTION pwlreject_Sampler(sampler,func,yval,xval,c) RESULT(rang)
!      CLASS(StochasticSamplingType),INTENT(INOUT) :: sampler
!      REAL(SDK),INTENT(IN) :: yval(:)
!      REAL(SDK),INTENT(IN) :: xval(:)
!      REAL(SDK),INTENT(IN),OPTIONAL :: c
!      REAL(SDK) :: rang
!      REAL(SDK) :: g, mult
!      REAL(SDK),ALLOCATABLE :: yscaled(:)
!      INTEGER(SIK) :: i,n
!
!      INTERFACE
!        FUNCTION func(x)
!          IMPORT :: SDK
!          REAL(SDK),INTENT(IN) :: x
!          REAL(SDK) :: func
!        ENDFUNCTION
!      ENDINTERFACE
!
!      n=SIZE(xval,DIM=1)-1
!      ALLOCATE(yscaled(n+1))
!
!      IF(PRESENT(c)) THEN
!        mult=c
!        yscaled=yval
!      ELSE
!        mult=0.0_SDK
!        DO i=1,n
!          mult=mult+(yval(i)+yval(i+1))/2*(xval(i+1)-xval(i))
!        ENDDO
!        yscaled=yval/mult
!      ENDIF
!
!      DO
!        rang=sampler%pwlinear(yscaled,xval)
!        DO i=1,n
!          IF(xval(i) >= rang) then
!            g=(yscaled(i)-yscaled(i-1))/(xval(i)-xval(i-1))*(rang-xval(i-1))+yscaled(i-1)
!            EXIT
!          ENDIF
!        ENDDO
!        write(*,*) sampler%uniform(0.0_SDK,mult*g)
!        IF (sampler%uniform(0.0_SDK,mult*g)<=func(rang)) RETURN
!      ENDDO
!    ENDFUNCTION pwlreject_Sampler
!
!-------------------------------------------------------------------------------
!> @brief Routine returns a random number from an arbitrary function func using
!> rejection sampling which is bound by a piece-wise linear function
!> @param sampler the type variable to sample from
!> @param func is the function which is sampled
!> @param yval is the y components of a piece-wise linear function bounding func
!> @param xval is the x components of a piece-wise linear function bounding func
!> @param c is an optional scalar which scales the piece-wise linear fucntion
!>        bounding func.  It is important to note that if c is present the pwl
!>        function is assumed to be normalized, if it is not present the pwl
!>        function is not scaled and used as is
!>
PURE FUNCTION RNskip(RNGdata,seed0,skip) RESULT(newseed)
  TYPE(RNGdataType),INTENT(IN) :: RNGdata
  INTEGER(SLK),INTENT(IN) :: seed0
  INTEGER(SLK),INTENT(IN) :: skip
  ! Local Variables
  INTEGER(SLK) :: newseed
  INTEGER(SLK) :: nskip,gen,g,inc,c,gp,rn,period,mask

  mask=ISHFT(NOT(0_SLK),RNGdata%RNlog2mod-64)
  IF(RNGdata%RNadd == 0) THEN
    period=ISHFT(1_SLK,RNGdata%RNlog2mod-2)
  ELSE
    period=ISHFT(1_SLK,RNGdata%RNlog2mod)
  ENDIF

  nskip=skip
  DO WHILE(nskip < 0_SLK)
    nskip=nskip+period
  ENDDO

  nskip=IAND(nskip,mask)
  gen=1
  g=RNGdata%RNmult
  inc=0
  c=RNGdata%RNadd
  DO WHILE(nskip > 0_SLK)
    IF(BTEST(nskip,0))  THEN
      gen=IAND(gen*g,mask)
      inc=IAND(inc*g,mask)
      inc=IAND(inc+c,mask)
    ENDIF
    gp=IAND(g+1,mask)
    g=IAND(g*g,mask)
    c=IAND(gp*c,mask)
    nskip=ISHFT(nskip,-1)
  ENDDO
  rn=IAND(gen*seed0,mask)
  rn=IAND(rn+inc,mask)
  newseed=rn
ENDFUNCTION RNskip
!
ENDMODULE StochasticSampling

max_stars_repo_name: picmc/Futility

max_issues_repo_name: picmc/Futility

max_forks_repo_name: wgurecky/Futility


text:
      INTEGER FUNCTION ignpoi(mu)
C**********************************************************************
C
C     INTEGER FUNCTION IGNPOI( MU )
C
C                    GENerate POIsson random deviate
C
C
C                              Function
C
C
C     Generates a single random deviate from a Poisson
C     distribution with mean MU.
C
C
C                              Arguments
C
C
C     MU --> The mean of the Poisson distribution from which
C            a random deviate is to be generated.
C                              DOUBLE PRECISION MU
C     JJV                    (MU >= 0.0)
C
C     IGNPOI <-- The random deviate.
C                              INTEGER IGNPOI (non-negative)
C
C
C                              Method
C
C
C     Renames KPOIS from TOMS as slightly modified by BWB to use RANF
C     instead of SUNIF.
C
C     For details see:
C
C               Ahrens, J.H. and Dieter, U.
C               Computer Generation of Poisson Deviates
C               From Modified Normal Distributions.
C               ACM Trans. Math. Software, 8, 2
C               (June 1982),163-179
C
C**********************************************************************
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C                                                                      C
C     P O I S S O N  DISTRIBUTION                                      C
C                                                                      C
C                                                                      C
C**********************************************************************C
C**********************************************************************C
C                                                                      C
C     FOR DETAILS SEE:                                                 C
C                                                                      C
C               AHRENS, J.H. AND DIETER, U.                            C
C               COMPUTER GENERATION OF POISSON DEVIATES                C
C               FROM MODIFIED NORMAL DISTRIBUTIONS.                    C
C               ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. C
C                                                                      C
C     (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)  C
C                                                                      C
C**********************************************************************C
C
C      INTEGER FUNCTION IGNPOI(IR,MU)
C
C     INPUT:  IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
C             MU=MEAN MU OF THE POISSON DISTRIBUTION
C     OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
C
C
C
C     MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR CASE B
C     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
C     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
C
C
C
C     SEPARATION OF CASES A AND B
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION mu
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION a0,a1,a2,a3,a4,a5,a6,a7,b1,b2,c,c0,c1,c2,c3,
     $     d,del,difmuk,e,
     $     fk,fx,fy,g,muold,muprev,omega,p,p0,px,py,q,s,t,u,v,x,xx
C     JJV I added a variable 'll' here - it is the 'l' for CASE A
      INTEGER j,k,kflag,l,ll,m
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION fact(10),pp(35)
C     ..
C     .. External Functions ..
      DOUBLE PRECISION ranf,sexpo,snorm
      EXTERNAL ranf,sexpo,snorm
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,log,exp,float,idint,max0,min0,dsign,sqrt
C     ..
C     JJV added this for case: mu unchanged
C     .. Save statement ..
      SAVE s, d, l, ll, omega, c3, c2, c1, c0, c, m, p, q, p0,
     +     a0, a1, a2, a3, a4, a5, a6, a7, fact, muprev, muold
C     ..
C     JJV end addition - I am including vars in Data statements
C     .. Data statements ..
C     JJV changed initial values of MUPREV and MUOLD to -1.0E37
C     JJV if no one calls IGNPOI with MU = -1.0E37 the first time,
C     JJV the code shouldn't break
      DATA muprev,muold/-1.0E37,-1.0E37/
      DATA a0,a1,a2,a3,a4,a5,a6,a7/-.5,.3333333,-.2500068,.2000118,
     +     -.1661269,.1421878,-.1384794,.1250060/
      DATA fact/1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880./
C     ..
C     .. Executable Statements ..

      IF (mu.EQ.muprev) GO TO 10
      IF (mu.LT.10.0) GO TO 120
C
C     C A S E  A. (RECALCULATION OF S,D,LL IF MU HAS CHANGED)
C
C     JJV This is the case where I changed 'l' to 'll'
C     JJV Here 'll' is set once and used in a comparison once

      muprev = mu
      s = sqrt(mu)
      d = 6.0*mu*mu
C
C             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
C             PROBABILITIES FK WHENEVER K >= M(MU). LL=IDINT(MU-1.1484)
C             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
C
      ll = idint(mu-1.1484)
C
C     STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
C
   10 g = mu + s*snorm()
      IF (g.LT.0.0) GO TO 20
      ignpoi = idint(g)
C
C     STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
C
      IF (ignpoi.GE.ll) RETURN
C
C     STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
C
      fk = dble(ignpoi)
      difmuk = mu - fk
      u = ranf()
      IF (d*u.GE.difmuk*difmuk*difmuk) RETURN
C
C     STEP P. PREPARATIONS FOR STEPS Q AND H.
C             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
C             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
C             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
C             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
C             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
C
   20 IF (mu.EQ.muold) GO TO 30
      muold = mu
      omega = .3989423/s
      b1 = .4166667E-1/mu
      b2 = .3*b1*b1
      c3 = .1428571*b1*b2
      c2 = b2 - 15.*c3
      c1 = b1 - 6.*b2 + 45.*c3
      c0 = 1. - b1 + 3.*b2 - 15.*c3
      c = .1069/mu
   30 IF (g.LT.0.0) GO TO 50
C
C             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
C
      kflag = 0
      GO TO 70
C
C     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
C
   40 IF (fy-u*fy.LE.py*exp(px-fx)) RETURN
C
C     STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
C             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
C             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
C
   50 e = sexpo()
      u = ranf()
      u = u + u - 1.0
      t = 1.8d0 + dsign(e,u)
      IF (t.LE. (-.6744)) GO TO 50
      ignpoi = idint(mu+s*t)
      fk = dble(ignpoi)
      difmuk = mu - fk
C
C             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
C
      kflag = 1
      GO TO 70
C
C     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
C
   60 IF (c*abs(u).GT.py*exp(px+e)-fy*exp(fx+e)) GO TO 50
      RETURN
C
C     STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
C             CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
C
   70 IF (ignpoi.GE.10) GO TO 80
      px = -mu
      py = mu**ignpoi/fact(ignpoi+1)
      GO TO 110
C
C             CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
C             A0-A7 FOR ACCURACY WHEN ADVISABLE
C             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)
C
   80 del = .8333333E-1/fk
      del = del - 4.8*del*del*del
      v = difmuk/fk
      IF (abs(v).LE.0.25) GO TO 90
      px = fk*log(1.0+v) - difmuk - del
      GO TO 100

   90 px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) -
     +     del
  100 py = .3989423/sqrt(fk)
  110 x = (0.5-difmuk)/s
      xx = x*x
      fx = -0.5*xx
      fy = omega* (((c3*xx+c2)*xx+c1)*xx+c0)
      IF (kflag) 40,40,60
C
C     C A S E  B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
C
C     JJV changed MUPREV assignment from 0.0 to initial value
  120 muprev = -1.0E37
C     Jpc 1999: the next lines seams to produce a bug
C     and I finaly commented them out
C     IF (mu.EQ.muold) GO TO 130
C     JJV added argument checker here
C     JJV added line label here
C 125  muold = mu
      m = max0(1,idint(mu))
      l = 0
      p = exp(-mu)
      q = p
      p0 = p
C
C     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
C
  130 u = ranf()
      ignpoi = 0
      IF (u.LE.p0) RETURN
C
C     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
C             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
C             (0.458=PP(9) FOR MU=10)
C
      IF (l.EQ.0) GO TO 150
      j = 1
      IF (u.GT.0.458) j = min0(l,m)
      DO 140 k = j,l
          IF (u.LE.pp(k)) GO TO 180
  140 CONTINUE
      IF (l.EQ.35) GO TO 130
C
C     STEP C. CREATION OF NEW POISSON PROBABILITIES P
C             AND THEIR CUMULATIVES Q=PP(K)
C
  150 l = l + 1
      DO 160 k = l,35
          p = p*mu/float(k)
          q = q + p
          pp(k) = q
          IF (u.LE.q) GO TO 170
  160 CONTINUE
      l = 35
      GO TO 130

  170 l = k
  180 ignpoi = k
      RETURN

      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx