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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!! DOUBLE PRECISION VERSION

!!!!! List of routines

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

!dcsr2csr_up,  zcsr2csr_up  :: extract upper triangular part of a full csr matrix
!dcsr2csr_low, zcsr2csr_low :: extract lower triangular part of a full csr matrix

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

!dcsr_transpose, zcsr_transpose :: transpose a csr matrix
!zcsr_htranspose                :: transpose conjugate a csr matrix

!dcsr_uplo_to_csr, zcsr_uplo_to_csr  :: convert a upper or lower csr matrix into a full csr one
!zhcsr_uplo_to_csr                   :: convert a upper or lower (Hermitian) csr matrix into a full csr one

!dcsr_convert_upper                 :: convert a upper,lower or full csr matrix into an upper csr matrix
!zcsr_convert_upper                 :: convert a upper,lower or full csr matrix into an upper csr matrix (complex symmetric)
!zhcsr_convert_full                 :: convert a upper,lower or full csr Hermitian matrix into a full csr matrix

!dcoo2csr :: transform coordinate format into csr format
!zcoo2csr :: transform coordinate format into csr format

!dsubmat                            :: extract csr submatrix from csr matrix
!zsubmat                            :: extract csr submatrix from csr matrix

subroutine daddcsr(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 double precision
  !
  !  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) DOUBLE PRECISION
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa DOUBLE PRECISION)
  !
  !  beta  (input) DOUBLE PRECISION
  !
  !  sb,isb,jsb (input) CSR format for the matrix B (sb DOUBLE PRECISION)
  !
  !  sc,isc,jsc (input/output) CSR format for the matrix C (sc DOUBLE PRECISION).
  !             see definition of 'opt' to determine input or output.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  double precision,dimension(*) :: sa,sb,sc
  integer,dimension(*) :: isa,isb,isc,jsa,jsb,jsc
  double precision :: 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 daddcsr

subroutine zaddcsr(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 double precision
  !
  !  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 DOUBLE PRECISION
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa COMPLEX DOUBLE PRECISION)
  !
  !  beta  (input) COMPLEX DOUBLE PRECISION
  !
  !  sb,isb,jsb (input) CSR format for the matrix B (sb COMPLEX DOUBLE PRECISION)
  !
  !  sc,isc,jsc (input/output) CSR format for the matrix C (sc COMPLEX DOUBLE PRECISION).
  !             see definition of 'opt' to determine input or output.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  complex(kind=kind(1.0d0)),dimension(*) :: sa,sb,sc
  integer,dimension(*) :: isa,isb,isc,jsa,jsb,jsc
  complex(kind=kind(1.0d0)) :: 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 zaddcsr

subroutine zdaddcsr(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
  !
  !  C is complex double precision, A and B are real double precision, (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 DOUBLE PRECISION
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa DOUBLE PRECISION)
  !
  !  beta  (input) COMPLEX DOUBLE PRECISION
  !
  !  sb,isb,jsb (input) CSR format for the matrix B (sb DOUBLE PRECISION)
  !
  !  sc,isc,jsc (input/output) CSR format for the matrix C (sc COMPLEX DOUBLE PRECISION).
  !             see definition of 'opt' to determine input or output.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  double precision,dimension(*) :: sa,sb
  complex(kind=kind(1.0d0)),dimension(*) :: sc
  integer,dimension(*) :: isa,isb,isc,jsa,jsb,jsc
  complex(kind=kind(1.0d0)) :: 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 zdaddcsr

subroutine zinc_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 double 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 DOUBLE PRECISION
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa COMPLEX DOUBLE PRECISION)
  !
  !
  !  sb,isb,jsb (input/output) CSR format for the matrix B (sb COMPLEX DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  complex(kind=kind(1.0d0)),dimension(*) :: sa,sb
  integer,dimension(*) :: isa,isb,jsa,jsb
  complex(kind=kind(1.0d0)) :: 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 zinc_addcsr

subroutine dcsr2csr_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 DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  double precision,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  double precision,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 dcsr2csr_up

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

subroutine zcsr2csr_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 DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex(kind=kind(1.0d0)),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 zcsr2csr_up

subroutine dcsr2csr_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 DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  double precision,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  double precision,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 dcsr2csr_low
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine zcsr2csr_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 DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex(kind=kind(1.0d0)),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 zcsr2csr_low

SUBROUTINE dcsrmm(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 DOUBLE 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) DOUBLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE PRECISION)
  !
  !  X     (input) DOUBLE PRECISION
  !        matrix of size (Mxrhs) or (Nxrhs) if TRANS='T'
  !
  !  beta  (input) DOUBLE PRECISION
  !
  !  B     (input/output) DOUBLE 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
  double precision :: alpha,beta
  double precision,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  double precision,dimension(*):: x
  double precision,dimension(*) ::b
!!!!!!!!

  integer ::i,j,k,s,n1
  double precision,parameter :: DZERO=0.0d0

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

!!!!!!!! Initialization

if (beta/=DZERO) then
call DSCAL(n1*rhs,beta,b,1)
else
b(1:n1*rhs)=DZERO
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)
              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 dcsrmm

SUBROUTINE zcsrmm(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 DOUBLE 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) DOUBLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE PRECISION)
  !
  !  X     (input) DOUBLE PRECISION
  !        matrix of size (Mxrhs) or (Nxrhs) if TRANS='T','C'.
  !
  !  beta  (input) DOUBLE PRECISION
  !
  !  B     (input/output) DOUBLE 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 (kind=kind(1.0d0)) :: alpha,beta
  complex (kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex (kind=kind(1.0d0)),dimension(*):: x
  complex (kind=kind(1.0d0)),dimension(*) ::b
!!!!!!!!

  integer ::i,j,k,s,n1
  complex (kind=kind(1.0d0)),parameter :: ZEROC=(0.0d0,0.0d0)

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

!!!!!!!! Initialization

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

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

     if (TRANS=='N') then
!call mkl_set_num_threads(1)
!        !$OMP PARALLEL DO SHARED(a,ia,ja,alpha,x,b,n,rhs) PRIVATE(i,k,j)
        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
              ! call ZAXPY(rhs,alpha*a(k),x(j,1),N,b(i,1),N) !!!<<< here
              !b(i,r)=b(i,r)+alpha*a(k)*x(j,r)
              !b(i,r)=b(i,r)+a(k)*x(j,r)
           end do
        end do
!        !$OMP END PARALLEL DO
!call mkl_set_num_threads(4)
     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)
                 enddo
           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)
              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)
              if (j/=i) b(i+(s-1)*N)=b(i+(s-1)*N)+alpha*conjg(a(k))*x(j+(s-1)*M)
              enddo
           end do
        end do

     end if

  end if

END SUBROUTINE zcsrmm

SUBROUTINE zhcsrmm(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 DOUBLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX DOUBLE PRECISION)
  !
  !  X     (input) COMPLEX DOUBLE PRECISION
  !        matrix of size (Mxrhs) or (Nxrhs) if TRANS='T','C'.
  !
  !  beta  (input) COMPLEX DOUBLE PRECISION
  !
  !  B     (input/output) COMPLEX DOUBLE 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(kind=kind(1.0d0)) :: alpha,beta
  complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex(kind=kind(1.0d0)),dimension(*) :: x
  complex(kind=kind(1.0d0)),dimension(*) :: b
!!!!!!!!

  integer ::i,j,k,s,n1
  double precision,parameter :: DZERO=0.0d0
  complex(kind=kind(1.0d0)),parameter :: ZEROC=(DZERO,DZERO)

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

!!!!!!!! Initialization

if (beta/=ZEROC) then
call ZSCAL(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)
                 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

     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)
                 enddo
           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)
        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*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)
              enddo
           end do
        end do

     end if

  end if

end SUBROUTINE zhcsrmm

SUBROUTINE dzcsrmm(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
  !
  !      DOUBLE PRECISION version
  !
  !      Note: same as zcsrmm 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 DOUBLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a DOUBLE PRECISION)
  !
  !  X     (input) COMPLEX DOUBLE PRECISION
  !        matrix of size (Mxrhs) or (Nxrhs) if TRANS='T','C'.
  !
  !  beta  (input) COMPLEX DOUBLE PRECISION
  !
  !  B     (input/output) COMPLEX DOUBLE 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 (kind=kind(1.0d0)) :: alpha,beta
  double precision,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex (kind=kind(1.0d0)),dimension(*):: x
  complex (kind=kind(1.0d0)),dimension(*) ::b
!!!!!!!!

  integer ::i,j,k,s,n1
  complex (kind=kind(1.0d0)),parameter :: ZEROC=(0.0d0,0.0d0)

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

!!!!!!!! Initialization

if (beta/=ZEROC) then
call ZSCAL(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)
                 enddo
           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)
                 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 dzcsrmm

subroutine dcsr_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 DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: N,M
  double precision,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  double precision,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 before
  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 dcsr_transpose

subroutine zcsr_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 DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: N,M
   complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
   complex(kind=kind(1.0d0)),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 zcsr_transpose

subroutine zcsr_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 DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: N,M
   complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
   complex(kind=kind(1.0d0)),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)))=conjg(a(k))
        tab(ja(k))=tab(ja(k))+1
     end do
  end do

deallocate(tab)

end subroutine zcsr_htranspose

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

subroutine dcsr_uplo_to_csr(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a double 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 DOUBLE PRECISION)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  !include 'f90_noruntime_interface.fi'
  integer :: N
  double precision,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  double precision,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,k,nnza,info,opt
  double precision,dimension(:),allocatable ::sb
  integer,dimension(:),allocatable :: sib
  integer,dimension(:),allocatable :: sjb
  double precision,parameter :: DONE=1.0d0,DZERO=0.0d0
!!!!
  nnza=ia(n+1)-1

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

!!! transpose
  call dcsr_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)=DZERO
     end do
  enddo
  !! add the two matrices to form the full matrix
  opt=0
  call daddcsr(N,N,opt,DONE,a,ia,ja,DONE,sb,sib,sjb,sa,sia,sja)

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

end subroutine dcsr_uplo_to_csr

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

subroutine zcsr_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 DOUBLE COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  !include 'f90_noruntime_interface.fi'
  integer :: N
  complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex(kind=kind(1.0d0)),dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,k,nnza,info,opt
  complex(kind=kind(1.0d0)),dimension(:),allocatable ::sb
  integer,dimension(:),allocatable :: sib
  integer,dimension(:),allocatable :: sjb
  double precision,parameter :: DONE=1.0d0,DZERO=0.0d0
  complex(kind=kind(1.0d0)),parameter :: ONEC=(DONE,DZERO),ZEROC=(DZERO,DZERO)
!!!!
  nnza=ia(n+1)-1

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

!!! transpose
  call zcsr_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 zaddcsr(N,N,opt,ONEC,a,ia,ja,ONEC,sb,sib,sjb,sa,sia,sja)

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

end subroutine zcsr_uplo_to_csr

subroutine zhcsr_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 DOUBLE COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa DOUBLE COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  !include 'f90_noruntime_interface.fi'
  integer :: N
  complex(kind=kind(1.0d0)),dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex(kind=kind(1.0d0)),dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,k,nnza,info,opt
  complex(kind=kind(1.0d0)),dimension(:),allocatable ::sb
  integer,dimension(:),allocatable :: sib
  integer,dimension(:),allocatable :: sjb
  double precision,parameter :: DONE=1.0d0,DZERO=0.0d0
  complex(kind=kind(1.0d0)),parameter :: ONEC=(DONE,DZERO),ZEROC=(DZERO,DZERO)
!!!!
  nnza=ia(n+1)-1

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

!!! transpose
  call zcsr_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 zaddcsr(N,N,opt,ONEC,a,ia,ja,ONEC,sb,sib,sjb,sa,sia,sja)

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

end subroutine zhcsr_uplo_to_csr

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

  subroutine dcsr_convert_upper(n,UPLO,sa,isa,jsa,ssa,sisa,sjsa)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !      Purpose: This subroutine convert a double real upper, loweror 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 DOUBLE PRECISION)
  !
  !  ssa,sisa,sjsa (output) CSR format for the matrix sA (sa  DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: n
  character(len=1) :: UPLO
  double precision,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 dcsr2csr_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 DCOPY(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 dcsr_transpose(N,N,sa,isa,jsa,ssa,sisa,sjsa)
    endif

end subroutine dcsr_convert_upper

  subroutine zcsr_convert_upper(n,UPLO,sa,isa,jsa,ssa,sisa,sjsa)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !      Purpose: This subroutine convert a double complex upper, loweror full CSR matrix A
  !      into a new CSR upper matrix sA (symmetric matrix)
    !
  !  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 COMPLEX DOUBLE PRECISION)
  !
  !  ssa,sisa,sjsa (output) CSR format for the matrix sA (sa  COMPLEX DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: n
  character(len=1) :: UPLO
  complex(kind=kind(1.0d0)),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 zcsr2csr_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 ZCOPY(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 zcsr_transpose(N,N,sa,isa,jsa,ssa,sisa,sjsa)
    endif

end subroutine zcsr_convert_upper

  subroutine zhcsr_convert_full(n,UPLO,sa,isa,jsa,ssa,sisa,sjsa)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a double 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 DOUBLE PRECISION)
  !
  !  ssa,sisa,sjsa (output) CSR format for the Hermitian matrix sA (sa COMPLEX DOUBLE PRECISION)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: n
  character(len=1) :: UPLO
  complex(kind=kind(1.0d0)),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 ZCOPY(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 zhcsr_uplo_to_csr(N,sa,isa,jsa,ssa,sisa,sjsa)
end if

end subroutine zhcsr_convert_full

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

subroutine dcoo2csr(n,nnz,ic,jc,c,isa,jsa,sa)
  implicit none
  integer :: n,nnz
  integer,dimension(*) :: ic,jc,isa,jsa
  double precision,dimension(*) :: c,sa
!!!
  integer :: k,k1,i,j,idum
  integer,dimension(n) :: iloc
  double precision :: 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 dcoo2csr

subroutine zcoo2csr(n,nnz,ic,jc,c,isa,jsa,sa)
  implicit none
  integer :: n,nnz
  integer,dimension(*) :: ic,jc,isa,jsa
  complex(kind=kind(1.0d0)),dimension(*) :: c,sa
!!!
  integer :: k,k1,i,j,idum
  integer,dimension(n) :: iloc
  complex(kind=kind(1.0d0)):: 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 zcoo2csr

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine dsubmat (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
    double precision,dimension(*):: a
    integer,dimension(*) :: ia,ja
    double precision,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 dsubmat
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine zsubmat (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(kind=kind(1.0d0)),dimension(*):: a
    integer,dimension(*) :: ia,ja
     complex(kind=kind(1.0d0)),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 zsubmat
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

max_stars_repo_name: lkusch/Kratos

max_issues_repo_name: lkusch/Kratos

max_forks_repo_name: lkusch/Kratos


text:
! From:
!
!     http://icl.cs.utk.edu/lapack-for-windows/lapack/
!

program main

	integer, parameter :: n=2, lda=2, nrhs=1, ldb=2

	double precision :: a(lda,n), b(ldb,nrhs)

	integer :: ipiv(n), info, i, j

	a(1,1) = 1
	a(1,2) = 2
	a(2,1) = 3
	a(2,2) = 4

	b(1,1) = 5
	b(2,1) = 6

	print *, "Hello world"

	! Solve linear system a*x = b, for x
	call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)

	print *, "info = ", info
	print *, "a (factored) = ", a

	! Expected output: x = -4.0  4.5
	print *, "x = ", b

	print *, "done main"

end program main

max_stars_repo_name: JeffIrwin/lapack-hello-world

max_issues_repo_name: JeffIrwin/lapack-hello-world

max_forks_repo_name: JeffIrwin/lapack-hello-world


text:
c { dg-do compile }
       PRINT 10, 2, 3
10     FORMAT (I1, X, I1) ! { dg-warning "Extension: X descriptor" "Extension: X descriptor"  }
       END

ext: f

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
C Copyright (C) 2002, Carnegie Mellon University and others.
C All Rights Reserved.
C This code is published under the Common Public License.

C $Id: main_win32.f,v 1.1 2002/05/02 18:52:13 andreasw Exp $
      program main
C
C     Main program for simulation
C
C     Authors:       Andreas Waechter    09-26-01
C
      implicit none
      character*40 MODNAM, STUBNAM
C
C     Read model and run name from file 'MODEL.DAT' in current directory
C
      open(1,file='MODEL.DAT',status='old',err = 999)
      read(1,'(a)',err=999) MODNAM
      read(1,'(a)',err=999) STUBNAM
      close(1)

      call FMAIN(MODNAM, STUBNAM)
      stop

 999  write(*,*) 'Error reading file MODEL.DAT'
      stop
      end

max_stars_repo_name: kant/CMMPPT

max_issues_repo_name: kant/CMMPPT

max_forks_repo_name: kant/CMMPPT


! { dg-do run }
!
! Tests the fix for PR79676 in which submod_test was private even to the
! submodule 'my_submod'.
!
! Contributed by Adam Hirst  <adam@aphirst.karoo.co.uk>
!
module my_mod
  private           ! This hid 'submod_test'.
  interface
    module subroutine submod_test(x)
      integer :: x
    end subroutine
  end interface
  integer answer
  public routine1, print_two, answer
contains
  subroutine routine1(x)
    integer :: x
    call submod_test(x)
  end subroutine
  subroutine print_two()
    integer, parameter :: two = 2
    answer = answer * two
  end subroutine
end module

module my_mod_2
  use my_mod
contains
  subroutine circular_dependency()
    call print_two()
  end subroutine
end module

submodule (my_mod) my_submod
  use my_mod_2
contains
module subroutine submod_test(x)
  integer :: x
  answer = x
  call circular_dependency()
end subroutine

end submodule

program hello
  use my_mod
  implicit none
  call routine1(2)
  if (answer .ne. 4) call abort
end program


text:
      SUBROUTINE SWAP
      RETURN
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C $Log: conlab.f,v $
C Revision 1.2  2009/03/25 12:36:43  gdsjaar
C Add copyright and license notice to all files.
C Permission to assert copyright has been granted; blot is now open source, BSD
C
C Revision 1.1  1994/04/07 19:57:10  gdsjaar
C Initial checkin of ACCESS/graphics/blotII2
C
c Revision 1.2  1990/12/14  08:48:58  gdsjaar
c Added RCS Id and Log to all files
c
C=======================================================================
      SUBROUTINE CONLAB (ICNTR, CNTR, NHIT, LABINC, LINSET,
     &   VARNP, XN, YN, ZN, IN2ELB, *)
C=======================================================================

C   --*** CONLAB *** (DETOUR) Label contour line
C   --   Written by Amy Gilkey - revised 06/01/87
C   --   D. P. Flanagan, 03/30/83
C   --
C   --CONLAB labels a contour line that intersects the given line with the
C   --contour letter.  If LABINC is greater than one, only every LABINCth
C   --crossing is labeled.
C   --
C   --Parameters:
C   --   ICNTR - IN - the contour label number
C   --   CNTR - IN - the contour value
C   --   NHIT - IN/OUT - the number of "hits" so far, initialize to zero
C   --   LABINC - IN - the number of "hits" to skip before labelling
C   --   LINSET - IN - the line set
C   --   VARNP - IN - the contour function values
C   --   XN, YN, ZN - IN - the nodal coordinates
C   --   IN2ELB - IN - the element block for each node;
C   --      <0 if not in any selected element block
C   --      =0 if in more than one selected element block
C   --   * - return statement if the cancel function is active
C   --
C   --Common Variables:
C   --   Uses IS3DIM of /D3NUMS/

      PARAMETER (DXO = .0015, DYO = .0015)

      COMMON /D3NUMS/ IS3DIM, NNPSUR, NUMNPF, LLNSET
      LOGICAL IS3DIM

      INTEGER LINSET(3)
      REAL VARNP(*)
      REAL XN(*), YN(*), ZN(*)
      INTEGER IN2ELB(*)

      LOGICAL INTERP_BL
      LOGICAL GRABRT
      LOGICAL EXISTS

      EXISTS (M) = (MOD(M,2) .NE. 0)

      IF (IS3DIM) THEN
C      --Skip invisible line
         IF (LINSET(3) .EQ. 0) GOTO 100
      END IF

      N1 = LINSET(1)
      N2 = LINSET(2)
      IF ((IN2ELB(N1) .GE. 0) .AND. (IN2ELB(N2) .GE. 0)) THEN

         IF (INTERP_BL (CNTR, VARNP(N1), VARNP(N2), PSI)) THEN
            NHIT = NHIT + 1
            IF (NHIT .GE. LABINC) THEN
               IF (GRABRT ()) RETURN 1
               IF (IS3DIM) THEN
C               --Replace invisible node with partial line node
                  IF (ABS (LINSET(3)) .GT. 1)
     &               N2 = ABS (LINSET(3))
               END IF
               X0 = XN(N1) * (1.-PSI) + XN(N2) * PSI
               Y0 = YN(N1) * (1.-PSI) + YN(N2) * PSI
               CALL MP2PT (1, X0, Y0, DX0, DY0, MASK)
               IF (EXISTS (MASK)) CALL PLTXTS (DX0+DXO, DY0+DYO,
     &            CHAR (ICHAR('A')-1+ICNTR))
               NHIT = 0
            END IF
         END IF
      END IF

  100 CONTINUE
      RETURN
      END

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


text:
*----------------------------------------------------------------------
*       S P I N T G R L
*----------------------------------------------------------------------
*       This routine reads the yint.lst data storing the concatenated
*       arrays for the different blocks

	SUBROUTINE SPINTGRL(cf_tot)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      include 'mpif.h'
      parameter (MAXPROC=100)
      common /MPI/ myid, nprocs, ierr, istat(MPI_STATUS_SIZE)
*
        PARAMETER (NWD=70,NOD=220,NOFFD=800)
        PARAMETER(MEIG=20,MTERM=20)
*
        CHARACTER EL*3,ATOM*6,TERM*6
        INTEGER       IN,OUT,PRI,OUC,OUF,OUH,OUD
        COMMON/INOUT/ IN,OUT,PRI,IUC,IUF,IUD,OUC,OUF,OUD,OUH,ISCW
        COMMON/LABEL/ EL(60),ATOM,TERM
        LOGICAL       FAIL,OMIT,EZERO,REL,ALL,TRACE
        COMMON/TEST/  FAIL,OMIT,EZERO,REL,ALL,TRACE
       COMMON/WAVE/EC,ED,AZD,PDE(NOD),IJE(noffd),EIJ(noffd),IPR
        COMMON/PARAM/  H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,MASS,
     :  NCFG,IB,IC,ID,
     :                D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,FINE,
     :                NSCF,NCLOSD,RMASS
        COMMON/RADIAL/R(NOD),RR(NOD),R2(NOD),YK(NOD),YR(NOD),X(NOD)
      POINTER (QWT,WT(1)),(QWP,WP(1)), (IQWPTR,W(1))
      COMMON /CFGS/ETOTAL,QWT,QWP,IQWPTR
*
      POINTER (pkval,kval(1)),(pvalue,value(1)),
     :        (pih,ih(1)),(pjh,jh(1)),
     :        (pcoeff,coeff(1)),(pnijptr,nijptr(1)),
     :        (pinptr,inptr(1))
      COMMON/STATE/INTPTR(6),IDIM,PKVAL,PVALUE,NZE,PIH,PJH,
     :             NCODIM,PCOEFF,PNIJPTR,PINPTR

*     .. Logical array passed from nonh to give the used integrals
      POINTER (plused,lused(1))
      COMMON/USED/plused
      LOGICAL lused

      POINTER (pico,ico(1))
      COMMON/ICOFF/pico

      POINTER (qjptr,jptr(1))
      COMMON/COLS/qjptr

      integer cf_tot(nblock), n_read, n_left;
      logical leigen(meig,mterm), lguess
      integer nume(mterm),iws(mterm),iiws(mterm),
     :        nze_bl(mterm), ncfg_bl(mterm), niv_bl(mterm)
      character term_bl(mterm)*3
      pointer   (qeigvec,eigvec(1)),(pen,en(1))
      common/st/leigen,nblock,nume,iws,iiws,ncfg_bl,nze_bl,
     :          term_bl, qeigvec, pen, lguess
*
      logical                   ::clst_disk, clst_memory
      common/use_disk/clst_disk, clst_memory

      CHARACTER*72 string

*     print *, 'Entering spintgrl'
*      read the begining of yint.lst
       rewind(iud)
       read(iud) ncl, mwf, nbb, lsd
       maxorb = nwf-nclosd
       if (ncl.ne.nclosd.or.mwf.ne.maxorb.or.nbb.ne.nblock) then
         write(0,*) 'Yint.lst data not consistent with cfg.h data'
         write(0,*) 'Yint.lst :', ncl, mwf, nbb
         write(0,*) 'cfg.h    :', nclosd, nwf, nblock
         stop
       end if
*      .. electron labels have already been read.
       read(iud) string

101    read(iud) string
       maxorb = maxorb - 24
       if (maxorb .gt. 0) go to 101
*      end read

       call alcsts(cf_tot);
       ico = 0
       ih = 0;
       jptr = 0;
       kval = 0;

       ib1 = 0
       ib2 = 0
*      read for each block
       do iblock = 1, nblock
         nij = 0
2        read(iud,iostat=ierr) n, (ih(i+nij),i=ib1+1,ib1+n),
     :                (ico(i+nij),i=ib1+1,ib1+n)
         nij = nij+n
         if (n .ge. lsd) go to 2
         read(iud,iostat=ierr) ncol,(jptr(i),i=ib2+1,(ib2+ncol))
       ib1 = ib1 + nij
       ib2 = ib2 + ncol
*       print '(6i12)',(ih(i),i=1,ib1)
*       print '(6i12)',(ico(i),i=1,ib1)
*       print '(6i12)',(jptr(i),i=1,ib2)
       end do

        nij = 0
C
C ***** READ  THE LIST OF INTEGRALS
C
        lasti = 0

*
        DO 10 INT = 1,4
*
*          ...F, G, R, or L integrals....
           icount = 1
           read(iud) itype, nint

*          .. ipackn of nonh is called kval here
*          .. integral pack number and logical variable indicating usage
           read(iud) (kval(i),i=lasti+1,nint),
     :               (lused(i),i=lasti+1,nint)
           lasti = nint
           icount = icount + 1
           intptr(int) = lasti
   10  CONTINUE

*      .. adjust for the fact that overlaps have been removed

       intptr(6) = intptr(4)
       intptr(5) = intptr(3)
       intptr(4) = intptr(2)
       intptr(3) = intptr(2)

*      ibe = 0;
*      n = ncodim;
*      do iblock = 1, nblock
*        do while (n == ncodim)
*        read(39) n, (coeff(j),j = ibe+1,ibe+n),
*     :                    (inptr(j), j = ibe+1,ibe+n)
*          ibe = ibe + n;
*        end do
*      end do

      if (clst_memory) then
        ibe = 0;
        do iblock = 1, nblock
          n_left = cf_tot(iblock)
          if (cf_tot(iblock)>ncodim) then
            n_read = ncodim;
          else
            n_read = cf_tot(iblock);
          end if
          do while (n_left > 0)
            if (n_left<ncodim) n_read = n_left;
            read(39) n, (coeff(j),j = ibe+1,ibe+n_read),
     :                    (inptr(j), j = ibe+1,ibe+n_read)
*          if (myid == 0) then
*            print*, n, (coeff(j),j = ibe+1,ibe+n_read),
*     :                    (inptr(j), j = ibe+1,ibe+n_read)
*            end if
            ibe = ibe + n_read;
            n_left = n_left - n_read
          end do
        end do
      end if
      return
      end


text:
program rsir
use my_init
use random, only: random_Poisson, random_unif, init_random_seed
use mcintegral, only: mc_int
use conf, only: mconf, rtf_n
use quartiles, only: quart

implicit none
!============= Define variables  ==================!
integer :: ni, ns, nr, nl, nd, sic, nod, t_recov, t_laten, &
           ntot, nsap, ro, t_infecc, hop, cop, traj, &
           ini, ins, inr, inl, ind, tsh, erst, opt, wop

real :: rot

integer, allocatable :: vr(:), vl(:), vd(:), vsic(:) ! memory vectors
!real, allocatable :: v_ran(:)     ! vector
integer, allocatable :: vnl(:), vns(:), vni(:), vnr(:), & !cummulant vectors
                        vnp(:), vnd(:)

integer, allocatable :: avnl(:), avns(:), avni(:), avnr(:), &
                        avnp(:), avnd(:), avl(:), avcumv(:) !av. cummulant vectors

real(wp), allocatable:: rt_n(:) !effective rep.num

! alpho:   parameter of Herd function
! hp:      lower limit to achieve stationarity in herd immunity
! gamma:   parameter to simulate gaussian confinement function
! io:      parameter to describe at what stage gaussian confinement begins

real(wp), parameter:: alpho=0.22_wp, hp=0.01_wp, gamma=50.0_wp, io=0.2_wp
real(wp) :: w, hf, cf, leth, un_leth

!dummy variables
integer :: j, i, k, lb, cumvl, fl
real(wp) :: lk, lkk, avarea
real(wp), allocatable :: cumv(:,:), cv(:), area(:), incv(:,:)
real :: start, finish

real(wp) :: wpar(3,10)
! =======================================================!

call cpu_time(start)
call init_random_seed() ! generate random seed

!======== read initial conditions and parameters ========!
open(unit=23, file='par.dat', status='old')
read(23,*) ini, ins, inr, inl, ind, t_recov, t_laten, ntot, &
           ro, hop, leth, cop, traj, tsh, opt, wop
close(unit=23)

open(unit=24,file='wpar.dat',status='old')
do i = 1 , 10
  read(24,*) wpar(3,i), wpar(2,i),  wpar(1,i)
enddo
close(unit=24)

t_infecc = t_recov - t_laten      !infective time
rot= float(ro)/float(t_infecc)    !infections per unit time

write(20 ,'(i8,6ES24.13E3)')  0 , dfloat(inl), dfloat(ins) , dfloat(ini) ,&
                             dfloat(inr), dfloat(ind), dfloat(ini)

! =======================================================!
! one colud in principal introduce fluctuations around rot to simulate
!super or subpar infectious events
! allocate( v_ran( (ini + ins + inr + inl + ind)*t_infecc ) )
! v_ran = rot

! =============================================================================!
!============== begining of trajectories calculation ==========================!
allocate(avnl(ntot + t_laten), avns(ntot), avni(ntot), avnr(ntot + t_recov),&
         avnp(ntot), avnd(ntot + t_laten), avl(ntot + t_laten), &
         cumv(traj,ntot), cv(ntot), area(traj), rt_n(ntot), &
         incv(traj,ntot))

avnl = 0
avns = 0
avni = 0
avnr = 0
avnd = 0
avl = 0

lk = rero
lkk = rero
rt_n = rtf_n(ntot, tsh, t_infecc, t_laten) * t_infecc/dfloat(ro)

do k = 1 , traj

  allocate( vsic(ntot), vl(ntot + t_laten), vr(ntot + t_recov), &
            vd(ntot + t_laten), vnl(ntot + t_laten),vns(ntot), &
            vni(ntot),vnr(ntot + t_recov), vnd(ntot  +t_laten) , vnp(ntot) )
  vsic = 0
  vl = 0
  vr = 0
  vd = 0
  vnl = 0
  vns = 0
  vni = 0
  vnr = 0
  vnd = 0
  vnp = 0

  lb = 1

  vnl(1) = inl
  vns(1) = ins
  vni(1) = ini
  vnr(1) = inr
  vnp(1) = vnl(1) + vns(1) + vni(1) + vnr(1)
  vnd(1) = ind

  nl = 0
  ns = 0
  ni = 0
  nr = 0
  nd = 0

  vr = 0
  vl = 0
  vd = 0
  cumvl = ini
  cf = one

  erst=0

!=================== Single trajectories calculation ==========================!
  do j = 1  , ntot

    if(j==1) then
      vr(t_recov) = vni(j)
      hf = one
      cf = one
    else

!=================== Herd immunity =====================!
      if (hop == 1) then
        hf = one !no contingency
      else if (hop == 2) then
        hf  =  (one + exp((dfloat(vni(1))/dfloat(vnp(1)))/alpho))/&
        (one + exp((dfloat(cumvl)/dfloat(vnp(1)))/alpho)) &
        * (one - hp) + hp
      endif

! ================= Confinement =========================!
      if (cop==1 .and. wop==0) then !no confinement
        cf = one
      else if (cop==2.and.wop==0) then  ! piecewise confinement
        cf = mconf(j,wpar)
      else if (cop == 3.and.wop==0 .and.  & !confinement upon the cumulative of te incidence
            dfloat(vni(j)) / dfloat(vnp(1)) >= io) then
        cf = exp(- gamma**2 * (dfloat(vni(j)) / dfloat(vnp(1)))**2)
      else if (cop == 4.and.wop==1) then  !confinement upon empirical R(t)
        if (ni >= 1000 .and. erst==0) then
          erst=1
          if (rt_n(j) > one) then
           cf = mconf(j,wpar)
          elseif ( rt_n(j) < one ) then
           cf = rt_n(j)
          end if
        else if (erst==1) then
          if (rt_n(j) > one) then
           cf = mconf(j,wpar)
          elseif ( rt_n(j) < one ) then
           cf = rt_n(j)
          end if
        else
          cf = mconf(j,wpar)
        endif
      endif

      if (wop==1) then
        w = cf
      else if (wop==0) then
        w = hf * cf !wise multiplication
      endif

! store weight function
      if (k==1) then
        if (j==1) then
          write(16,'(i8,3ES24.13E3)') j , one , cf , hf
        else
          write(16,'(i8,3ES24.13E3)') j , w , cf , hf
        endif
      else if (traj/=1.and.k==traj/2) then
        if (j==1) then
          write(17,'(i8,3ES24.13E3)') j , one , cf , hf
        else
          write(17,'(i8,3ES24.13E3)') j , w , cf , hf
        endif
      else if (traj>=3.and.k==traj) then
        if (j==1) then
          write(18,'(i8,3ES24.13E3)') j , one , cf , hf
        else
          write(18,'(i8,3ES24.13E3)') j , w , cf , hf
        endif
      endif

! =============== Incidence generation ==================!
      if (ns<=0) then
        sic = 0
        w = rero
      else ! ----------- > total number susceptible in contagious domains
        sic =  0
        do i = 1, ni
          ! nsap = random_Poisson( v_ran(lb) * real(w) ,.true.)
            nsap = random_Poisson( rot * real(w) ,.true.)
          lb = lb + 1
          sic = sic + nsap
        ! we give here certain conditions to avoid having negative numbers
          if (ns - sic <= 0) then ! --- >  last number of sic is remaining ns
            ns = sic
            exit
          endif
        enddo
      endif

      vsic(j) = sic
      vl(j + t_laten) = vsic(j)
      vr(j + t_recov) = vsic(j)

    endif

! ===================== Deads ===========================!

    if (j==1) then
      nod = 0
      do i = 1 , vni(1)
        un_leth = random_unif(rero , one)
        if (un_leth > one - leth) then
          nod = nod + 1
        end if
      enddo
    else if (j>= t_laten + 1.and.vl(j) /= 0) then
      nod = 0
      do i = 1 , vl(j)
        un_leth = random_unif(rero , one)
        if (un_leth > one - leth) then
          nod = nod + 1
        end if
      enddo
    else
      un_leth = rero
      nod = 0
    endif

    vd(j+ t_laten ) = nod

! =================== Cummulative =======================!
    cumvl = cumvl + vl(j)
    cumv(k,j) = cumvl  ! write cumulant vector
    incv(k,j) = dfloat(vl(j))

! ================ SLIRD evolution =======================!
    vnl(j+1) = vnl(j) + vsic(j) - vl(j)
    vns(j+1) = vns(j) - vsic(j)
    vni(j+1) = vni(j) + vl(j) - vr(j)
    vnr(j+1) = vnr(j) + vr(j) - vd(j)
    vnd(j+1) = vnd(j) + vd(j)
    vnp(j+1) = vnl(j+1) + vns(j+1) + vni(j+1) + vnr(j+1)

    nl = vnl(j+1)
    ns = vns(j+1)
    ni = vni(j+1)
    nr = vnr(j+1)
    nd = vnd(j+1)

    if (ns<=0) then
      ns = 0
    endif

! ============ > write single trajectories
    if (k==1) then
      write(20 ,'(i8,7ES24.13E3)')  j , dfloat(nl), dfloat(ns) , dfloat(ni) ,&
      dfloat(nr), dfloat(nd), dfloat(cumvl)
      write(25,'(7i10)') j , vsic(j), vl(j), vr(j), vni(j), vd(j), nod
    endif

    if (traj>1.and.k==traj/2) then
      write(21 ,'(i8,7ES24.13E3)')  j , dfloat(nl), dfloat(ns) , dfloat(ni) ,&
      dfloat(nr), dfloat(nd), dfloat(cumvl)
      write(26,'(7i10)') j , vsic(j), vl(j), vr(j), vni(j), vd(j), nod
    endif

! =========== > Ommit non-propagated events
    if (j == ntot.and.ns/=ins) then
      lkk = lkk + one
      fl=1
    else if (j>1.and.ns==ins) then
      fl=0
    endif

  enddo !------------- > End of single trajectory calculation

!=============================================================================!

  if (fl==1) then
    avnl = avnl + vnl
    avns = avns + vns
    avni = avni + vni
    avnr = avnr + vnr
    avnd = avnd + vnd
    avnp = avnp + vnp
    avl = avl + vl
  endif

!get minimum and maximum areas under the cummulative curves
!we compare the area of single trajectories and store only the
!smallets and the greatest

  cv = cumv(k,:)
  if (fl == 1) then
   area (k) = mc_int(cv,ntot)
  endif

  deallocate( vsic, vl, vr, vd, vnl, vns, vni, vnr, vnd, vnp )
  lk = lk + one

print*, int(lk), int(lkk)

end do !------------- > End of trajectories calculation
! =============================================================================!

!===================== Averaging over trajectories ============================!
cumvl = ini

if (traj==1) then

  do i = 1 , ntot
    cumvl = cumvl + avl(i)

    write(15, '(i8,7ES24.13E3)') i, dfloat(avl(i)), dfloat(avnl(i)), &
                                 dfloat(avns(i)), dfloat(avni(i)),&
                                 dfloat(avnr(i)), dfloat(avnd(i)), &
                                 dfloat(cumvl)
  enddo

else

  do i = 1 , ntot

    cumvl = cumvl + avl(i)

    if (i == 1) then
      write(15, '(i8,12ES24.13E3)') i, dfloat(avl(i))/lkk, dfloat(inl), &
                         dfloat(ins), dfloat(ini), dfloat(inr), dfloat(ind), &
                         dfloat(ini), cumv(minloc(area),1),cumv(maxloc(area),1), &
                         incv(minloc(area),1), incv(maxloc(area),1)
    else
      write(15, '(i8,12ES24.13E3)') i, dfloat(avl(i))/lkk, dfloat(avnl(i))/lkk, &
                    dfloat(avns(i))/lkk, dfloat(avni(i))/lkk,&
                    dfloat(avnr(i))/lkk, dfloat(avnd(i))/lkk, &
                    dfloat(cumvl)/lkk, cumv(minloc(area),i),cumv(maxloc(area),i),&
                    incv(minloc(area),i), incv(maxloc(area),i)
    endif

  enddo
endif

! compute quartiles
call quart(cumv, traj, ntot)

! =============================================================================!
!write all trajectories
if ( opt == 0 ) then
  do k=1 , traj
    do i = 1 , ntot
      write(100,'(i8,2ES24.13E3)') i, incv(k,i), cumv(k,i)
    enddo
    write(100,*)
  enddo
endif

call cpu_time(finish)
print*, start/60., finish/60., (finish-start)/60.

end program

max_stars_repo_name: RenatoSalArrDu/StochasticSLIRD

max_issues_repo_name: RenatoSalArrDu/StochasticSLIRD

max_forks_repo_name: RenatoSalArrDu/StochasticSLIRD


text: #redirect Sugar Daddies

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE SMUTEX(BUF,LBUF,IBUF,SCR,LSCR)
!     SMUTEX IS THE EXECUTION PHASE OF THE SEISMIC PROCESS MUTE, WHICH MUTES OR
!  ZEROS A PORTION OF A TRACE.  THE MUTE PARAMETERS MUST BE ON DISC FILE MUNIT
!  AND THE TRACE (WITH HEADER) MUST BE IN MEMORY LOCATION BUF.  THE MUTE TIMES
!  FOR TRACES BETWEEN THOSE SHOTS OR RPS DESCRIBED BY THE USER ARE CALCULATED
!  BY LINEAR INTERPOLATION.  LIKEWISE, TIMES FOR TRACES NOT SPECIFIED BY THE
!  USER OBTAINED THROUGH EXTRAPOLATION OR INTERPOLATION.  DEEP WATER DELAYS ARE
!  HONORED - I.E. MUTE TIMES ARE RELATIVE TO TIME ZERO REGARDLESS OF THE TIME
!  OF THE FIRST DATA SAMPLE.
!     SUBROUTINE MUTEED CONTAINS THE EXPLAINATION OF THE USER PARAMETERS AND THE
!  ORDER OF THE USER PARAMETERS ON DISC.
!
!  ARGUMENTS:
!  BUF    - THE TRACE TO BE MUTED, INCLUDING THE TRACE HEADER.  THE FIRST
!           DATA SAMPLE MUST BE AT TIME DELAY.  THIS IS THE FLOATING
!           POINT (REAL) TRACE ARRAY.
!  LBUF   - THE LONG INTEGER TRACE ARRAY.  THIS IS REALLY THE SAME AS BUF, BUT
!           PRIME FORTRAN DOESN'T ALLOW EQUIVALENCING ANYTHING TO AN ARGUMENT.
!  IBUF   - THE SHORT INTEGER TRACE ARRAY.  NEEDED FOR 16 BIT TRACE HEADER
!           ADDRESSES.
!  SCR    - A SCRATCH ARRAY FOR READING THE PARAMETERS.  THEREFORE, SCR MUST
!           BE AT LEAST 56 32BIT WORDS BIG.  SCR MAY BE DESTROYED BY THE CALLING
!           ROUTINE.
!  LSCR   - THE SAME SCRATCH ARRAY BECAUSE OF THE EQUIVALENCING PROBLEM.
!
!
!  WRITTEN AND COPYRIGHTED BY:
!  PAUL HENKART, SCRIPPS INSTITUTION OF OCEANOGRAPHY, 3 MAY 1982
!  ALL RIGHTS ARE RESERVED BY THE AUTHOR.  PERMISSION TO COPY OR REPRODUCE THIS
!  SUBROUTINE, BY COMPUTER OR OTHER MEANS, MAY BE OBTAINED ONLY FROM THE AUTHOR.
!
!
!  10 Dec 98 - Bizarre situation where default time of 99999 and a VERY
!              small sample interval converted to a negative index. geez
!  23 Aug 91 - make it work (smuted was changed sometime or other and it broke smutex)
!                   - make maxset = 300
!                   - make interp work
!  27 Aug 95 - set stime to delay if delay > stime
!  March 1999 - set the SEG-Y trace header start/end times in mils
!               (short words 56 & 57)
!  21 March 00 - Spatial interpolation was incorrect.
!  4 Oct 02 - Allow fno 0
!  20 Apr 07 - Make iaddwb a multiplier
!  12 Apr 10 - Eliminate writing the mute times to the segy header.
!
      PARAMETER (MAXSET = 300)                                          ! THE MAXIMUM NUMBER OF ELEMENTS THE USER ARRAY CAN BE
      PARAMETER (NWRDS=MAXSET+9)                                        ! THE NUMBER OF WORDS IN EACH DISC PAR LIST
      DIMENSION BUF(111),LBUF(111),IBUF(111),SCR(1111),LSCR(1111)
      INTEGER*2 IBUF
      DIMENSION CURSET(MAXSET)
      COMMON /SMUTER/ MUNIT,NLISTS
      COMMON /SIOAP/ IASGND,IRELSE,IN,IOUT,NEXTAD,LAPSIZ,IFREE,IUSEAP
      COMMON /APMEM/A(32766)
      COMMON /READT/ ILUN,NUMHDR
      INTEGER FNO, fno2
      REAL newset(maxset), oldset(maxset)
      LOGICAL FIRST
      SAVE
      DATA FIRST /.TRUE./, fno2/0/
!****
!****     FIND THE PARAMETER LIST (ON DISC) FOR THIS SHOT (RP)
!****
      IF(IBUF(15).EQ.2) RETURN                                          ! IS IT A DEAD TRACE
      ISIG=0
   10 IF( first ) THEN
          first = .FALSE.
          mlists = 1
          CALL podisc( munit, 1, 0 )                                    ! rewind
          CALL rddisc( munit, lscr, nwrds, istat )                      ! read the first parameter list
          fno = lscr(1)
          lno = lscr(2)
          iaddwb = lscr(3)
          nsets = lscr(4)
          ltype = lscr(5)
          lprint = lscr(6)
          interp = lscr(7)
          DO 15 i = 1, nsets
             newset(i) = scr(i+9)
             curset(i) = scr(i+9)
   15     CONTINUE
      ENDIF
      IF( lbuf(7) .EQ. 0 ) THEN
          lnum = lbuf(3)                                                ! shot sorted
          trno = lbuf(4)
      ELSE
          lnum = lbuf(6)                                                ! rp sorted
          trno = lbuf(7)
      ENDIF
   20 IF( lnum .LT. fno ) THEN                                          ! assume shot numbers always increase
          IF( mlists .EQ. 1 ) THEN
              IF( interp .EQ. 0 ) RETURN
              GOTO 500
          ENDIF
      ENDIF
      IF( lnum .LE. lno ) GOTO 500                                      ! is the shot in this list
      IF( mlists .GT. nlists ) THEN
          IF( interp .EQ. 0 ) RETURN
          GOTO 500
      ENDIF
      IF( lnum .GE. fno2 ) THEN
          IF( mlists .GT. 1 ) THEN
              fno = fno2
              lno = lno2
              nsets = nsets2
          ENDIF
          DO 90 i = 1, nsets
             curset(i) = newset(i)
             oldset(i) = newset(i)
   90     CONTINUE
          IF( mlists .LT. nlists ) THEN
              CALL rddisc( munit, lscr, nwrds, istat )                          ! read the next parameter list
              fno2 = lscr(1)
              lno2 = lscr(2)
              nsets2 = lscr(4)
              ltype = lscr(5)
              DO 100 i = 1, maxset
                 oldset(i) = newset(i)
                 newset(i) = scr(i+9)
  100         CONTINUE
          ENDIF
          mlists = mlists + 1
          GOTO 20
      ENDIF
      IF( interp .EQ. 0 ) RETURN
!****
!****   Do the spatial interpolation
!****
      IF( nsets .NE. nsets2 ) THEN
          PRINT *,' ***  ERROR  ***  The number of mutes must be the ',
     &      'same in order to perform spatial variation.'
          STOP
      ENDIF
      DO 300 i = 1, nsets, 3
!         IF( oldset(i) .NE. newset(i) ) THEN
!             PRINT *,' ***  ERROR  ***  The trace numbers or ranges ',
!     &         ' must be the same in each list when spatially varying.'
!             STOP
!         ENDIF
!         curset(i) = oldset(i)
         temp =  (FLOAT(lnum) -FLOAT(lno)) / (FLOAT(fno2)-FLOAT(lno))
         curset(i) = temp * (newset(i) - oldset(i)) + oldset(i)
         curset(i+1) = temp * (newset(i+1) - oldset(i+1)) + oldset(i+1)
         curset(i+2) = temp * (newset(i+2) - oldset(i+2)) + oldset(i+2)
  300 CONTINUE
!****
!****    NOW FIND THE MUTE TIME FOR THE  TRACE NUMBER (OR RANGE)
!****
  500 I=1
      IF( LTYPE .EQ. 1 ) THEN
          TEMP=LBUF(10)                                                 ! TRNO=ABS(LBUF(10) DOESN'T WORK!!!
          TRNO=ABS(TEMP)                                                ! SET THE TRACE NO TO THE RANGE IF XTP WAS GIVEN
      ENDIF
  510 IF( trno .LT. curset(i) ) THEN
          IF( interp .EQ. 0 ) RETURN
          time1 = curset(i+1)
          time2 = curset(i+2)
          GOTO 2000
      ENDIF
      IF( trno .EQ. curset(i) ) THEN
          time1 = curset(i+1)
          time2 = curset(i+2)
          GOTO 2000
      ENDIF
      IF( i+3 .GE. nsets ) THEN                                         ! any more smuts?
          IF( interp .EQ. 0 ) RETURN                                    ! no more smut
          time1 = curset(i+1)
          time2 = curset(i+2)
          GOTO 2000
      ENDIF
      IF( trno .GE. curset(i+3) ) THEN
          i = i + 3
          GOTO 510
      ENDIF
      IF( interp .EQ. 0 ) RETURN                                        ! it's bigger than the last given and no extrapolation
      temp = (trno - curset(i)) / (curset(i+3) - curset(i) )
      time1 = temp * (curset(i+4) - curset(i+1)) + curset(i+1)
      time2 = temp * (curset(i+5) - curset(i+2)) + curset(i+2)
!****
!****   Apply the mute !
!****
 2000 CONTINUE
      SI=BUF(49)                                                        ! THE SAMPLE INTERVAL IN SECONDS
      DELAY=BUF(46)                                                     ! THE DEEP WATER DELAY IN SECONDS
      IF( IADDWB .GT. 0 ) THEN                                          ! ADD IN THE WATER BOTTOM TIME
          time1 = time1 + buf(50) * iaddwb
          time2 = time2 + buf(50) * iaddwb
      ENDIF
      IS = NINT((TIME1-DELAY)/SI) + 1
      IF( is .LT. 1 ) is = 1                                            ! watch out for delay > time1
      IE = NINT((TIME2-DELAY)/SI) + 1                                   ! THE END INDEX OF THE MUTE
      NSAMPS=IBUF(58)                                                   ! 16 BIT NUMBER OF SAMPLES PER TRACE
      IF( TIME2 .EQ. 99999. ) IE = NSAMPS
      IF(IE.LT.IS) IE=IS
      IF(IS.GE.IE) RETURN                                               ! DON'T DO ANY MUTING IF NONE TO DO!
      IF(IE.GT.NSAMPS) IE=NSAMPS
!      IF(IN.NE.0.AND.IUSEAP.NE.0) GO TO 2100                            ! DATA IS IN THE AP IF IN<>0
      IF(IASGND.NE.0.AND.IUSEAP.NE.0) GO TO 2100
      IF(IUSEAP.EQ.0.AND.IN.NE.0) GO TO 2020
      CALL MUTE(BUF(numhdr+1),IS,IE,SI,NSAMPS)                             ! THE IN MEMORY MUTE
      GO TO 2200
 2020 CALL MUTE(A(IN),IS,IE,SI,NSAMPS)
      GO TO 2200
 2100 CONTINUE                                                          ! THE TRACE BELONGS IN THE AP
      CALL INAP(BUF(NUMHDR+1),NSAMPS)                                   ! PUT THE DATA IN THE AP
      IF(IS.NE.1.OR.IE.LT.NSAMPS) GO TO 2150                            ! ARE WE KILLING THE WHOLE TRACE?
      CALL VCLR(IN,1,NSAMPS)                                            ! ZERO THE WHOLE TRACE
      IBUF(15)=2                                                        ! SET THE DEAD TRACE FLAG
      GO TO 2200
 2150 N=5                                                               !  THE LENGTH OF THE TAPER
      CALL APWR                                                         ! MAKE SURE WE DON'T CLOBBER ANYTHING IN PROGRESS IN NEXTAD
      TEMP=.2                                                           ! THE INCREMENT OF THE TAPER BETWEEN POINTS OF THE TAPER
      IF(IE.LE.N) N=0
      IF(IE-IS.LT.N+N) N=0                                              ! WATCH OUT FOR VERY SHORT MUTES
      CALL APPUT(TEMP,NEXTAD,1,2)                                       ! PUT THE INCREMENT IN AP LOCATION NEXTAD
      IIS=IS+IN-1                                                       ! ADD IN THE AP ADDRESS OF THE BEGINNING OF THE TRACE
      IIE=IE+IN-1
      CALL APWD                                                         ! WAIT FOR THE AP DATA TRANSFER TO COMPLETE
!      CALL MUTEAP(IIS,IIE,N,NEXTAD,NSAMPS)
 2200 CONTINUE                                                          !  NOW PUT THE MUTE TIMES IN THE TRACE HEADER
      IF(IS.EQ.1.AND.IE.GE.NSAMPS) IBUF(15)=2                           ! SET THE DEAD TRACE FLAG
!      BUF(52)=(IS-1)*SI+DELAY
!      ibuf(56) = NINT(buf(52)*1000.)                                    ! mute start time in mils
!      BUF(53)=(IE-1)*SI+DELAY
!      ibuf(57) = NINT(buf(53)*1000.)                                    ! mute end time in mils
      IF(IAND(LPRINT,2).EQ.2) PRINT 2210,LNUM,TRNO,(IS-1)*SI+DELAY,
     *   (IE-1)*SI+DELAY,IS,IE
 2210 FORMAT(' SHOT/RP ',I6,' TRACE ',F6.1,' IS MUTED FROM ',F10.4,
     &   ' TO ',F10.4,' (',2I8,')')
      RETURN
      end


text:
      MODULE vmec_persistent
      USE stel_kinds, ONLY: rprec1 => rprec
      IMPLICIT NONE
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER, DIMENSION(:), ALLOCATABLE :: ixm, jmin3
      REAL(rprec1), DIMENSION(:,:), ALLOCATABLE :: cosmu, sinmu,
     1   cosmum, sinmum, cosmui, cosmumi, sinmui, sinmumi,
     2   cosnv, sinnv, cosnvn, sinnvn
      REAL(rprec1), DIMENSION(:), ALLOCATABLE ::
     1      xm, xn, xm_nyq, xn_nyq, cos01, sin01
c-----------------------------------------------
      END MODULE vmec_persistent