psblas3/base/modules/serial/psb_serial_mod.f90

!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006-2018
!        Salvatore Filippone
!        Alfredo Buttari
!
!    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. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
module psb_serial_mod
  use psb_const_mod
  use psb_error_mod
  use psb_realloc_mod
  use psb_string_mod
  use psb_sort_mod

  use psi_serial_mod, &
       & psb_gth => psi_gth,&
       & psb_sct => psi_sct

  use psb_s_serial_mod
  use psb_d_serial_mod
  use psb_c_serial_mod
  use psb_z_serial_mod

  interface psb_nrm1
    module procedure psb_snrm1, psb_dnrm1, psb_cnrm1, psb_znrm1
  end interface psb_nrm1

  interface psb_minreal
    module procedure psb_sminreal, psb_dminreal, psb_cminreal, psb_zminreal
  end interface psb_minreal


  interface psb_nspaxpby
    subroutine psb_d_nspaxpby(nz,iz,z,alpha, nx, ix, x, beta, ny,iy,y,info)
      use psb_const_mod, only : psb_ipk_, psb_spk_, psb_dpk_
      integer(psb_ipk_), intent(out)              :: nz
      integer(psb_ipk_), intent(out)              :: iz(:)
      real(psb_dpk_), intent (out)      :: z(:)
      integer(psb_ipk_), intent(in)               :: nx, ny
      integer(psb_ipk_), intent(in)               :: ix(:), iy(:)
      real(psb_dpk_), intent (in)       :: x(:), y(:)
      real(psb_dpk_), intent (in)       :: alpha, beta
      integer(psb_ipk_), intent(out)              :: info
    end subroutine psb_d_nspaxpby
    subroutine psb_s_nspaxpby(nz,iz,z,alpha, nx, ix, x, beta, ny,iy,y,info)
      use psb_const_mod
      integer(psb_ipk_), intent(out)              :: nz
      integer(psb_ipk_), intent(out)              :: iz(:)
      real(psb_spk_), intent (out)      :: z(:)
      integer(psb_ipk_), intent(in)               :: nx, ny
      integer(psb_ipk_), intent(in)               :: ix(:), iy(:)
      real(psb_spk_), intent (in)       :: x(:), y(:)
      real(psb_spk_), intent (in)       :: alpha, beta
      integer(psb_ipk_), intent(out)              :: info
    end subroutine psb_s_nspaxpby
    subroutine psb_c_nspaxpby(nz,iz,z,alpha, nx, ix, x, beta, ny,iy,y,info)
      use psb_const_mod
      integer(psb_ipk_), intent(out)              :: nz
      integer(psb_ipk_), intent(out)              :: iz(:)
      complex(psb_spk_), intent (out)      :: z(:)
      integer(psb_ipk_), intent(in)               :: nx, ny
      integer(psb_ipk_), intent(in)               :: ix(:), iy(:)
      complex(psb_spk_), intent (in)       :: x(:), y(:)
      complex(psb_spk_), intent (in)       :: alpha, beta
      integer(psb_ipk_), intent(out)              :: info
    end subroutine psb_c_nspaxpby
    subroutine psb_z_nspaxpby(nz,iz,z,alpha, nx, ix, x, beta, ny,iy,y,info)
      use psb_const_mod
      integer(psb_ipk_), intent(out)              :: nz
      integer(psb_ipk_), intent(out)              :: iz(:)
      complex(psb_dpk_), intent (out)      :: z(:)
      integer(psb_ipk_), intent(in)               :: nx, ny
      integer(psb_ipk_), intent(in)               :: ix(:), iy(:)
      complex(psb_dpk_), intent (in)       :: x(:), y(:)
      complex(psb_dpk_), intent (in)       :: alpha, beta
      integer(psb_ipk_), intent(out)              :: info
    end subroutine psb_z_nspaxpby
  end interface psb_nspaxpby

  interface
    subroutine symbmm (n, m, l, ia, ja, diaga, &
         & ib, jb, diagb, ic, jc, diagc, index)
      import :: psb_ipk_
      integer(psb_ipk_) :: n,m,l,  ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
           & diagc,  index(*)
      integer(psb_ipk_), allocatable :: ic(:),jc(:)
    end subroutine symbmm
    subroutine lsymbmm (n, m, l, ia, ja, diaga, &
         & ib, jb, diagb, ic, jc, diagc, index)
      import :: psb_ipk_, psb_lpk_
      integer(psb_lpk_) :: n,m,l,  ia(*), ja(*), diaga, ib(*), jb(*), diagb,&
           & diagc,  index(*)
      integer(psb_lpk_), allocatable :: ic(:),jc(:)
    end subroutine lsymbmm
  end interface


contains

  elemental function psb_snrm1(x) result(res)
    real(psb_spk_), intent(in)  :: x
    real(psb_spk_)              :: res
    res = abs( x )
  end function psb_snrm1

  elemental function psb_dnrm1(x) result(res)
    real(psb_dpk_), intent(in) :: x
    real(psb_dpk_)             :: res
    res = abs(  x )
  end function psb_dnrm1

  elemental function psb_cnrm1(x) result(res)
    complex(psb_spk_), intent(in)  :: x
    real(psb_spk_)                 :: res
    res = abs( real( x ) ) + abs( aimag( x ) )
  end function psb_cnrm1

  elemental function psb_znrm1(x) result(res)
    complex(psb_dpk_), intent(in)  :: x
    real(psb_dpk_)                 :: res
    res = abs( real( x ) ) + abs( aimag( x ) )
  end function psb_znrm1

  elemental function psb_sminreal(x) result(res)
    real(psb_spk_), intent(in)  :: x
    real(psb_spk_)              :: res
    res = ( x )
  end function psb_sminreal

  elemental function psb_dminreal(x) result(res)
    real(psb_dpk_), intent(in) :: x
    real(psb_dpk_)             :: res
    res = (  x )
  end function psb_dminreal

  elemental function psb_cminreal(x) result(res)
    complex(psb_spk_), intent(in)  :: x
    real(psb_spk_)                 :: res
    res = min( real( x ) , aimag( x ) )
  end function psb_cminreal

  elemental function psb_zminreal(x) result(res)
    complex(psb_dpk_), intent(in)  :: x
    real(psb_dpk_)                 :: res
    res = min( real( x ) , aimag( x ) )
  end function psb_zminreal


  subroutine crot( n, cx, incx, cy, incy, c, s )
    !
    !  -- lapack auxiliary routine (version 3.0) --
    !     univ. of tennessee, univ. of california berkeley, nag ltd.,
    !     courant institute, argonne national lab, and rice university
    !     october 31, 1992
    !
    !     .. scalar arguments ..
    integer(psb_mpk_) :: incx, incy, n
    real(psb_spk_)    c
    complex(psb_spk_)   s
    !     ..
    !     .. array arguments ..
    complex(psb_spk_) cx( * ), cy( * )
    !     ..
    !
    !  purpose
    !  == = ====
    !
    !  zrot   applies a plane rotation, where the cos (c) is real and the
    !  sin (s) is complex, and the vectors cx and cy are complex.
    !
    !  arguments
    !  == = ======
    !
    !  n       (input) integer
    !          the number of elements in the vectors cx and cy.
    !
    !  cx      (input/output) complex*16 array, dimension (n)
    !          on input, the vector x.
    !          on output, cx is overwritten with c*x + s*y.
    !
    !  incx    (input) integer
    !          the increment between successive values of cy.  incx <> 0.
    !
    !  cy      (input/output) complex*16 array, dimension (n)
    !          on input, the vector y.
    !          on output, cy is overwritten with -conjg(s)*x + c*y.
    !
    !  incy    (input) integer
    !          the increment between successive values of cy.  incx <> 0.
    !
    !  c       (input) double precision
    !  s       (input) complex*16
    !          c and s define a rotation
    !             [  c          s  ]
    !             [ -conjg(s)   c  ]
    !          where c*c + s*conjg(s) = 1.0.
    !
    ! == = ==================================================================
    !
    !     .. local scalars ..
    integer(psb_mpk_) :: i, ix, iy
    complex(psb_spk_)         stemp
    !     ..
    !     .. intrinsic functions ..
    !     ..
    !     .. executable statements ..
    !
    if( n <= 0 ) return
    if( incx == 1 .and. incy == 1 ) then
      !
      !     code for both increments equal to 1
      !
      do  i = 1, n
        stemp = c*cx(i) + s*cy(i)
        cy(i) = c*cy(i) - conjg(s)*cx(i)
        cx(i) = stemp
      end do
    else
      !
      !     code for unequal increments or equal increments not equal to 1
      !
      ix = 1
      iy = 1
      if( incx < 0 )ix = ( -n+1 )*incx + 1
      if( incy < 0 )iy = ( -n+1 )*incy + 1
      do  i = 1, n
        stemp  = c*cx(ix) + s*cy(iy)
        cy(iy) = c*cy(iy) - conjg(s)*cx(ix)
        cx(ix) = stemp
        ix = ix + incx
        iy = iy + incy
      end do
    end if
    return
  end subroutine crot
  !
  !
  subroutine crotg(ca,cb,c,s)
    complex(psb_spk_) ca,cb,s
    real(psb_spk_) c
    real(psb_spk_) norm,scale
    complex(psb_spk_) alpha
    !
    if (cabs(ca) == 0.0) then
      !
      c = 0.0d0
      s = (1.0,0.0)
      ca = cb
      return
    end if
    !

    scale = cabs(ca) + cabs(cb)
    norm = scale*sqrt((cabs(ca/cmplx(scale,0.0)))**2 +&
         &   (cabs(cb/cmplx(scale,0.0)))**2)
    alpha = ca /cabs(ca)
    c = cabs(ca) / norm
    s = alpha * conjg(cb) / norm
    ca = alpha * norm
    !

    return
  end subroutine crotg



  subroutine zrot( n, cx, incx, cy, incy, c, s )
    !
    !  -- lapack auxiliary routine (version 3.0) --
    !     univ. of tennessee, univ. of california berkeley, nag ltd.,
    !     courant institute, argonne national lab, and rice university
    !     october 31, 1992
    !
    !     .. scalar arguments ..
    integer(psb_mpk_) :: incx, incy, n
    real(psb_dpk_)    c
    complex(psb_dpk_)   s
    !     ..
    !     .. array arguments ..
    complex(psb_dpk_) cx( * ), cy( * )
    !     ..
    !
    !  purpose
    !  == = ====
    !
    !  zrot   applies a plane rotation, where the cos (c) is real and the
    !  sin (s) is complex, and the vectors cx and cy are complex.
    !
    !  arguments
    !  == = ======
    !
    !  n       (input) integer
    !          the number of elements in the vectors cx and cy.
    !
    !  cx      (input/output) complex*16 array, dimension (n)
    !          on input, the vector x.
    !          on output, cx is overwritten with c*x + s*y.
    !
    !  incx    (input) integer
    !          the increment between successive values of cy.  incx <> 0.
    !
    !  cy      (input/output) complex*16 array, dimension (n)
    !          on input, the vector y.
    !          on output, cy is overwritten with -conjg(s)*x + c*y.
    !
    !  incy    (input) integer
    !          the increment between successive values of cy.  incx <> 0.
    !
    !  c       (input) double precision
    !  s       (input) complex*16
    !          c and s define a rotation
    !             [  c          s  ]
    !             [ -conjg(s)   c  ]
    !          where c*c + s*conjg(s) = 1.0.
    !
    ! == = ==================================================================
    !
    !     .. local scalars ..
    integer(psb_mpk_) :: i, ix, iy
    complex(psb_dpk_)         stemp
    !     ..
    !     .. intrinsic functions ..
    intrinsic          dconjg
    !     ..
    !     .. executable statements ..
    !
    if( n <= 0 ) return
    if( incx == 1 .and. incy == 1 ) then
      !
      !     code for both increments equal to 1
      !
      do  i = 1, n
        stemp = c*cx(i) + s*cy(i)
        cy(i) = c*cy(i) - dconjg(s)*cx(i)
        cx(i) = stemp
      end do
    else
      !
      !     code for unequal increments or equal increments not equal to 1
      !
      ix = 1
      iy = 1
      if( incx < 0 )ix = ( -n+1 )*incx + 1
      if( incy < 0 )iy = ( -n+1 )*incy + 1
      do  i = 1, n
        stemp  = c*cx(ix) + s*cy(iy)
        cy(iy) = c*cy(iy) - dconjg(s)*cx(ix)
        cx(ix) = stemp
        ix = ix + incx
        iy = iy + incy
      end do
    end if
    return
  end subroutine zrot
  !
  !
  subroutine zrotg(ca,cb,c,s)
    complex(psb_dpk_) ca,cb,s
    real(psb_dpk_) c
    real(psb_dpk_) norm,scale
    complex(psb_dpk_) alpha
    !
    if (cdabs(ca) == 0.0d0) then
      !
      c = 0.0d0
      s = (1.0d0,0.0d0)
      ca = cb
      return
    end if
    !

    scale = cdabs(ca) + cdabs(cb)
    norm = scale*dsqrt((cdabs(ca/cmplx(scale,0.0d0,kind=psb_dpk_)))**2 +&
         &   (cdabs(cb/cmplx(scale,0.0d0,kind=psb_dpk_)))**2)
    alpha = ca /cdabs(ca)
    c = cdabs(ca) / norm
    s = alpha * conjg(cb) / norm
    ca = alpha * norm
    !

    return
  end subroutine zrotg


end module psb_serial_mod