base/Makefile

base/Makefile

@@ -6,6 +6,7 @@ LIBNAME=$(BASELIBNAME)
 LIBMOD=psb_base_mod$(.mod)
 lib:
 	(cd modules; make lib LIBNAME=$(BASELIBNAME) F90=$(MPF90) F90COPT="$(F90COPT) $(MPI_OPT)")
+	(cd newserial; make lib LIBNAME=$(BASELIBNAME) )
 	(cd comm; make lib LIBNAME=$(BASELIBNAME))
 	(cd internals; make lib LIBNAME=$(BASELIBNAME))
 	(cd tools; make lib LIBNAME=$(BASELIBNAME))

base/modules/Makefile

@@ -41,6 +41,7 @@ psb_hash_mod.o: psb_const_mod.o psb_realloc_mod.o
 
 psb_base_mod.o: $(MODULES)
 
+
 blacsmod:
 	(make psb_blacs_mod.o F90COPT="$(F90COPT) $(EXTRA_OPT)")
 

base/modules/psb_error_mod.F90

@@ -471,6 +471,9 @@ contains
       write (0,'("Argument M is required when argument PARTS is specified")')
     case(600)
       write (0,'("Sparse Matrix and descriptors are in an invalid state for this subroutine call: ",i0)')i_e_d(1)
+    case(700)
+      write (0,'("Base version has been called: the actual derived type is incomplete!")')
+
     case (1122)
       write (0,'("Invalid state for communication descriptor")')
     case (1123)

base/newserial/Makefile

@@ -0,0 +1,25 @@
+include ../../Make.inc
+
+MODULES = psbn_base_mat_mod.o psbn_coo_mat.o psbn_csr_mat.o
+
+
+LIBMOD=
+OBJS = 
+LIBDIR=..
+CINCLUDES=-I. 
+FINCLUDES=$(FMFLAG)$(LIBDIR) $(FMFLAG). $(FIFLAG).
+
+
+lib:  $(MODULES)  $(OBJS)  $(LIBMOD)
+	$(AR) $(LIBDIR)/$(LIBNAME) $(MODULES) $(OBJS) $(MPFOBJS)
+	$(RANLIB) $(LIBDIR)/$(LIBNAME) #	/bin/cp -p $(LIBMOD) $(LIBDIR)
+	/bin/cp -p *$(.mod) $(LIBDIR)
+
+psbn_coo_mat.o psbn_csr_mat.o:	psbn_base_mat_mod.o
+
+
+clean:
+	/bin/rm -f $(MODULES) $(OBJS) $(MPFOBJS) *$(.mod)
+
+veryclean: clean
+

base/newserial/psbn_base_mat_mod.f03

@@ -0,0 +1,457 @@
+module psbn_base_mat_mod
+
+  integer, parameter :: psbn_spmat_null_=0, psbn_spmat_bld_=1
+  integer, parameter :: psbn_spmat_asb_=2, psbn_spmat_upd_=4
+
+  integer, parameter :: psbn_ireg_flgs_=10, psbn_ip2_=0
+  integer, parameter :: psbn_iflag_=2, psbn_ichk_=3
+  integer, parameter :: psbn_nnzt_=4, psbn_zero_=5,psbn_ipc_=6
+  ! Duplicate coefficients handling
+  ! These are usually set while calling spcnv as one of its
+  ! optional arugments.
+  integer, parameter :: psbn_dupl_ovwrt_ = 0
+  integer, parameter :: psbn_dupl_add_   = 1
+  integer, parameter :: psbn_dupl_err_   = 2
+  integer, parameter :: psbn_dupl_def_   = psbn_dupl_ovwrt_
+  ! Matrix update mode
+  integer, parameter :: psbn_upd_srch_   = 98764
+  integer, parameter :: psbn_upd_perm_   = 98765
+  integer, parameter :: psbn_upd_dflt_   = psbn_upd_srch_
+  ! Mark a COO matrix with sorted entries.
+  integer, parameter :: psbn_isrtdcoo_   = 98761
+  integer, parameter :: psbn_maxjdrows_=8, psbn_minjdrows_=4
+  integer, parameter :: psbn_dbleint_=2
+  character(len=5)   :: psbn_fidef_='CSR'
+
+
+  type  :: psbn_base_sparse_mat
+    integer              :: m, n
+    integer, private     :: state 
+    logical, private     :: triangle, unitd, upper
+  contains 
+    procedure, pass(a) :: base_get_nrows
+    procedure, pass(a) :: base_get_ncols
+    procedure, pass(a) :: base_get_nzeros
+    procedure, pass(a) :: base_get_size
+    procedure, pass(a) :: base_get_state
+    procedure, pass(a) :: base_is_bld
+    procedure, pass(a) :: base_is_upd
+    procedure, pass(a) :: base_is_asb
+    procedure, pass(a) :: base_is_upper
+    procedure, pass(a) :: base_is_lower
+    procedure, pass(a) :: base_is_triangle
+    procedure, pass(a) :: base_is_unit
+    procedure, pass(a) :: base_get_neigh
+    procedure, pass(a) :: base_allocate_mn
+    procedure, pass(a) :: base_allocate_mnnz
+    procedure, pass(a) :: base_free
+    generic, public    :: allocate => base_allocate_mn, base_allocate_mnnz
+    generic, public    :: get_nrows => base_get_nrows
+    generic, public    :: get_ncols => base_get_ncols
+    generic, public    :: get_nzeros => base_get_nzeros
+    generic, public    :: get_size  => base_get_size
+    generic, public    :: get_state => base_get_state
+    generic, public    :: is_triangle => base_is_triangle
+    generic, public    :: is_unit => base_is_unit
+    generic, public    :: is_upper => base_is_upper
+    generic, public    :: is_lower => base_is_lower
+    generic, public    :: is_bld => base_is_bld
+    generic, public    :: is_upd => base_is_upd
+    generic, public    :: is_asb => base_is_asb
+    generic, public    :: get_neigh => base_get_neigh
+    generic, public    :: free => base_free
+    
+  end type psbn_base_sparse_mat
+  
+contains
+ 
+  function base_get_state(a) result(res)
+    class(psbn_base_sparse_mat), intent(in) :: a
+    integer :: res
+    res = a%state
+  end function base_get_state
+ 
+  function base_get_nrows(a) result(res)
+    class(psbn_base_sparse_mat), intent(in) :: a
+    integer :: res
+    res = a%m
+  end function base_get_nrows
+
+  function base_get_ncols(a) result(res)
+    class(psbn_base_sparse_mat), intent(in) :: a
+    integer :: res
+    res = a%n
+  end function base_get_ncols
+
+  function base_is_triangle(a) result(res)
+    class(psbn_base_sparse_mat), intent(in) :: a
+    logical :: res
+    res = a%triangle
+  end function base_is_triangle
+
+  function base_is_unit(a) result(res)
+    class(psbn_base_sparse_mat), intent(in) :: a
+    logical :: res
+    res = a%unitd
+  end function base_is_unit
+
+  function base_is_upper(a) result(res)
+    class(psbn_base_sparse_mat), intent(in) :: a
+    logical :: res
+    res = a%upper
+  end function base_is_upper
+
+  function base_is_lower(a) result(res)
+    class(psbn_base_sparse_mat), intent(in) :: a
+    logical :: res
+    res = .not.a%upper
+  end function base_is_lower
+
+  function base_is_bld(a) result(res)
+    class(psbn_base_sparse_mat), intent(in) :: a
+    logical :: res
+    res = (a%state == psbn_spmat_bld_)
+  end function base_is_bld
+
+  function base_is_upd(a) result(res)
+    class(psbn_base_sparse_mat), intent(in) :: a
+    logical :: res
+    res = (a%state == psbn_spmat_upd_)
+  end function base_is_upd
+
+  function base_is_asb(a) result(res)
+    class(psbn_base_sparse_mat), intent(in) :: a
+    logical :: res
+    res = (a%state == psbn_spmat_asb_)
+  end function base_is_asb
+
+
+  function base_get_nzeros(a) result(res)
+    use psb_error_mod
+    class(psbn_base_sparse_mat), intent(in) :: a
+    integer :: res
+    
+    Integer :: err_act
+    character(len=20)  :: name='base_get_nzeros'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    res = -1
+    ! This is the base version. If we get here
+    ! it means the derived class is incomplete,
+    ! so we throw an error.
+    call psb_errpush(700,name)
+          
+    if (err_act /= psb_act_ret_) then
+      call psb_error()
+    end if
+    return
+
+  end function base_get_nzeros
+
+  function base_get_size(a) result(res)
+    use psb_error_mod
+    class(psbn_base_sparse_mat), intent(in) :: a
+    integer :: res
+    
+    Integer :: err_act
+    character(len=20)  :: name='base_get_size'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    res = -1
+    ! This is the base version. If we get here
+    ! it means the derived class is incomplete,
+    ! so we throw an error.
+    call psb_errpush(700,name)
+          
+    if (err_act /= psb_act_ret_) then
+      call psb_error()
+    end if
+    return
+
+  end function base_get_size
+
+
+  subroutine base_get_neigh(a,idx,neigh,n,info,lev)
+    use psb_error_mod
+    class(psbn_base_sparse_mat), intent(in) :: a   
+    integer, intent(in)                :: idx 
+    integer, intent(out)               :: n   
+    integer, allocatable, intent(out)  :: neigh(:)
+    integer, intent(out)               :: info
+    integer, optional, intent(in)      :: lev 
+    
+    Integer :: err_act
+    character(len=20)  :: name='base_get_neigh'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    info = 700
+    ! This is the base version. If we get here
+    ! it means the derived class is incomplete,
+    ! so we throw an error.
+    call psb_errpush(700,name)
+
+    if (err_act /= psb_act_ret_) then
+      call psb_error()
+    end if
+    return
+
+  end subroutine base_get_neigh
+
+  subroutine  base_allocate_mn(m,n,a) 
+    use psb_error_mod
+    integer, intent(in) :: m,n
+    class(psbn_base_sparse_mat), intent(inout) :: a
+
+    Integer :: err_act
+    character(len=20)  :: name='base_allocate_mn'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    ! This is the base version. If we get here
+    ! it means the derived class is incomplete,
+    ! so we throw an error.
+    call psb_errpush(700,name)
+          
+    if (err_act /= psb_act_ret_) then
+      call psb_error()
+    end if
+    return
+
+  end subroutine base_allocate_mn
+
+  subroutine  base_allocate_mnnz(m,n,nz,a) 
+    use psb_error_mod
+    integer, intent(in) :: m,n,nz
+    class(psbn_base_sparse_mat), intent(inout) :: a
+    Integer :: err_act
+    character(len=20)  :: name='base_allocate_mnz'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    ! This is the base version. If we get here
+    ! it means the derived class is incomplete,
+    ! so we throw an error.
+    call psb_errpush(700,name)
+          
+    if (err_act /= psb_act_ret_) then
+      call psb_error()
+    end if
+    return
+
+  end subroutine base_allocate_mnnz
+
+  subroutine  base_free(a) 
+    use psb_error_mod
+    class(psbn_base_sparse_mat), intent(inout) :: a
+    Integer :: err_act
+    character(len=20)  :: name='base_free'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    ! This is the base version. If we get here
+    ! it means the derived class is incomplete,
+    ! so we throw an error.
+    call psb_errpush(700,name)
+          
+    if (err_act /= psb_act_ret_) then
+      call psb_error()
+    end if
+    return
+
+  end subroutine base_free
+
+end module psbn_base_mat_mod
+
+module psbn_d_base_mat_mod
+
+  use psbn_base_mat_mod
+  type, extends(psbn_base_sparse_mat) :: psbn_d_base_sparse_mat
+  contains
+    procedure, pass(a) :: d_base_csmv
+    procedure, pass(a) :: d_base_csmm
+    generic, public    :: psbn_csmm => d_base_csmm, d_base_csmv
+    procedure, pass(a) :: d_base_cssv
+    procedure, pass(a) :: d_base_cssm
+    generic, public    :: psbn_cssm => d_base_cssm, d_base_cssv
+    
+  end type psbn_d_base_sparse_mat
+
+contains 
+
+  subroutine d_base_csmm(alpha,a,x,beta,y,info,trans) 
+    use psb_error_mod
+    class(psbn_d_base_sparse_mat), intent(in) :: a
+    real(kind(1.d0)), intent(in)    :: alpha, beta, x(:,:)
+    real(kind(1.d0)), intent(inout) :: y(:,:)
+    integer, intent(out)            :: info
+    character, optional, intent(in) :: trans
+
+    Integer :: err_act
+    character(len=20)  :: name='d_base_csmm'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    ! This is the base version. If we get here
+    ! it means the derived class is incomplete,
+    ! so we throw an error.
+    call psb_errpush(700,name)
+          
+    if (err_act /= psb_act_ret_) then
+      call psb_error()
+    end if
+    return
+
+  end subroutine d_base_csmm
+
+  subroutine d_base_csmv(alpha,a,x,beta,y,info,trans) 
+    use psb_error_mod
+    class(psbn_d_base_sparse_mat), intent(in) :: a
+    real(kind(1.d0)), intent(in)    :: alpha, beta, x(:)
+    real(kind(1.d0)), intent(inout) :: y(:)
+    integer, intent(out)            :: info
+    character, optional, intent(in) :: trans
+
+    Integer :: err_act
+    character(len=20)  :: name='d_base_csmv'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    ! This is the base version. If we get here
+    ! it means the derived class is incomplete,
+    ! so we throw an error.
+    call psb_errpush(700,name)
+          
+    if (err_act /= psb_act_ret_) then
+      call psb_error()
+    end if
+    return
+
+
+  end subroutine d_base_csmv
+
+  subroutine d_base_cssm(alpha,a,x,beta,y,info,trans) 
+    use psb_error_mod
+    class(psbn_d_base_sparse_mat), intent(in) :: a
+    real(kind(1.d0)), intent(in)    :: alpha, beta, x(:,:)
+    real(kind(1.d0)), intent(inout) :: y(:,:)
+    integer, intent(out)            :: info
+    character, optional, intent(in) :: trans
+    
+    Integer :: err_act
+    character(len=20)  :: name='d_base_cssm'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    info = 700
+    ! This is the base version. If we get here
+    ! it means the derived class is incomplete,
+    ! so we throw an error.
+    call psb_errpush(700,name)
+          
+    if (err_act /= psb_act_ret_) then
+      call psb_error()
+    end if
+    return
+
+  end subroutine d_base_cssm
+
+  subroutine d_base_cssv(alpha,a,x,beta,y,info,trans) 
+    use psb_error_mod
+    class(psbn_d_base_sparse_mat), intent(in) :: a
+    real(kind(1.d0)), intent(in)    :: alpha, beta, x(:)
+    real(kind(1.d0)), intent(inout) :: y(:)
+    integer, intent(out)            :: info
+    character, optional, intent(in) :: trans
+
+    Integer :: err_act
+    character(len=20)  :: name='d_base_cssv'
+    logical, parameter :: debug=.false.
+
+    call psb_erractionsave(err_act)
+    info = 700
+    ! This is the base version. If we get here
+    ! it means the derived class is incomplete,
+    ! so we throw an error.
+    call psb_errpush(700,name)
+          
+    if (err_act /= psb_act_ret_) then
+      call psb_error()
+    end if
+    return
+
+
+  end subroutine d_base_cssv
+
+end module psbn_d_base_mat_mod
+
+module psbn_d_mat_mod
+
+  use psbn_d_base_mat_mod
+  
+  type :: psbn_d_sparse_mat
+
+    class(psbn_d_base_sparse_mat), allocatable  :: a 
+    
+  contains
+    
+    procedure, pass(a) :: d_csmv
+    procedure, pass(a) :: d_csmm
+    generic, public    :: psbn_csmm => d_csmm, d_csmv
+
+    procedure, pass(a) :: d_cssv
+    procedure, pass(a) :: d_cssm
+    generic, public    :: psbn_cssm => d_cssm, d_cssv
+    
+  end type psbn_d_sparse_mat
+
+contains 
+
+  subroutine d_csmm(alpha,a,x,beta,y,info,trans) 
+    class(psbn_d_sparse_mat), intent(in) :: a
+    real(kind(1.d0)), intent(in)    :: alpha, beta, x(:,:)
+    real(kind(1.d0)), intent(inout) :: y(:,:)
+    integer, intent(out)            :: info
+    character, optional, intent(in) :: trans
+
+    call a%a%psbn_csmm(alpha,x,beta,y,info,trans) 
+
+  end subroutine d_csmm
+
+  subroutine d_csmv(alpha,a,x,beta,y,info,trans) 
+    class(psbn_d_sparse_mat), intent(in) :: a
+    real(kind(1.d0)), intent(in)    :: alpha, beta, x(:)
+    real(kind(1.d0)), intent(inout) :: y(:)
+    integer, intent(out)            :: info
+    character, optional, intent(in) :: trans
+
+    call a%a%psbn_csmm(alpha,x,beta,y,info,trans) 
+
+  end subroutine d_csmv
+
+  subroutine d_cssm(alpha,a,x,beta,y,info,trans) 
+    class(psbn_d_sparse_mat), intent(in) :: a
+    real(kind(1.d0)), intent(in)    :: alpha, beta, x(:,:)
+    real(kind(1.d0)), intent(inout) :: y(:,:)
+    integer, intent(out)            :: info
+    character, optional, intent(in) :: trans
+    
+    call a%a%psbn_cssm(alpha,x,beta,y,info,trans) 
+
+  end subroutine d_cssm
+
+  subroutine d_cssv(alpha,a,x,beta,y,info,trans) 
+    class(psbn_d_sparse_mat), intent(in) :: a
+    real(kind(1.d0)), intent(in)    :: alpha, beta, x(:)
+    real(kind(1.d0)), intent(inout) :: y(:)
+    integer, intent(out)            :: info
+    character, optional, intent(in) :: trans
+
+    call a%a%psbn_cssm(alpha,x,beta,y,info,trans) 
+
+  end subroutine d_cssv
+
+end module psbn_d_mat_mod
+

base/newserial/psbn_coo_mat.f03

@@ -0,0 +1,289 @@
+module psbn_d_coo_sparse_mat_mod
+  
+  use psbn_d_base_mat_mod
+
+  type, extends(psbn_d_base_sparse_mat) :: psbn_d_coo_sparse_mat
+    integer              :: nnz,  state
+    logical              :: sorted
+    
+    integer, allocatable :: ia(:), ja(:)
+    real(kind(1.d0)), allocatable :: val(:)
+  contains
+    procedure, pass(a)  :: d_coo_get_nzeros
+    procedure, pass(a)  :: d_base_csmm => d_coo_csmm
+    procedure, pass(a)  :: d_base_csmv => d_coo_csmv
+    generic, public     :: base_get_nzeros => d_coo_get_nzeros
+
+
+  end type psbn_d_coo_sparse_mat
+
+contains 
+
+  function d_coo_get_nzeros(a) result(res)
+    class(psbn_d_coo_sparse_mat), intent(in) :: a
+    integer :: res
+    res = a%nnz
+  end function d_coo_get_nzeros
+  
+
+  subroutine d_coo_csmv(alpha,a,x,beta,y,info,trans) 
+    use psb_const_mod
+    class(psbn_d_coo_sparse_mat), intent(in) :: a
+    real(psb_dpk_), intent(in)          :: alpha, beta, x(:)
+    real(psb_dpk_), intent(inout)       :: y(:)
+    integer, intent(out)                :: info
+    character, optional, intent(in)     :: trans
+    
+    character :: trans_
+    integer   :: i,j,k,m,n, nnz, ir, jc
+    real(psb_dpk_) :: acc
+    logical   :: tra
+
+    if (present(trans)) then
+      trans_ = trans
+    else
+      trans_ = 'N'
+    end if
+    
+    tra = ((trans_=='T').or.(trans_=='t'))
+
+    if (tra) then 
+      m = a%get_ncols()
+      n = a%get_nrows()
+    else
+      n = a%get_ncols()
+      m = a%get_nrows()
+    end if
+    nnz = a%get_nzeros()
+
+    if (alpha == dzero) then
+      if (beta == dzero) then
+        do i = 1, m
+          y(i) = dzero
+        enddo
+      else
+        do  i = 1, m
+          y(i) = beta*y(i)
+        end do
+      endif
+      return
+    else 
+      if (.not.a%is_unit()) then 
+        if (beta == dzero) then
+          do i = 1, m
+            y(i) = dzero
+          enddo
+        else
+          do  i = 1, m
+            y(i) = beta*y(i)
+          end do
+        endif
+        
+      else  if (a%is_unit()) then 
+        if (beta == dzero) then
+          do i = 1, min(m,n)
+            y(i) = alpha*x(i)
+          enddo
+          do i = min(m,n)+1, m
+            y(i) = dzero
+          enddo
+        else
+          do  i = 1, min(m,n) 
+            y(i) = beta*y(i) + alpha*x(i)
+          end do
+          do i = min(m,n)+1, m
+            y(i) = beta*y(i)
+          enddo
+        endif
+        
+      endif
+
+    end if
+
+    if (.not.tra) then 
+      i    = 1
+      j    = i
+      if (nnz > 0) then 
+        ir   = a%ia(1) 
+        acc  = zero
+        do 
+          if (i>nnz) then 
+            y(ir) = y(ir) + alpha * acc
+            exit
+          endif
+          if (ia(i) /= ir) then 
+            y(ir) = y(ir) + alpha * acc
+            ir    = ia(i) 
+            acc   = zero
+          endif
+          acc     = acc + a%val(i) * x(a%ja(i))
+          i       = i + 1               
+        enddo
+      end if
+
+    else if (tra) then 
+      if (alpha.eq.done) then
+        i    = 1
+        do i=1,nnz
+          ir = a%ja(i)
+          jc = a%ia(i)
+          y(ir) = y(ir) +  a%val(i)*x(jc)
+        enddo
+        
+      else if (alpha.eq.-done) then
+        
+        do i=1,nnz
+          ir = a%ja(i)
+          jc = a%ia(i)
+          y(ir) = y(ir) - a%val(i)*x(jc)
+        enddo
+        
+      else                    
+        
+        do i=1,nnz
+          ir = ja(i)
+          jc = ia(i)
+          y(ir) = y(ir) + alpha*a%val(i)*x(jc)
+        enddo
+        
+      end if                  !.....end testing on alpha
+      
+    endif
+
+  end subroutine d_coo_csmv
+
+  subroutine d_coo_csmm(alpha,a,x,beta,y,info,trans) 
+    use psb_const_mod
+    class(psbn_d_coo_sparse_mat), intent(in) :: a
+    real(psb_dpk_), intent(in)          :: alpha, beta, x(:,:)
+    real(psb_dpk_), intent(inout)       :: y(:,:)
+    integer, intent(out)                :: info
+    character, optional, intent(in)     :: trans
+    
+    character :: trans_
+    integer   :: i,j,k,m,n, nnz, ir, jc, nc
+    real(psb_dpk_), allocatable  :: acc(:)
+    logical   :: tra
+
+    if (present(trans)) then
+      trans_ = trans
+    else
+      trans_ = 'N'
+    end if
+    
+    tra = ((trans_=='T').or.(trans_=='t'))
+
+    if (tra) then 
+      m = a%get_ncols()
+      n = a%get_nrows()
+    else
+      n = a%get_ncols()
+      m = a%get_nrows()
+    end if
+    nnz = a%get_nzeros()
+
+    nc = size(x,2) 
+    if (nc /= size(y,2)) then 
+      write(0,*) 'Mismatch in column sizes!!'
+      return
+    end if
+    allocate(acc(nc))
+
+    if (alpha == dzero) then
+      if (beta == dzero) then
+        do i = 1, m
+          y(i,:) = dzero
+        enddo
+      else
+        do  i = 1, m
+          y(i,:) = beta*y(i,:)
+        end do
+      endif
+      return
+    else 
+      if (.not.a%is_unit()) then 
+        if (beta == dzero) then
+          do i = 1, m
+            y(i,:) = dzero
+          enddo
+        else
+          do  i = 1, m
+            y(i,:) = beta*y(i,:)
+          end do
+        endif
+        
+      else  if (a%is_unit()) then 
+        if (beta == dzero) then
+          do i = 1, min(m,n)
+            y(i,:) = alpha*x(i,:)
+          enddo
+          do i = min(m,n)+1, m
+            y(i,:) = dzero
+          enddo
+        else
+          do  i = 1, min(m,n) 
+            y(i,:) = beta*y(i,:) + alpha*x(i,:)
+          end do
+          do i = min(m,n)+1, m
+            y(i,:) = beta*y(i,:)
+          enddo
+        endif
+        
+      endif
+
+    end if
+
+    if (.not.tra) then 
+      i    = 1
+      j    = i
+      if (nnz > 0) then 
+        ir   = a%ia(1) 
+        acc  = zero
+        do 
+          if (i>nnz) then 
+            y(ir,:) = y(ir,:) + alpha * acc
+            exit
+          endif
+          if (ia(i) /= ir) then 
+            y(ir,:) = y(ir,:) + alpha * acc
+            ir    = ia(i) 
+            acc   = zero
+          endif
+          acc     = acc + a%val(i) * x(a%ja(i),:)
+          i       = i + 1               
+        enddo
+      end if
+
+    else if (tra) then 
+      if (alpha.eq.done) then
+        i    = 1
+        do i=1,nnz
+          ir = a%ja(i)
+          jc = a%ia(i)
+          y(ir,:) = y(ir,:) +  a%val(i)*x(jc,:)
+        enddo
+        
+      else if (alpha.eq.-done) then
+        
+        do i=1,nnz
+          ir = a%ja(i)
+          jc = a%ia(i)
+          y(ir,:) = y(ir,:) - a%val(i)*x(jc,:)
+        enddo
+        
+      else                    
+        
+        do i=1,nnz
+          ir = ja(i)
+          jc = ia(i)
+          y(ir,:) = y(ir,:) + alpha*a%val(i)*x(jc,:)
+        enddo
+        
+      end if                  !.....end testing on alpha
+      
+    endif
+
+  end subroutine d_coo_csmm
+  
+end module psbn_d_coo_sparse_mat_mod
+

base/newserial/psbn_csr_mat.f03

@@ -0,0 +1,874 @@
+module psbn_d_csr_sparse_mat_mod
+
+  use psbn_d_base_mat_mod
+
+  type, extends(psbn_d_base_sparse_mat) :: psbn_d_csr_sparse_mat
+    logical              :: sorted
+    integer, allocatable :: irp(:), ja(:)
+    real(kind(1.d0)), allocatable :: val(:)
+  contains
+    procedure, pass(a)  :: d_csr_get_nzeros
+    procedure, pass(a)  :: d_base_csmm => d_csr_csmm
+    procedure, pass(a)  :: d_base_csmv => d_csr_csmv
+    generic, public     :: base_get_nzeros => d_csr_get_nzeros
+    procedure, pass(a)  :: d_base_cssm => d_csr_cssm
+    procedure, pass(a)  :: d_base_cssv => d_csr_cssv
+
+
+  end type psbn_d_csr_sparse_mat
+
+contains 
+
+  function d_csr_get_nzeros(a) result(res)
+    class(psbn_d_csr_sparse_mat), intent(in) :: a
+    integer :: res
+    res = a%irp(a%m+1)-1
+  end function d_csr_get_nzeros
+  
+
+  subroutine d_csr_csmv(alpha,a,x,beta,y,info,trans) 
+    use psb_const_mod
+    class(psbn_d_csr_sparse_mat), intent(in) :: a
+    real(psb_dpk_), intent(in)          :: alpha, beta, x(:)
+    real(psb_dpk_), intent(inout)       :: y(:)
+    integer, intent(out)                :: info
+    character, optional, intent(in)     :: trans
+    
+    character :: trans_
+    integer   :: i,j,k,m,n, nnz, ir, jc
+    real(psb_dpk_) :: acc
+    logical   :: tra
+
+    if (present(trans)) then
+      trans_ = trans
+    else
+      trans_ = 'N'
+    end if
+    
+    if (.not.a%is_asb()) then 
+      write(0,*) 'Error: csmv called on an unassembled mat'
+    end if
+    
+    tra = ((trans_=='T').or.(trans_=='t'))
+
+    if (tra) then 
+      m = a%get_ncols()
+      n = a%get_nrows()
+    else
+      n = a%get_ncols()
+      m = a%get_nrows()
+    end if
+    
+    
+    if (alpha == dzero) then
+      if (beta == dzero) then
+        do i = 1, m
+          y(i) = dzero
+        enddo
+      else
+        do  i = 1, m
+          y(i) = beta*y(i)
+        end do
+      endif
+      return
+    end if
+
+    if (.not.tra) then 
+
+      if (beta == dzero) then 
+
+        if (alpha == done) then 
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j))          
+            enddo
+            y(i) = acc
+          end do
+        
+        else if (alpha == -done) then 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j))          
+            enddo
+            y(i) = -acc
+          end do
+      
+        else 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j))          
+            enddo
+            y(i) = alpha*acc
+          end do
+
+        end if
+        
+        
+      else if (beta == done) then 
+
+        if (alpha == done) then 
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j))          
+            enddo
+            y(i) = y(i) + acc
+          end do
+        
+        else if (alpha == -done) then 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j))          
+            enddo
+            y(i) = y(i) -acc
+          end do
+      
+        else 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j))          
+            enddo
+            y(i) = y(i) + alpha*acc
+          end do
+
+        end if
+        
+      else if (beta == -done) then 
+
+        if (alpha == done) then 
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j))          
+            enddo
+            y(i) = -y(i) + acc
+          end do
+        
+        else if (alpha == -done) then 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j))          
+            enddo
+            y(i) = -y(i) -acc
+          end do
+      
+        else 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j))          
+            enddo
+            y(i) = -y(i) + alpha*acc
+          end do
+
+        end if        
+
+      else 
+
+        if (alpha == done) then 
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j))          
+            enddo
+            y(i) = beta*y(i) + acc
+          end do
+        
+        else if (alpha == -done) then 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j))          
+            enddo
+            y(i) = beta*y(i) - acc
+          end do
+      
+        else 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j))          
+            enddo
+            y(i) = beta*y(i) + alpha*acc
+          end do
+
+        end if
+                
+      end if
+
+    else if (tra) then 
+
+      if (beta == dzero) then 
+        do i=1, m
+          y(i) = dzero
+        end do
+      else if (beta == done) then 
+        ! Do nothing
+      else if (beta == -done) then 
+        do i=1, m
+          y(i) = -y(i) 
+        end do
+      else
+        do i=1, m
+          y(i) = beta*y(i) 
+        end do
+      end if
+
+      if (alpha.eq.done) then
+
+        do i=1,n
+          do j=a%irp(i), a%irp(i+1)-1
+            ir = a%ja(j)
+            y(ir) = y(ir) +  a%val(j)*x(i)
+          end do
+        enddo
+        
+      else if (alpha.eq.-done) then
+        
+        do i=1,n
+          do j=a%irp(i), a%irp(i+1)-1
+            ir = a%ja(j)
+            y(ir) = y(ir) -  a%val(j)*x(i)
+          end do
+        enddo
+        
+      else                    
+        
+        do i=1,n
+          do j=a%irp(i), a%irp(i+1)-1
+            ir = a%ja(j)
+            y(ir) = y(ir) + alpha*a%val(j)*x(i)
+          end do
+        enddo
+        
+      end if               
+      
+    endif
+      
+    if (a%is_unit()) then 
+      do i=1, min(m,n)
+        y(i) = y(i) + alpha*x(i)
+      end do
+    end if
+
+
+  end subroutine d_csr_csmv
+
+  subroutine d_csr_csmm(alpha,a,x,beta,y,info,trans) 
+    use psb_const_mod
+    class(psbn_d_csr_sparse_mat), intent(in) :: a
+    real(psb_dpk_), intent(in)          :: alpha, beta, x(:,:)
+    real(psb_dpk_), intent(inout)       :: y(:,:)
+    integer, intent(out)                :: info
+    character, optional, intent(in)     :: trans
+    
+    character :: trans_
+    integer   :: i,j,k,m,n, nnz, ir, jc, nc
+    real(psb_dpk_), allocatable  :: acc(:)
+    logical   :: tra
+
+    if (present(trans)) then
+      trans_ = trans
+    else
+      trans_ = 'N'
+    end if
+    
+    tra = ((trans_=='T').or.(trans_=='t'))
+    if (.not.a%is_asb()) then 
+      write(0,*) 'Error: csmv called on an unassembled mat'
+    end if
+    
+    if (tra) then 
+      m = a%get_ncols()
+      n = a%get_nrows()
+    else
+      n = a%get_ncols()
+      m = a%get_nrows()
+    end if
+
+    nc = size(x,2) 
+    if (nc /= size(y,2)) then 
+      write(0,*) 'Mismatch in column sizes!!'
+      return
+    end if
+    allocate(acc(nc))
+
+   if (alpha == dzero) then
+      if (beta == dzero) then
+        do i = 1, m
+          y(i,:) = dzero
+        enddo
+      else
+        do  i = 1, m
+          y(i,:) = beta*y(i,:)
+        end do
+      endif
+      return
+    end if
+
+    if (.not.tra) then 
+
+      if (beta == dzero) then 
+
+        if (alpha == done) then 
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j),:)          
+            enddo
+            y(i,:) = acc
+          end do
+        
+        else if (alpha == -done) then 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j),:)          
+            enddo
+            y(i,:) = -acc
+          end do
+      
+        else 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j),:)          
+            enddo
+            y(i,:) = alpha*acc
+          end do
+
+        end if
+        
+        
+      else if (beta == done) then 
+
+        if (alpha == done) then 
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j),:)          
+            enddo
+            y(i,:) = y(i,:) + acc
+          end do
+        
+        else if (alpha == -done) then 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j),:)          
+            enddo
+            y(i,:) = y(i,:) -acc
+          end do
+      
+        else 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j),:)          
+            enddo
+            y(i,:) = y(i,:) + alpha*acc
+          end do
+
+        end if
+        
+      else if (beta == -done) then 
+
+        if (alpha == done) then 
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j),:)          
+            enddo
+            y(i,:) = -y(i,:) + acc
+          end do
+        
+        else if (alpha == -done) then 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j),:)          
+            enddo
+            y(i,:) = -y(i,:) -acc
+          end do
+      
+        else 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j),:)          
+            enddo
+            y(i,:) = -y(i,:) + alpha*acc
+          end do
+
+        end if        
+
+      else 
+
+        if (alpha == done) then 
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j),:)          
+            enddo
+            y(i,:) = beta*y(i,:) + acc
+          end do
+        
+        else if (alpha == -done) then 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j),:)          
+            enddo
+            y(i,:) = beta*y(i,:) - acc
+          end do
+      
+        else 
+
+          do i=1,m 
+            acc  = zero
+            do j=a%irp(i), a%irp(i+1)-1
+              acc  = acc + a%val(j) * x(a%ja(j),:)          
+            enddo
+            y(i,:) = beta*y(i,:) + alpha*acc
+          end do
+
+        end if
+                
+      end if
+
+    else if (tra) then 
+
+      if (beta == dzero) then 
+        do i=1, m
+          y(i,:) = dzero
+        end do
+      else if (beta == done) then 
+        ! Do nothing
+      else if (beta == -done) then 
+        do i=1, m
+          y(i,:) = -y(i,:) 
+        end do
+      else
+        do i=1, m
+          y(i,:) = beta*y(i,:) 
+        end do
+      end if
+
+      if (alpha.eq.done) then
+
+        do i=1,n
+          do j=a%irp(i), a%irp(i+1)-1
+            ir = a%ja(j)
+            y(ir,:) = y(ir,:) +  a%val(j)*x(i,:)
+          end do
+        enddo
+        
+      else if (alpha.eq.-done) then
+        
+        do i=1,n
+          do j=a%irp(i), a%irp(i+1)-1
+            ir = a%ja(j)
+            y(ir,:) = y(ir,:) -  a%val(j)*x(i,:)
+          end do
+        enddo
+        
+      else                    
+        
+        do i=1,n
+          do j=a%irp(i), a%irp(i+1)-1
+            ir = a%ja(j)
+            y(ir,:) = y(ir,:) + alpha*a%val(j)*x(i,:)
+          end do
+        enddo
+        
+      end if               
+      
+    endif
+      
+    if (a%is_unit()) then 
+      do i=1, min(m,n)
+        y(i,:) = y(i,:) + alpha*x(i,:)
+      end do
+    end if
+    
+
+  end subroutine d_csr_csmm
+
+
+  subroutine d_csr_cssv(alpha,a,x,beta,y,info,trans) 
+    use psb_const_mod
+    class(psbn_d_csr_sparse_mat), intent(in) :: a
+    real(psb_dpk_), intent(in)          :: alpha, beta, x(:)
+    real(psb_dpk_), intent(inout)       :: y(:)
+    integer, intent(out)                :: info
+    character, optional, intent(in)     :: trans
+    
+    character :: trans_
+    integer   :: i,j,k,m,n, nnz, ir, jc
+    real(psb_dpk_) :: acc
+    real(psb_dpk_), allocatable :: tmp(:)
+    logical   :: tra
+
+    if (present(trans)) then
+      trans_ = trans
+    else
+      trans_ = 'N'
+    end if
+    if (.not.a%is_asb()) then 
+      write(0,*) 'Error: cssv called on an unassembled mat'
+    end if
+        
+    tra = ((trans_=='T').or.(trans_=='t'))
+    m = a%get_nrows()
+    
+    if (.not. (a%is_triangle())) then 
+      write(0,*) 'Called SV on a non-triangular mat!'
+    end if
+
+
+    if (alpha == dzero) then
+      if (beta == dzero) then
+        do i = 1, m
+          y(i) = dzero
+        enddo
+      else
+        do  i = 1, m
+          y(i) = beta*y(i)
+        end do
+      endif
+      return
+    end if
+
+    if (beta == dzero) then 
+      call inner_csrsv(tra,a,x,y)
+      do  i = 1, m
+        y(i) = alpha*y(i)
+      end do
+    else 
+      allocate(tmp(m), stat=info) 
+      if (info /= 0) then 
+        write(0,*) 'Memory allocation error in CSRSV '
+        return
+      end if
+      tmp(1:m) = x(1:m)
+      call inner_csrsv(tra,a,tmp,y)
+      do  i = 1, m
+        y(i) = alpha*tmp(i) + beta*y(i)
+      end do
+    end if
+    
+
+  contains 
+
+    subroutine inner_csrsv(tra,a,x,y) 
+      use psb_const_mod
+      logical, intent(in)                 :: tra  
+      class(psbn_d_csr_sparse_mat), intent(in) :: a
+      real(psb_dpk_), intent(in)          :: x(:)
+      real(psb_dpk_), intent(out)         :: y(:)
+
+      integer :: i,j,k,m, ir, jc
+      real(psb_dpk_) :: acc
+      
+      if (.not.tra) then 
+
+        if (a%is_lower()) then 
+          if (a%is_unit()) then 
+            do i=1, a%get_nrows()
+              acc = dzero 
+              do j=a%irp(i), a%irp(i+1)-1
+                acc = acc + a%val(j)*x(a%ja(j))
+              end do
+              y(i) = x(i) - acc
+            end do
+          else if (.not.a%is_unit()) then 
+            do i=1, a%get_nrows()
+              acc = dzero 
+              do j=a%irp(i), a%irp(i+1)-2
+                acc = acc + a%val(j)*x(a%ja(j))
+              end do
+              y(i) = (x(i) - acc)/a%val(a%irp(i+1)-1)
+            end do
+          end if
+        else if (a%is_upper()) then 
+
+          if (a%is_unit()) then 
+            do i=a%get_nrows(), 1, -1 
+              acc = dzero 
+              do j=a%irp(i), a%irp(i+1)-1
+                acc = acc + a%val(j)*x(a%ja(j))
+              end do
+              y(i) = x(i) - acc
+            end do
+          else if (.not.a%is_unit()) then 
+            do i=a%get_nrows(), 1, -1 
+              acc = dzero 
+              do j=a%irp(i)+1, a%irp(i+1)-1
+                acc = acc + a%val(j)*x(a%ja(j))
+              end do
+              y(i) = (x(i) - acc)/a%val(a%irp(i))
+            end do
+          end if
+          
+        end if
+
+      else if (tra) then 
+
+        do i=1, a%get_nrows()
+          y(i) = x(i)
+        end do
+
+        if (a%is_lower()) then 
+          if (a%is_unit()) then 
+            do i=a%get_nrows(), 1, -1
+              y(i) = y(i)/a%val(a%irp(i+1)-1)
+              acc = y(i) 
+              do j=a%irp(i), a%irp(i+1)-1
+                jc    = a%ja(j)
+                y(jc) = y(jc) - a%val(j)*acc 
+              end do
+            end do
+          else if (.not.a%is_unit()) then 
+            do i=a%get_nrows(), 1, -1
+              y(i) = y(i)/a%val(a%irp(i+1)-1)
+              acc  = y(i) 
+              do j=a%irp(i), a%irp(i+1)-2
+                jc    = a%ja(j)
+                y(jc) = y(jc) - a%val(j)*acc 
+              end do
+            end do
+          end if
+        else if (a%is_upper()) then 
+
+          if (a%is_unit()) then 
+            do i=1, a%get_nrows()
+              acc  = y(i) 
+              do j=a%irp(i), a%irp(i+1)-1
+                jc    = a%ja(j)
+                y(jc) = y(jc) - a%val(j)*acc 
+              end do
+            end do
+          else if (.not.a%is_unit()) then 
+            do i=1, a%get_nrows()
+              y(i) = y(i)/a%val(a%irp(i))
+              acc  = y(i) 
+              do j=a%irp(i)+1, a%irp(i+1)-1
+                jc    = a%ja(j)
+                y(jc) = y(jc) - a%val(j)*acc 
+              end do
+            end do
+          end if
+          
+        end if
+      end if
+    end subroutine inner_csrsv
+
+  end subroutine d_csr_cssv
+
+
+
+  subroutine d_csr_cssm(alpha,a,x,beta,y,info,trans) 
+    use psb_const_mod
+    class(psbn_d_csr_sparse_mat), intent(in) :: a
+    real(psb_dpk_), intent(in)          :: alpha, beta, x(:,:)
+    real(psb_dpk_), intent(inout)       :: y(:,:)
+    integer, intent(out)                :: info
+    character, optional, intent(in)     :: trans
+    
+    character :: trans_
+    integer   :: i,j,k,m,n, nnz, ir, jc, nc
+    real(psb_dpk_) :: acc
+    real(psb_dpk_), allocatable :: tmp(:,:)
+    logical   :: tra
+
+    if (present(trans)) then
+      trans_ = trans
+    else
+      trans_ = 'N'
+    end if
+    if (.not.a%is_asb()) then 
+      write(0,*) 'Error: cssm called on an unassembled mat'
+    end if
+        
+    tra = ((trans_=='T').or.(trans_=='t'))
+    m   = a%get_nrows()
+    nc  = size(x,2) 
+    if (nc /= size(y,2)) then 
+      write(0,*) 'Mismatch in column sizes!!'
+      return
+    end if
+    
+    if (.not. (a%is_triangle())) then 
+      write(0,*) 'Called SM on a non-triangular mat!'
+    end if
+
+
+    if (alpha == dzero) then
+      if (beta == dzero) then
+        do i = 1, m
+          y(i,:) = dzero
+        enddo
+      else
+        do  i = 1, m
+          y(i,:) = beta*y(i,:)
+        end do
+      endif
+      return
+    end if
+
+    if (beta == dzero) then 
+      call inner_csrsm(tra,a,x,y,info)
+      do  i = 1, m
+        y(i,:) = alpha*y(i,:)
+      end do
+    else 
+      allocate(tmp(m,nc), stat=info) 
+      if (info /= 0) then 
+        write(0,*) 'Memory allocation error in CSRSM '
+        return
+      end if
+      tmp(1:m,:) = x(1:m,:)
+      call inner_csrsm(tra,a,tmp,y,info)
+      do  i = 1, m
+        y(i,:) = alpha*tmp(i,:) + beta*y(i,:)
+      end do
+    end if
+    
+
+  contains 
+
+    subroutine inner_csrsm(tra,a,x,y,info) 
+      use psb_const_mod
+      logical, intent(in)                 :: tra  
+      class(psbn_d_csr_sparse_mat), intent(in) :: a
+      real(psb_dpk_), intent(in)          :: x(:,:)
+      real(psb_dpk_), intent(out)         :: y(:,:)
+      integer, intent(out)                :: info
+      integer :: i,j,k,m, ir, jc
+      real(psb_dpk_), allocatable  :: acc(:)
+      
+      allocate(acc(size(x,2)), stat=info)
+      if (info /= 0) then 
+        write(0,*) 'Memory allocation error in CSRSM '
+        return
+      end if
+      
+      
+      if (.not.tra) then 
+
+        if (a%is_lower()) then 
+          if (a%is_unit()) then 
+            do i=1, a%get_nrows()
+              acc = dzero 
+              do j=a%irp(i), a%irp(i+1)-1
+                acc = acc + a%val(j)*x(a%ja(j),:)
+              end do
+              y(i,:) = x(i,:) - acc
+            end do
+          else if (.not.a%is_unit()) then 
+            do i=1, a%get_nrows()
+              acc = dzero 
+              do j=a%irp(i), a%irp(i+1)-2
+                acc = acc + a%val(j)*x(a%ja(j),:)
+              end do
+              y(i,:) = (x(i,:) - acc)/a%val(a%irp(i+1)-1)
+            end do
+          end if
+        else if (a%is_upper()) then 
+
+          if (a%is_unit()) then 
+            do i=a%get_nrows(), 1, -1 
+              acc = dzero 
+              do j=a%irp(i), a%irp(i+1)-1
+                acc = acc + a%val(j)*x(a%ja(j),:)
+              end do
+              y(i,:) = x(i,:) - acc
+            end do
+          else if (.not.a%is_unit()) then 
+            do i=a%get_nrows(), 1, -1 
+              acc = dzero 
+              do j=a%irp(i)+1, a%irp(i+1)-1
+                acc = acc + a%val(j)*x(a%ja(j),:)
+              end do
+              y(i,:) = (x(i,:) - acc)/a%val(a%irp(i))
+            end do
+          end if
+          
+        end if
+
+      else if (tra) then 
+
+        do i=1, a%get_nrows()
+          y(i,:) = x(i,:)
+        end do
+
+        if (a%is_lower()) then 
+          if (a%is_unit()) then 
+            do i=a%get_nrows(), 1, -1
+              y(i,:) = y(i,:)/a%val(a%irp(i+1)-1)
+              acc = y(i,:) 
+              do j=a%irp(i), a%irp(i+1)-1
+                jc    = a%ja(j)
+                y(jc,:) = y(jc,:) - a%val(j)*acc 
+              end do
+            end do
+          else if (.not.a%is_unit()) then 
+            do i=a%get_nrows(), 1, -1
+              y(i,:) = y(i,:)/a%val(a%irp(i+1)-1)
+              acc  = y(i,:) 
+              do j=a%irp(i), a%irp(i+1)-2
+                jc    = a%ja(j)
+                y(jc,:) = y(jc,:) - a%val(j)*acc 
+              end do
+            end do
+          end if
+        else if (a%is_upper()) then 
+
+          if (a%is_unit()) then 
+            do i=1, a%get_nrows()
+              acc  = y(i,:) 
+              do j=a%irp(i), a%irp(i+1)-1
+                jc    = a%ja(j)
+                y(jc,:) = y(jc,:) - a%val(j)*acc 
+              end do
+            end do
+          else if (.not.a%is_unit()) then 
+            do i=1, a%get_nrows()
+              y(i,:) = y(i,:)/a%val(a%irp(i))
+              acc    = y(i,:) 
+              do j=a%irp(i)+1, a%irp(i+1)-1
+                jc      = a%ja(j)
+                y(jc,:) = y(jc,:) - a%val(j)*acc 
+              end do
+            end do
+          end if
+          
+        end if
+      end if
+    end subroutine inner_csrsm
+
+  end subroutine d_csr_cssm
+
+end module psbn_d_csr_sparse_mat_mod
+

krylov/psb_zrgmres.f90

@@ -348,7 +348,7 @@ Subroutine psb_zrgmres(a,prec,b,x,eps,desc_a,info,itmax,iter,err,itrace,irst,ist
       rti  = h(i,i)
       rti1 = h(i+1,i) 
       call zrotg(rti,rti1,tmp,s(i))
-      c(i) = cmplx(tmp,dzero)
+      c(i) = cmplx(tmp,dzero,kind=psb_dpk_)
       call zrot(1,h(i,i),1,h(i+1,i),1,real(c(i)),s(i))
       h(i+1,i) = zzero
       call zrot(1,rs(i),1,rs(i+1),1,real(c(i)),s(i))
@@ -574,11 +574,11 @@ contains
     !
 
     scale = cdabs(ca) + cdabs(cb)
-    norm = scale*dsqrt((cdabs(ca/dcmplx(scale,0.0d0)))**2 +&
+    norm = scale*dsqrt((cdabs(ca/cmplx(scale,0.0d0,kind=psb_dpk_)))**2 +&
-         &   (cdabs(cb/dcmplx(scale,0.0d0)))**2)
+         &   (cdabs(cb/cmplx(scale,0.0d0,kind=psb_dpk_)))**2)
     alpha = ca /cdabs(ca)
     c = cdabs(ca) / norm
-    s = alpha * dconjg(cb) / norm
+    s = alpha * conjg(cb) / norm
     ca = alpha * norm
     !
 

util/psb_mmio_mod.f90

@@ -240,7 +240,7 @@ contains
       do j=1, ncol
         do i=1, nrow
           read(infile,fmt=*,end=902) bre,bim
-          b(i,j) = cmplx(bre,bim)
+          b(i,j) = cmplx(bre,bim,kind=psb_spk_)
         end do
       end do
       
@@ -317,7 +317,7 @@ contains
       do j=1, ncol
         do i=1, nrow
           read(infile,fmt=*,end=902) bre,bim
-          b(i,j) = cmplx(bre,bim)
+          b(i,j) = cmplx(bre,bim,kind=psb_dpk_)
         end do
       end do
       
@@ -689,7 +689,7 @@ contains
       a%descra = 'G'      
       do i=1,nnzero
         read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),are,aim
-        a%aspk(i) = cmplx(are,aim)
+        a%aspk(i) = cmplx(are,aim,kind=psb_spk_)
       end do
       a%infoa(psb_nnz_) = nnzero
       
@@ -704,7 +704,7 @@ contains
       a%descra = 'G'      
       do i=1,nnzero
         read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),are,aim
-        a%aspk(i) = cmplx(are,aim)
+        a%aspk(i) = cmplx(are,aim,kind=psb_spk_)
       end do
 
       nzr = nnzero
@@ -728,7 +728,7 @@ contains
       a%descra = 'G'      
       do i=1,nnzero
         read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),are,aim
-        a%aspk(i) = cmplx(are,aim)
+        a%aspk(i) = cmplx(are,aim,kind=psb_spk_)
       end do
 
       nzr = nnzero
@@ -864,7 +864,7 @@ contains
       a%descra = 'G'      
       do i=1,nnzero
         read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),are,aim
-        a%aspk(i) = cmplx(are,aim)
+        a%aspk(i) = cmplx(are,aim,kind=psb_dpk_)
       end do
       a%infoa(psb_nnz_) = nnzero
       
@@ -879,7 +879,7 @@ contains
       a%descra = 'G'      
       do i=1,nnzero
         read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),are,aim
-        a%aspk(i) = cmplx(are,aim)
+        a%aspk(i) = cmplx(are,aim,kind=psb_dpk_)
       end do
 
       nzr = nnzero
@@ -903,7 +903,7 @@ contains
       a%descra = 'G'      
       do i=1,nnzero
         read(infile,fmt=*,end=902) a%ia1(i),a%ia2(i),are,aim
-        a%aspk(i) = cmplx(are,aim)
+        a%aspk(i) = cmplx(are,aim,kind=psb_dpk_)
       end do
 
       nzr = nnzero