From 2d7fef55ed41de8bb949e1fbc33681c836c8b1b7 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Tue, 26 Mar 2013 15:45:30 +0000 Subject: [PATCH] psblas-3.0-maint: test/serial/d_matgen.F90 test/serial/psb_d_cxx_impl.f90 test/serial/psb_d_cxx_mat_mod.f90 test/serial/psb_d_cyy_impl.f90 test/serial/psb_d_cyy_mat_mod.f90 Fixed example to compile properly. --- test/serial/d_matgen.F90 | 36 +- test/serial/psb_d_cxx_impl.f90 | 733 ++++++-- test/serial/psb_d_cxx_mat_mod.f90 | 214 ++- test/serial/psb_d_cyy_impl.f90 | 2712 ++++++++++++++++++++++++++++- test/serial/psb_d_cyy_mat_mod.f90 | 271 ++- 5 files changed, 3643 insertions(+), 323 deletions(-) diff --git a/test/serial/d_matgen.F90 b/test/serial/d_matgen.F90 index 18a20c0df..ff2658d85 100644 --- a/test/serial/d_matgen.F90 +++ b/test/serial/d_matgen.F90 @@ -1,3 +1,38 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.0 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ +! +! This is a simple test program to show how to plug a new matrix format +! into the library framework. +! ! program d_matgen use psb_base_mod @@ -17,7 +52,6 @@ program d_matgen ! sparse matrix and preconditioner type(psb_dspmat_type) :: a -!!$ type(psb_dprec_type) :: prec ! descriptor type(psb_desc_type) :: desc_a ! dense matrices diff --git a/test/serial/psb_d_cxx_impl.f90 b/test/serial/psb_d_cxx_impl.f90 index 50d7b3c0e..0dad379aa 100644 --- a/test/serial/psb_d_cxx_impl.f90 +++ b/test/serial/psb_d_cxx_impl.f90 @@ -1,3 +1,34 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.0 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ ! == =================================== ! @@ -13,7 +44,8 @@ ! == =================================== subroutine psb_d_cxx_csmv(alpha,a,x,beta,y,info,trans) - use psb_base_mod + use psb_error_mod + use psb_string_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_csmv implicit none class(psb_d_cxx_sparse_mat), intent(in) :: a @@ -25,8 +57,9 @@ subroutine psb_d_cxx_csmv(alpha,a,x,beta,y,info,trans) character :: trans_ integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc real(psb_dpk_) :: acc - logical :: tra + logical :: tra, ctra integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cxx_csmv' logical, parameter :: debug=.false. @@ -46,9 +79,10 @@ subroutine psb_d_cxx_csmv(alpha,a,x,beta,y,info,trans) endif - tra = (psb_toupper(trans_) == 'T').or.(psb_toupper(trans_)=='C') + tra = (psb_toupper(trans_) == 'T') + ctra = (psb_toupper(trans_) == 'C') - if (tra) then + if (tra.or.ctra) then m = a%get_ncols() n = a%get_nrows() else @@ -57,21 +91,23 @@ subroutine psb_d_cxx_csmv(alpha,a,x,beta,y,info,trans) end if if (size(x,1) psb_d_cxx_csmm implicit none class(psb_d_cxx_sparse_mat), intent(in) :: a @@ -322,8 +406,9 @@ subroutine psb_d_cxx_csmm(alpha,a,x,beta,y,info,trans) character :: trans_ integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc, nc real(psb_dpk_), allocatable :: acc(:) - logical :: tra + logical :: tra, ctra integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cxx_csmm' logical, parameter :: debug=.false. @@ -341,9 +426,10 @@ subroutine psb_d_cxx_csmm(alpha,a,x,beta,y,info,trans) call psb_errpush(info,name) goto 9999 endif - tra = (psb_toupper(trans_) == 'T').or.(psb_toupper(trans_)=='C') + tra = (psb_toupper(trans_) == 'T') + ctra = (psb_toupper(trans_) == 'C') - if (tra) then + if (tra.or.ctra) then m = a%get_ncols() n = a%get_nrows() else @@ -352,14 +438,16 @@ subroutine psb_d_cxx_csmm(alpha,a,x,beta,y,info,trans) end if if (size(x,1) psb_d_cxx_cssv implicit none class(psb_d_cxx_sparse_mat), intent(in) :: a @@ -624,8 +759,9 @@ subroutine psb_d_cxx_cssv(alpha,a,x,beta,y,info,trans) integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: tmp(:) - logical :: tra + logical :: tra,ctra integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cxx_cssv' logical, parameter :: debug=.false. @@ -642,7 +778,8 @@ subroutine psb_d_cxx_cssv(alpha,a,x,beta,y,info,trans) goto 9999 endif - tra = (psb_toupper(trans_) == 'T').or.(psb_toupper(trans_)=='C') + tra = (psb_toupper(trans_) == 'T') + ctra = (psb_toupper(trans_) == 'C') m = a%get_nrows() if (.not. (a%is_triangle())) then @@ -652,14 +789,16 @@ subroutine psb_d_cxx_cssv(alpha,a,x,beta,y,info,trans) end if if (size(x) psb_d_cxx_cssm implicit none class(psb_d_cxx_sparse_mat), intent(in) :: a @@ -839,8 +1027,9 @@ subroutine psb_d_cxx_cssm(alpha,a,x,beta,y,info,trans) integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc, nc real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: tmp(:,:) - logical :: tra + logical :: tra, ctra integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cxx_cssm' logical, parameter :: debug=.false. @@ -859,7 +1048,8 @@ subroutine psb_d_cxx_cssm(alpha,a,x,beta,y,info,trans) endif - tra = (psb_toupper(trans_) == 'T').or.(psb_toupper(trans_)=='C') + tra = (psb_toupper(trans_) == 'T') + ctra = (psb_toupper(trans_) == 'C') m = a%get_nrows() nc = min(size(x,2) , size(y,2)) @@ -885,8 +1075,8 @@ subroutine psb_d_cxx_cssm(alpha,a,x,beta,y,info,trans) end if if (beta == dzero) then - call inner_cxxsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,& - & a%irp,a%ja,a%val,x,size(x,1),y,size(y,1),info) + call inner_cxxsm(tra,ctra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,& + & a%irp,a%ja,a%val,x,size(x,1,kind=psb_ipk_),y,size(y,1,kind=psb_ipk_),info) do i = 1, m y(i,1:nc) = alpha*y(i,1:nc) end do @@ -898,8 +1088,8 @@ subroutine psb_d_cxx_cssm(alpha,a,x,beta,y,info,trans) goto 9999 end if - call inner_cxxsm(tra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,& - & a%irp,a%ja,a%val,x,size(x,1),tmp,size(tmp,1),info) + call inner_cxxsm(tra,ctra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,& + & a%irp,a%ja,a%val,x,size(x,1,kind=psb_ipk_),tmp,size(tmp,1,kind=psb_ipk_),info) do i = 1, m y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc) end do @@ -927,14 +1117,14 @@ subroutine psb_d_cxx_cssm(alpha,a,x,beta,y,info,trans) contains - subroutine inner_cxxsm(tra,lower,unit,nr,nc,& + subroutine inner_cxxsm(tra,ctra,lower,unit,nr,nc,& & irp,ja,val,x,ldx,y,ldy,info) implicit none - logical, intent(in) :: tra,lower,unit - integer(psb_ipk_), intent(in) :: nr,nc,ldx,ldy,irp(*),ja(*) - real(psb_dpk_), intent(in) :: val(*), x(ldx,*) - real(psb_dpk_), intent(out) :: y(ldy,*) - integer(psb_ipk_), intent(out) :: info + logical, intent(in) :: tra,ctra,lower,unit + integer(psb_ipk_), intent(in) :: nr,nc,ldx,ldy,irp(*),ja(*) + real(psb_dpk_), intent(in) :: val(*), x(ldx,*) + real(psb_dpk_), intent(out) :: y(ldy,*) + integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: i,j,k,m, ir, jc real(psb_dpk_), allocatable :: acc(:) @@ -946,8 +1136,7 @@ contains end if - if (.not.tra) then - + if ((.not.tra).and.(.not.ctra)) then if (lower) then if (unit) then do i=1, nr @@ -1035,13 +1224,82 @@ contains end if end if + + else if (ctra) then + + do i=1, nr + y(i,1:nc) = x(i,1:nc) + end do + + if (lower) then + if (unit) then + do i=nr, 1, -1 + acc = y(i,1:nc) + do j=irp(i), irp(i+1)-1 + jc = ja(j) + y(jc,1:nc) = y(jc,1:nc) - (val(j))*acc + end do + end do + else if (.not.unit) then + do i=nr, 1, -1 + y(i,1:nc) = y(i,1:nc)/(val(irp(i+1)-1)) + acc = y(i,1:nc) + do j=irp(i), irp(i+1)-2 + jc = ja(j) + y(jc,1:nc) = y(jc,1:nc) - (val(j))*acc + end do + end do + end if + else if (.not.lower) then + + if (unit) then + do i=1, nr + acc = y(i,1:nc) + do j=irp(i), irp(i+1)-1 + jc = ja(j) + y(jc,1:nc) = y(jc,1:nc) - (val(j))*acc + end do + end do + else if (.not.unit) then + do i=1, nr + y(i,1:nc) = y(i,1:nc)/(val(irp(i))) + acc = y(i,1:nc) + do j=irp(i)+1, irp(i+1)-1 + jc = ja(j) + y(jc,1:nc) = y(jc,1:nc) - (val(j))*acc + end do + end do + end if + + end if end if end subroutine inner_cxxsm end subroutine psb_d_cxx_cssm +function psb_d_cxx_maxval(a) result(res) + use psb_error_mod + use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_maxval + implicit none + class(psb_d_cxx_sparse_mat), intent(in) :: a + real(psb_dpk_) :: res + + integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc, nc, info + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='d_cxx_maxval' + logical, parameter :: debug=.false. + + + res = dzero + nnz = a%get_nzeros() + if (allocated(a%val)) then + nnz = min(nnz,size(a%val)) + res = maxval(abs(a%val(1:nnz))) + end if +end function psb_d_cxx_maxval + function psb_d_cxx_csnmi(a) result(res) - use psb_base_mod + use psb_error_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_csnmi implicit none class(psb_d_cxx_sparse_mat), intent(in) :: a @@ -1051,11 +1309,12 @@ function psb_d_cxx_csnmi(a) result(res) real(psb_dpk_) :: acc logical :: tra integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_csnmi' logical, parameter :: debug=.false. - res = dzero + res = dzero do i = 1, a%get_nrows() acc = dzero @@ -1067,9 +1326,9 @@ function psb_d_cxx_csnmi(a) result(res) end function psb_d_cxx_csnmi - function psb_d_cxx_csnm1(a) result(res) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_csnm1 implicit none @@ -1081,11 +1340,12 @@ function psb_d_cxx_csnm1(a) result(res) real(psb_dpk_), allocatable :: vt(:) logical :: tra integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cxx_csnm1' logical, parameter :: debug=.false. - res = -done + res = dzero nnz = a%get_nzeros() m = a%get_nrows() n = a%get_ncols() @@ -1095,7 +1355,7 @@ function psb_d_cxx_csnm1(a) result(res) do i=1, m do j=a%irp(i),a%irp(i+1)-1 k = a%ja(j) - vt(k) = vt(k) + abs(a%val(k)) + vt(k) = vt(k) + abs(a%val(j)) end do end do res = maxval(vt(1:n)) @@ -1106,7 +1366,8 @@ function psb_d_cxx_csnm1(a) result(res) end function psb_d_cxx_csnm1 subroutine psb_d_cxx_rowsum(d,a) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_rowsum class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) @@ -1115,7 +1376,8 @@ subroutine psb_d_cxx_rowsum(d,a) real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: vt(:) logical :: tra - integer(psb_ipk_) :: err_act, info, int_err(5) + integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='rowsum' logical, parameter :: debug=.false. @@ -1124,10 +1386,8 @@ subroutine psb_d_cxx_rowsum(d,a) m = a%get_nrows() if (size(d) < m) then info=psb_err_input_asize_small_i_ - int_err(1) = 1 - int_err(2) = size(d) - int_err(3) = m - call psb_errpush(info,name,i_err=int_err) + ierr(1) = 1; ierr(2) = size(d); ierr(3) = m + call psb_errpush(info,name,i_err=ierr) goto 9999 end if @@ -1154,7 +1414,8 @@ subroutine psb_d_cxx_rowsum(d,a) end subroutine psb_d_cxx_rowsum subroutine psb_d_cxx_arwsum(d,a) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_arwsum class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) @@ -1163,7 +1424,8 @@ subroutine psb_d_cxx_arwsum(d,a) real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: vt(:) logical :: tra - integer(psb_ipk_) :: err_act, info, int_err(5) + integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='rowsum' logical, parameter :: debug=.false. @@ -1172,10 +1434,8 @@ subroutine psb_d_cxx_arwsum(d,a) m = a%get_nrows() if (size(d) < m) then info=psb_err_input_asize_small_i_ - int_err(1) = 1 - int_err(2) = size(d) - int_err(3) = m - call psb_errpush(info,name,i_err=int_err) + ierr(1) = 1; ierr(2) = size(d); ierr(3) = m + call psb_errpush(info,name,i_err=ierr) goto 9999 end if @@ -1202,7 +1462,8 @@ subroutine psb_d_cxx_arwsum(d,a) end subroutine psb_d_cxx_arwsum subroutine psb_d_cxx_colsum(d,a) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_colsum class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) @@ -1211,7 +1472,8 @@ subroutine psb_d_cxx_colsum(d,a) real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: vt(:) logical :: tra - integer(psb_ipk_) :: err_act, info, int_err(5) + integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='colsum' logical, parameter :: debug=.false. @@ -1221,10 +1483,8 @@ subroutine psb_d_cxx_colsum(d,a) n = a%get_ncols() if (size(d) < n) then info=psb_err_input_asize_small_i_ - int_err(1) = 1 - int_err(2) = size(d) - int_err(3) = n - call psb_errpush(info,name,i_err=int_err) + ierr(1) = 1; ierr(2) = size(d); ierr(3) = n + call psb_errpush(info,name,i_err=ierr) goto 9999 end if @@ -1233,7 +1493,7 @@ subroutine psb_d_cxx_colsum(d,a) do i=1, m do j=a%irp(i),a%irp(i+1)-1 k = a%ja(j) - d(k) = d(k) + (a%val(k)) + d(k) = d(k) + (a%val(j)) end do end do @@ -1253,7 +1513,8 @@ subroutine psb_d_cxx_colsum(d,a) end subroutine psb_d_cxx_colsum subroutine psb_d_cxx_aclsum(d,a) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_aclsum class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) @@ -1262,7 +1523,8 @@ subroutine psb_d_cxx_aclsum(d,a) real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: vt(:) logical :: tra - integer(psb_ipk_) :: err_act, info, int_err(5) + integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='aclsum' logical, parameter :: debug=.false. @@ -1272,10 +1534,8 @@ subroutine psb_d_cxx_aclsum(d,a) n = a%get_ncols() if (size(d) < n) then info=psb_err_input_asize_small_i_ - int_err(1) = 1 - int_err(2) = size(d) - int_err(3) = n - call psb_errpush(info,name,i_err=int_err) + ierr(1) = 1; ierr(2) = size(d); ierr(3) = n + call psb_errpush(info,name,i_err=ierr) goto 9999 end if @@ -1284,7 +1544,7 @@ subroutine psb_d_cxx_aclsum(d,a) do i=1, m do j=a%irp(i),a%irp(i+1)-1 k = a%ja(j) - d(k) = d(k) + abs(a%val(k)) + d(k) = d(k) + abs(a%val(j)) end do end do @@ -1303,9 +1563,9 @@ subroutine psb_d_cxx_aclsum(d,a) end subroutine psb_d_cxx_aclsum - subroutine psb_d_cxx_get_diag(a,d,info) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_get_diag implicit none class(psb_d_cxx_sparse_mat), intent(in) :: a @@ -1313,6 +1573,7 @@ subroutine psb_d_cxx_get_diag(a,d,info) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act, mnm, i, j, k + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='get_diag' logical, parameter :: debug=.false. @@ -1322,22 +1583,29 @@ subroutine psb_d_cxx_get_diag(a,d,info) mnm = min(a%get_nrows(),a%get_ncols()) if (size(d) < mnm) then info=psb_err_input_asize_invalid_i_ - call psb_errpush(info,name,i_err=(/2,size(d),0,0,0/)) + ierr(1) = 2; ierr(2) = size(d); + call psb_errpush(info,name,i_err=ierr) goto 9999 end if - do i=1, mnm - do k=a%irp(i),a%irp(i+1)-1 - j=a%ja(k) - if ((j == i) .and.(j <= mnm )) then - d(i) = a%val(k) - endif - enddo - end do + if (a%is_triangle().and.a%is_unit()) then + d(1:mnm) = done + else + do i=1, mnm + d(i) = dzero + do k=a%irp(i),a%irp(i+1)-1 + j=a%ja(k) + if ((j == i) .and.(j <= mnm )) then + d(i) = a%val(k) + endif + enddo + end do + end if do i=mnm+1,size(d) d(i) = dzero end do + call psb_erractionrestore(err_act) return @@ -1352,33 +1620,64 @@ subroutine psb_d_cxx_get_diag(a,d,info) end subroutine psb_d_cxx_get_diag -subroutine psb_d_cxx_scal(d,a,info) - use psb_base_mod +subroutine psb_d_cxx_scal(d,a,info,side) + use psb_error_mod + use psb_const_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_scal + use psb_string_mod implicit none class(psb_d_cxx_sparse_mat), intent(inout) :: a real(psb_dpk_), intent(in) :: d(:) integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: side integer(psb_ipk_) :: err_act,mnm, i, j, m + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='scal' + character :: side_ + logical :: left logical, parameter :: debug=.false. info = psb_success_ call psb_erractionsave(err_act) - m = a%get_nrows() - if (size(d) < m) then - info=psb_err_input_asize_invalid_i_ - call psb_errpush(info,name,i_err=(/2,size(d),0,0,0/)) - goto 9999 + side_ = 'L' + if (present(side)) then + side_ = psb_toupper(side) end if - do i=1, m - do j = a%irp(i), a%irp(i+1) -1 - a%val(j) = a%val(j) * d(i) - end do - enddo + left = (side_ == 'L') + + if (left) then + m = a%get_nrows() + if (size(d) < m) then + info=psb_err_input_asize_invalid_i_ + ierr(1) = 2; ierr(2) = size(d); + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + do i=1, m + do j = a%irp(i), a%irp(i+1) -1 + a%val(j) = a%val(j) * d(i) + end do + enddo + else + m = a%get_ncols() + if (size(d) < m) then + info=psb_err_input_asize_invalid_i_ + ierr(1) = 2; ierr(2) = size(d); + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + do i=1,a%get_nzeros() + j = a%ja(i) + a%val(i) = a%val(i) * d(j) + enddo + end if + + call psb_erractionrestore(err_act) return @@ -1395,7 +1694,8 @@ end subroutine psb_d_cxx_scal subroutine psb_d_cxx_scals(d,a,info) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_scals implicit none class(psb_d_cxx_sparse_mat), intent(inout) :: a @@ -1403,6 +1703,7 @@ subroutine psb_d_cxx_scals(d,a,info) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act,mnm, i, j, m + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='scal' logical, parameter :: debug=.false. @@ -1444,12 +1745,14 @@ end subroutine psb_d_cxx_scals subroutine psb_d_cxx_reallocate_nz(nz,a) - use psb_base_mod + use psb_error_mod + use psb_realloc_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_reallocate_nz implicit none integer(psb_ipk_), intent(in) :: nz class(psb_d_cxx_sparse_mat), intent(inout) :: a integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cxx_reallocate_nz' logical, parameter :: debug=.false. @@ -1479,13 +1782,14 @@ subroutine psb_d_cxx_reallocate_nz(nz,a) end subroutine psb_d_cxx_reallocate_nz subroutine psb_d_cxx_mold(a,b,info) - use psb_base_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_mold + use psb_error_mod implicit none class(psb_d_cxx_sparse_mat), intent(in) :: a class(psb_d_base_sparse_mat), intent(out), allocatable :: b integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='reallocate_nz' logical, parameter :: debug=.false. @@ -1508,13 +1812,15 @@ subroutine psb_d_cxx_mold(a,b,info) end subroutine psb_d_cxx_mold subroutine psb_d_cxx_allocate_mnnz(m,n,a,nz) - use psb_base_mod + use psb_error_mod + use psb_realloc_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_allocate_mnnz implicit none integer(psb_ipk_), intent(in) :: m,n class(psb_d_cxx_sparse_mat), intent(inout) :: a integer(psb_ipk_), intent(in), optional :: nz integer(psb_ipk_) :: err_act, info, nz_ + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='allocate_mnz' logical, parameter :: debug=.false. @@ -1522,12 +1828,14 @@ subroutine psb_d_cxx_allocate_mnnz(m,n,a,nz) info = psb_success_ if (m < 0) then info = psb_err_iarg_neg_ - call psb_errpush(info,name,i_err=(/1,0,0,0,0/)) + ierr(1) = ione; ierr(2) = izero; + call psb_errpush(info,name,i_err=ierr) goto 9999 endif if (n < 0) then info = psb_err_iarg_neg_ - call psb_errpush(info,name,i_err=(/2,0,0,0,0/)) + ierr(1) = 2; ierr(2) = izero; + call psb_errpush(info,name,i_err=ierr) goto 9999 endif if (present(nz)) then @@ -1537,7 +1845,8 @@ subroutine psb_d_cxx_allocate_mnnz(m,n,a,nz) end if if (nz_ < 0) then info = psb_err_iarg_neg_ - call psb_errpush(info,name,i_err=(/3,0,0,0,0/)) + ierr(1) = 3; ierr(2) = izero; + call psb_errpush(info,name,i_err=ierr) goto 9999 endif @@ -1571,7 +1880,11 @@ end subroutine psb_d_cxx_allocate_mnnz subroutine psb_d_cxx_csgetptn(imin,imax,a,nz,ia,ja,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) - use psb_base_mod + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_error_mod + use psb_d_base_mat_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_csgetptn implicit none @@ -1587,6 +1900,7 @@ subroutine psb_d_cxx_csgetptn(imin,imax,a,nz,ia,ja,info,& logical :: append_, rscale_, cscale_ integer(psb_ipk_) :: nzin_, jmin_, jmax_, err_act, i + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='csget' logical, parameter :: debug=.false. @@ -1667,6 +1981,10 @@ contains subroutine cxx_getptn(imin,imax,jmin,jmax,a,nz,ia,ja,nzin,append,info,& & iren) + use psb_const_mod + use psb_error_mod + use psb_realloc_mod + use psb_sort_mod implicit none class(psb_d_cxx_sparse_mat), intent(in) :: a @@ -1738,7 +2056,11 @@ end subroutine psb_d_cxx_csgetptn subroutine psb_d_cxx_csgetrow(imin,imax,a,nz,ia,ja,val,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) - use psb_base_mod + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_error_mod + use psb_d_base_mat_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_csgetrow implicit none @@ -1755,6 +2077,7 @@ subroutine psb_d_cxx_csgetrow(imin,imax,a,nz,ia,ja,val,info,& logical :: append_, rscale_, cscale_ integer(psb_ipk_) :: nzin_, jmin_, jmax_, err_act, i + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='csget' logical, parameter :: debug=.false. @@ -1836,6 +2159,10 @@ contains subroutine cxx_getrow(imin,imax,jmin,jmax,a,nz,ia,ja,val,nzin,append,info,& & iren) + use psb_const_mod + use psb_error_mod + use psb_realloc_mod + use psb_sort_mod implicit none class(psb_d_cxx_sparse_mat), intent(in) :: a @@ -1910,7 +2237,9 @@ end subroutine psb_d_cxx_csgetrow subroutine psb_d_cxx_csgetblk(imin,imax,a,b,info,& & jmin,jmax,iren,append,rscale,cscale) - use psb_base_mod + ! Output is always in COO format + use psb_error_mod + use psb_const_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_csgetblk implicit none @@ -1923,6 +2252,7 @@ subroutine psb_d_cxx_csgetblk(imin,imax,a,b,info,& integer(psb_ipk_), intent(in), optional :: jmin,jmax logical, intent(in), optional :: rscale,cscale integer(psb_ipk_) :: err_act, nzin, nzout + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='csget' logical :: append_ logical, parameter :: debug=.false. @@ -1968,7 +2298,8 @@ end subroutine psb_d_cxx_csgetblk subroutine psb_d_cxx_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) - use psb_base_mod + use psb_error_mod + use psb_realloc_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_csput implicit none @@ -1980,9 +2311,10 @@ subroutine psb_d_cxx_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cxx_csput' logical, parameter :: debug=.false. - integer(psb_ipk_) :: nza, i,j,k, nzl, isza, int_err(5) + integer(psb_ipk_) :: nza, i,j,k, nzl, isza call psb_erractionsave(err_act) @@ -1990,27 +2322,27 @@ subroutine psb_d_cxx_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) if (nz <= 0) then info = psb_err_iarg_neg_ - int_err(1)=1 - call psb_errpush(info,name,i_err=int_err) + ierr(1)=1 + call psb_errpush(info,name,i_err=ierr) goto 9999 end if if (size(ia) < nz) then info = psb_err_input_asize_invalid_i_ - int_err(1)=2 - call psb_errpush(info,name,i_err=int_err) + ierr(1)=2 + call psb_errpush(info,name,i_err=ierr) goto 9999 end if if (size(ja) < nz) then info = psb_err_input_asize_invalid_i_ - int_err(1)=3 - call psb_errpush(info,name,i_err=int_err) + ierr(1)=3 + call psb_errpush(info,name,i_err=ierr) goto 9999 end if if (size(val) < nz) then info = psb_err_input_asize_invalid_i_ - int_err(1)=4 - call psb_errpush(info,name,i_err=int_err) + ierr(1)=4 + call psb_errpush(info,name,i_err=ierr) goto 9999 end if @@ -2058,6 +2390,10 @@ contains subroutine psb_d_cxx_srch_upd(nz,ia,ja,val,a,& & imin,imax,jmin,jmax,info,gtl) + use psb_const_mod + use psb_realloc_mod + use psb_string_mod + use psb_sort_mod implicit none class(psb_d_cxx_sparse_mat), intent(inout) :: a @@ -2250,7 +2586,7 @@ end subroutine psb_d_cxx_csput subroutine psb_d_cxx_reinit(a,clear) - use psb_base_mod + use psb_error_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_reinit implicit none @@ -2258,6 +2594,7 @@ subroutine psb_d_cxx_reinit(a,clear) logical, intent(in), optional :: clear integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='reinit' logical :: clear_ logical, parameter :: debug=.false. @@ -2299,11 +2636,13 @@ subroutine psb_d_cxx_reinit(a,clear) end subroutine psb_d_cxx_reinit subroutine psb_d_cxx_trim(a) - use psb_base_mod + use psb_realloc_mod + use psb_error_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_trim implicit none class(psb_d_cxx_sparse_mat), intent(inout) :: a integer(psb_ipk_) :: err_act, info, nz, m + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='trim' logical, parameter :: debug=.false. @@ -2332,7 +2671,7 @@ subroutine psb_d_cxx_trim(a) end subroutine psb_d_cxx_trim subroutine psb_d_cxx_print(iout,a,iv,head,ivr,ivc) - use psb_base_mod + use psb_string_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_print implicit none @@ -2343,9 +2682,10 @@ subroutine psb_d_cxx_print(iout,a,iv,head,ivr,ivc) integer(psb_ipk_), intent(in), optional :: ivr(:), ivc(:) integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cxx_print' logical, parameter :: debug=.false. - + character(len=*), parameter :: datatype='real' character(len=80) :: frmtv integer(psb_ipk_) :: irs,ics,i,j, nmx, ni, nr, nc, nz @@ -2362,7 +2702,11 @@ subroutine psb_d_cxx_print(iout,a,iv,head,ivr,ivc) nmx = max(nr,nc,1) ni = floor(log10(1.0*nmx)) + 1 - write(frmtv,'(a,i3.3,a,i3.3,a)') '(2(i',ni,',1x),es26.18,1x,2(i',ni,',1x))' + if (datatype=='real') then + write(frmtv,'(a,i3.3,a,i3.3,a)') '(2(i',ni,',1x),es26.18,1x,2(i',ni,',1x))' + else + write(frmtv,'(a,i3.3,a,i3.3,a)') '(2(i',ni,',1x),2(es26.18,1x),2(i',ni,',1x))' + end if write(iout,*) nr, nc, nz if(present(iv)) then do i=1, nr @@ -2402,7 +2746,9 @@ end subroutine psb_d_cxx_print subroutine psb_d_cp_cxx_from_coo(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_realloc_mod + use psb_d_base_mat_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cp_cxx_from_coo implicit none @@ -2414,22 +2760,23 @@ subroutine psb_d_cp_cxx_from_coo(a,b,info) integer(psb_ipk_), allocatable :: itemp(:) !locals logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, i,j,irw, idl,err_act, nc + integer(psb_ipk_) :: nza, nr, i,j,irw, err_act, nc integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name - info = psb_success_ - ! This is to have fix_coo called behind the scenes - call b%cp_to_coo(tmp,info) - if (info == psb_success_) call a%mv_from_coo(tmp,info) + info = psb_success_ + ! This is to have fix_coo called behind the scenes + call tmp%cp_from_coo(b,info) + if (info == psb_success_) call a%mv_from_coo(tmp,info) end subroutine psb_d_cp_cxx_from_coo subroutine psb_d_cp_cxx_to_coo(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_d_base_mat_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cp_cxx_to_coo implicit none @@ -2440,7 +2787,7 @@ subroutine psb_d_cp_cxx_to_coo(a,b,info) integer(psb_ipk_), allocatable :: itemp(:) !locals logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, nc,i,j,irw, idl,err_act + integer(psb_ipk_) :: nza, nr, nc,i,j,irw, err_act integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name @@ -2469,7 +2816,9 @@ end subroutine psb_d_cp_cxx_to_coo subroutine psb_d_mv_cxx_to_coo(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_realloc_mod + use psb_d_base_mat_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_mv_cxx_to_coo implicit none @@ -2480,7 +2829,7 @@ subroutine psb_d_mv_cxx_to_coo(a,b,info) integer(psb_ipk_), allocatable :: itemp(:) !locals logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, nc,i,j,irw, idl,err_act + integer(psb_ipk_) :: nza, nr, nc,i,j,irw, err_act integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name @@ -2511,7 +2860,10 @@ end subroutine psb_d_mv_cxx_to_coo subroutine psb_d_mv_cxx_from_coo(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_realloc_mod + use psb_error_mod + use psb_d_base_mat_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_mv_cxx_from_coo implicit none @@ -2522,18 +2874,23 @@ subroutine psb_d_mv_cxx_from_coo(a,b,info) integer(psb_ipk_), allocatable :: itemp(:) !locals logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, i,j,irw, idl,err_act, nc + integer(psb_ipk_) :: nza, nr, i,j,irw, err_act, nc integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit - character(len=20) :: name + character(len=20) :: name='mv_from_coo' info = psb_success_ + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + call b%fix(info) if (info /= psb_success_) return + nr = b%get_nrows() nc = b%get_ncols() nza = b%get_nzeros() + call a%psb_d_base_sparse_mat%mv_from(b%psb_d_base_sparse_mat) ! Dirty trick: call move_alloc to have the new data allocated just once. @@ -2541,9 +2898,8 @@ subroutine psb_d_mv_cxx_from_coo(a,b,info) call move_alloc(b%ja,a%ja) call move_alloc(b%val,a%val) call psb_realloc(max(nr+1,nc+1),a%irp,info) - call b%free() - if (info /= psb_success_) return + if (nza <= 0) then a%irp(:) = 1 else @@ -2563,7 +2919,7 @@ subroutine psb_d_mv_cxx_from_coo(a,b,info) if (i >= irw) exit inner if (i>nr) then write(debug_unit,*) trim(name),& - & 'Strange situation: i>nr ',i,nr,j,nza,irw,idl + & 'Strange situation: i>nr ',i,nr,j,nza,irw exit outer end if a%irp(i+1) = a%irp(i) @@ -2598,7 +2954,8 @@ end subroutine psb_d_mv_cxx_from_coo subroutine psb_d_mv_cxx_to_fmt(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_d_base_mat_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_mv_cxx_to_fmt implicit none @@ -2609,7 +2966,7 @@ subroutine psb_d_mv_cxx_to_fmt(a,b,info) !locals type(psb_d_coo_sparse_mat) :: tmp logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, i,j,irw, idl,err_act, nc + integer(psb_ipk_) :: nza, nr, i,j,irw, err_act, nc integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name @@ -2636,7 +2993,9 @@ end subroutine psb_d_mv_cxx_to_fmt subroutine psb_d_cp_cxx_to_fmt(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_d_base_mat_mod + use psb_realloc_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cp_cxx_to_fmt implicit none @@ -2647,7 +3006,7 @@ subroutine psb_d_cp_cxx_to_fmt(a,b,info) !locals type(psb_d_coo_sparse_mat) :: tmp logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, i,j,irw, idl,err_act, nc + integer(psb_ipk_) :: nza, nr, i,j,irw, err_act, nc integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name @@ -2674,7 +3033,8 @@ end subroutine psb_d_cp_cxx_to_fmt subroutine psb_d_mv_cxx_from_fmt(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_d_base_mat_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_mv_cxx_from_fmt implicit none @@ -2685,7 +3045,7 @@ subroutine psb_d_mv_cxx_from_fmt(a,b,info) !locals type(psb_d_coo_sparse_mat) :: tmp logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, i,j,irw, idl,err_act, nc + integer(psb_ipk_) :: nza, nr, i,j,irw, err_act, nc integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name @@ -2713,7 +3073,9 @@ end subroutine psb_d_mv_cxx_from_fmt subroutine psb_d_cp_cxx_from_fmt(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_d_base_mat_mod + use psb_realloc_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cp_cxx_from_fmt implicit none @@ -2724,7 +3086,7 @@ subroutine psb_d_cp_cxx_from_fmt(a,b,info) !locals type(psb_d_coo_sparse_mat) :: tmp logical :: rwshr_ - integer(psb_ipk_) :: nz, nr, i,j,irw, idl,err_act, nc + integer(psb_ipk_) :: nz, nr, i,j,irw, err_act, nc integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name @@ -2749,7 +3111,8 @@ end subroutine psb_d_cp_cxx_from_fmt subroutine psb_d_cxx_cp_from(a,b) - use psb_base_mod + use psb_error_mod + use psb_realloc_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_cp_from implicit none @@ -2758,6 +3121,7 @@ subroutine psb_d_cxx_cp_from(a,b) integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='cp_from' logical, parameter :: debug=.false. @@ -2767,9 +3131,9 @@ subroutine psb_d_cxx_cp_from(a,b) call a%allocate(b%get_nrows(),b%get_ncols(),b%get_nzeros()) call a%psb_d_base_sparse_mat%cp_from(b%psb_d_base_sparse_mat) - a%irp = b%irp - a%ja = b%ja - a%val = b%val + call psb_safe_cpy( b%irp, a%irp , info) + call psb_safe_cpy( b%ja , a%ja , info) + call psb_safe_cpy( b%val, a%val , info) if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) @@ -2788,7 +3152,7 @@ subroutine psb_d_cxx_cp_from(a,b) end subroutine psb_d_cxx_cp_from subroutine psb_d_cxx_mv_from(a,b) - use psb_base_mod + use psb_error_mod use psb_d_cxx_mat_mod, psb_protect_name => psb_d_cxx_mv_from implicit none @@ -2797,6 +3161,7 @@ subroutine psb_d_cxx_mv_from(a,b) integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='mv_from' logical, parameter :: debug=.false. diff --git a/test/serial/psb_d_cxx_mat_mod.f90 b/test/serial/psb_d_cxx_mat_mod.f90 index 9d4bca6a6..23076abbd 100644 --- a/test/serial/psb_d_cxx_mat_mod.f90 +++ b/test/serial/psb_d_cxx_mat_mod.f90 @@ -1,10 +1,58 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.0 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ +! +! package: psb_d_cxx_mat_mod +! +! This module contains the definition of the psb_d_cxx_sparse_mat type +! which implements an actual storage format. It is a fictitious example, +! identical to CSR except in name, to show how to extend the library. +! module psb_d_cxx_mat_mod use psb_d_base_mat_mod + !> \namespace psb_base_mod \class psb_d_cxx_sparse_mat + !! \extends psb_d_base_mat_mod::psb_d_base_sparse_mat + !! + !! psb_d_cxx_sparse_mat type and the related methods. + !! This is a very common storage type, and is the default for assembled + !! matrices in our library type, extends(psb_d_base_sparse_mat) :: psb_d_cxx_sparse_mat - integer(psb_ipk_), allocatable :: irp(:), ja(:) + !> Pointers to beginning of rows in JA and VAL. + integer(psb_ipk_), allocatable :: irp(:) + !> Column indices. + integer(psb_ipk_), allocatable :: ja(:) + !> Coefficient values. real(psb_dpk_), allocatable :: val(:) contains @@ -18,6 +66,7 @@ module psb_d_cxx_mat_mod procedure, pass(a) :: d_inner_cssv => psb_d_cxx_cssv procedure, pass(a) :: d_scals => psb_d_cxx_scals procedure, pass(a) :: d_scal => psb_d_cxx_scal + procedure, pass(a) :: maxval => psb_d_cxx_maxval procedure, pass(a) :: csnmi => psb_d_cxx_csnmi procedure, pass(a) :: csnm1 => psb_d_cxx_csnm1 procedure, pass(a) :: rowsum => psb_d_cxx_rowsum @@ -54,50 +103,65 @@ module psb_d_cxx_mat_mod private :: d_cxx_get_nzeros, d_cxx_free, d_cxx_get_fmt, & & d_cxx_get_size, d_cxx_sizeof, d_cxx_get_nz_row + !> \memberof psb_d_cxx_sparse_mat + !| \see psb_base_mat_mod::psb_base_reallocate_nz interface subroutine psb_d_cxx_reallocate_nz(nz,a) - import :: psb_d_cxx_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat integer(psb_ipk_), intent(in) :: nz class(psb_d_cxx_sparse_mat), intent(inout) :: a end subroutine psb_d_cxx_reallocate_nz end interface + !> \memberof psb_d_cxx_sparse_mat + !| \see psb_base_mat_mod::psb_base_reinit interface subroutine psb_d_cxx_reinit(a,clear) - import :: psb_d_cxx_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat class(psb_d_cxx_sparse_mat), intent(inout) :: a logical, intent(in), optional :: clear end subroutine psb_d_cxx_reinit end interface + !> \memberof psb_d_cxx_sparse_mat + !| \see psb_base_mat_mod::psb_base_trim interface subroutine psb_d_cxx_trim(a) - import :: psb_d_cxx_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat class(psb_d_cxx_sparse_mat), intent(inout) :: a end subroutine psb_d_cxx_trim end interface + + !> \memberof psb_d_cxx_sparse_mat + !| \see psb_base_mat_mod::psb_base_mold interface subroutine psb_d_cxx_mold(a,b,info) - import :: psb_d_cxx_sparse_mat, psb_d_base_sparse_mat, psb_long_int_k_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_d_base_sparse_mat, psb_long_int_k_ class(psb_d_cxx_sparse_mat), intent(in) :: a class(psb_d_base_sparse_mat), intent(out), allocatable :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_cxx_mold end interface + + !> \memberof psb_d_cxx_sparse_mat + !| \see psb_base_mat_mod::psb_base_allocate_mnnz interface subroutine psb_d_cxx_allocate_mnnz(m,n,a,nz) - import :: psb_d_cxx_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat integer(psb_ipk_), intent(in) :: m,n class(psb_d_cxx_sparse_mat), intent(inout) :: a integer(psb_ipk_), intent(in), optional :: nz end subroutine psb_d_cxx_allocate_mnnz end interface + + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_print interface subroutine psb_d_cxx_print(iout,a,iv,head,ivr,ivc) - import :: psb_d_cxx_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat integer(psb_ipk_), intent(in) :: iout class(psb_d_cxx_sparse_mat), intent(in) :: a integer(psb_ipk_), intent(in), optional :: iv(:) @@ -106,98 +170,120 @@ module psb_d_cxx_mat_mod end subroutine psb_d_cxx_print end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cp_to_coo interface subroutine psb_d_cp_cxx_to_coo(a,b,info) - import :: psb_d_coo_sparse_mat, psb_d_cxx_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_coo_sparse_mat, psb_d_cxx_sparse_mat class(psb_d_cxx_sparse_mat), intent(in) :: a class(psb_d_coo_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_cp_cxx_to_coo end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cp_from_coo interface subroutine psb_d_cp_cxx_from_coo(a,b,info) - import :: psb_d_cxx_sparse_mat, psb_d_coo_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_d_coo_sparse_mat class(psb_d_cxx_sparse_mat), intent(inout) :: a class(psb_d_coo_sparse_mat), intent(in) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_cp_cxx_from_coo end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cp_to_fmt interface subroutine psb_d_cp_cxx_to_fmt(a,b,info) - import :: psb_d_cxx_sparse_mat, psb_d_base_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_d_base_sparse_mat class(psb_d_cxx_sparse_mat), intent(in) :: a class(psb_d_base_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_cp_cxx_to_fmt end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cp_from_fmt interface subroutine psb_d_cp_cxx_from_fmt(a,b,info) - import :: psb_d_cxx_sparse_mat, psb_d_base_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_d_base_sparse_mat class(psb_d_cxx_sparse_mat), intent(inout) :: a class(psb_d_base_sparse_mat), intent(in) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_cp_cxx_from_fmt end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_mv_to_coo interface subroutine psb_d_mv_cxx_to_coo(a,b,info) - import :: psb_d_cxx_sparse_mat, psb_d_coo_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_d_coo_sparse_mat class(psb_d_cxx_sparse_mat), intent(inout) :: a class(psb_d_coo_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_mv_cxx_to_coo end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_mv_from_coo interface subroutine psb_d_mv_cxx_from_coo(a,b,info) - import :: psb_d_cxx_sparse_mat, psb_d_coo_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_d_coo_sparse_mat class(psb_d_cxx_sparse_mat), intent(inout) :: a class(psb_d_coo_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_mv_cxx_from_coo end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_mv_to_fmt interface subroutine psb_d_mv_cxx_to_fmt(a,b,info) - import :: psb_d_cxx_sparse_mat, psb_d_base_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_d_base_sparse_mat class(psb_d_cxx_sparse_mat), intent(inout) :: a class(psb_d_base_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_mv_cxx_to_fmt end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_mv_from_fmt interface subroutine psb_d_mv_cxx_from_fmt(a,b,info) - import :: psb_d_cxx_sparse_mat, psb_d_base_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_d_base_sparse_mat class(psb_d_cxx_sparse_mat), intent(inout) :: a class(psb_d_base_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_mv_cxx_from_fmt end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cp_from interface subroutine psb_d_cxx_cp_from(a,b) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(inout) :: a type(psb_d_cxx_sparse_mat), intent(in) :: b end subroutine psb_d_cxx_cp_from end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_mv_from interface subroutine psb_d_cxx_mv_from(a,b) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(inout) :: a type(psb_d_cxx_sparse_mat), intent(inout) :: b end subroutine psb_d_cxx_mv_from end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csput interface subroutine psb_d_cxx_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(inout) :: a real(psb_dpk_), intent(in) :: val(:) integer(psb_ipk_), intent(in) :: nz,ia(:), ja(:),& @@ -207,10 +293,12 @@ module psb_d_cxx_mat_mod end subroutine psb_d_cxx_csput end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_base_mat_mod::psb_base_csgetptn interface subroutine psb_d_cxx_csgetptn(imin,imax,a,nz,ia,ja,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a integer(psb_ipk_), intent(in) :: imin,imax integer(psb_ipk_), intent(out) :: nz @@ -223,10 +311,12 @@ module psb_d_cxx_mat_mod end subroutine psb_d_cxx_csgetptn end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csgetrow interface subroutine psb_d_cxx_csgetrow(imin,imax,a,nz,ia,ja,val,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a integer(psb_ipk_), intent(in) :: imin,imax integer(psb_ipk_), intent(out) :: nz @@ -240,10 +330,12 @@ module psb_d_cxx_mat_mod end subroutine psb_d_cxx_csgetrow end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csgetblk interface subroutine psb_d_cxx_csgetblk(imin,imax,a,b,info,& & jmin,jmax,iren,append,rscale,cscale) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_d_coo_sparse_mat, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_, psb_d_coo_sparse_mat class(psb_d_cxx_sparse_mat), intent(in) :: a class(psb_d_coo_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(in) :: imin,imax @@ -255,17 +347,23 @@ module psb_d_cxx_mat_mod end subroutine psb_d_cxx_csgetblk end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cssv interface subroutine psb_d_cxx_cssv(alpha,a,x,beta,y,info,trans) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(in) :: alpha, beta, x(:) real(psb_dpk_), intent(inout) :: y(:) integer(psb_ipk_), intent(out) :: info character, optional, intent(in) :: trans end subroutine psb_d_cxx_cssv + end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cssm + interface subroutine psb_d_cxx_cssm(alpha,a,x,beta,y,info,trans) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) real(psb_dpk_), intent(inout) :: y(:,:) @@ -274,17 +372,24 @@ module psb_d_cxx_mat_mod end subroutine psb_d_cxx_cssm end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csmv interface subroutine psb_d_cxx_csmv(alpha,a,x,beta,y,info,trans) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(in) :: alpha, beta, x(:) real(psb_dpk_), intent(inout) :: y(:) integer(psb_ipk_), intent(out) :: info character, optional, intent(in) :: trans end subroutine psb_d_cxx_csmv + end interface + + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csmm + interface subroutine psb_d_cxx_csmm(alpha,a,x,beta,y,info,trans) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) real(psb_dpk_), intent(inout) :: y(:,:) @@ -294,75 +399,104 @@ module psb_d_cxx_mat_mod end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_maxval + interface + function psb_d_cxx_maxval(a) result(res) + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ + class(psb_d_cxx_sparse_mat), intent(in) :: a + real(psb_dpk_) :: res + end function psb_d_cxx_maxval + end interface + + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csnmi interface function psb_d_cxx_csnmi(a) result(res) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_) :: res end function psb_d_cxx_csnmi end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csnm1 interface function psb_d_cxx_csnm1(a) result(res) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_) :: res end function psb_d_cxx_csnm1 end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_rowsum interface subroutine psb_d_cxx_rowsum(d,a) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) end subroutine psb_d_cxx_rowsum end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_arwsum interface subroutine psb_d_cxx_arwsum(d,a) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) end subroutine psb_d_cxx_arwsum end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_colsum interface subroutine psb_d_cxx_colsum(d,a) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) end subroutine psb_d_cxx_colsum end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_aclsum interface subroutine psb_d_cxx_aclsum(d,a) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) end subroutine psb_d_cxx_aclsum end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_get_diag interface subroutine psb_d_cxx_get_diag(a,d,info) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) integer(psb_ipk_), intent(out) :: info end subroutine psb_d_cxx_get_diag end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_scal interface - subroutine psb_d_cxx_scal(d,a,info) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + subroutine psb_d_cxx_scal(d,a,info,side) + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(inout) :: a real(psb_dpk_), intent(in) :: d(:) integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: side end subroutine psb_d_cxx_scal end interface + !> \memberof psb_d_cxx_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_scals interface subroutine psb_d_cxx_scals(d,a,info) - import :: psb_d_cxx_sparse_mat, psb_dpk_, psb_ipk_ + import :: psb_ipk_, psb_d_cxx_sparse_mat, psb_dpk_ class(psb_d_cxx_sparse_mat), intent(inout) :: a real(psb_dpk_), intent(in) :: d integer(psb_ipk_), intent(out) :: info @@ -418,11 +552,7 @@ contains res = -1 if (allocated(a%ja)) then - if (res >= 0) then - res = min(res,size(a%ja)) - else - res = size(a%ja) - end if + res = size(a%ja) end if if (allocated(a%val)) then if (res >= 0) then @@ -475,8 +605,8 @@ contains if (allocated(a%ja)) deallocate(a%ja) if (allocated(a%val)) deallocate(a%val) call a%set_null() - call a%set_nrows(0) - call a%set_ncols(0) + call a%set_nrows(izero) + call a%set_ncols(izero) return diff --git a/test/serial/psb_d_cyy_impl.f90 b/test/serial/psb_d_cyy_impl.f90 index 07867bb13..1e3a5662b 100644 --- a/test/serial/psb_d_cyy_impl.f90 +++ b/test/serial/psb_d_cyy_impl.f90 @@ -1,8 +1,51 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.0 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ + +! == =================================== ! -! A minimal, non functional implementation of a matrix type module. -! +! +! +! Computational routines +! +! +! +! +! +! +! == =================================== + subroutine psb_d_cyy_csmv(alpha,a,x,beta,y,info,trans) - use psb_base_mod + use psb_error_mod + use psb_string_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_csmv implicit none class(psb_d_cyy_sparse_mat), intent(in) :: a @@ -14,18 +57,344 @@ subroutine psb_d_cyy_csmv(alpha,a,x,beta,y,info,trans) character :: trans_ integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc real(psb_dpk_) :: acc - logical :: tra + logical :: tra, ctra integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cyy_csmv' logical, parameter :: debug=.false. call psb_erractionsave(err_act) info = psb_success_ + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + + if (.not.a%is_asb()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + tra = (psb_toupper(trans_) == 'T') + ctra = (psb_toupper(trans_) == 'C') + + if (tra.or.ctra) then + m = a%get_ncols() + n = a%get_nrows() + else + n = a%get_ncols() + m = a%get_nrows() + end if + + if (size(x,1) psb_d_cyy_csmm implicit none class(psb_d_cyy_sparse_mat), intent(in) :: a @@ -37,16 +406,347 @@ subroutine psb_d_cyy_csmm(alpha,a,x,beta,y,info,trans) character :: trans_ integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc, nc real(psb_dpk_), allocatable :: acc(:) - logical :: tra + logical :: tra, ctra integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cyy_csmm' logical, parameter :: debug=.false. info = psb_success_ + call psb_erractionsave(err_act) + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + + if (.not.a%is_asb()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + tra = (psb_toupper(trans_) == 'T') + ctra = (psb_toupper(trans_) == 'C') + + if (tra.or.ctra) then + m = a%get_ncols() + n = a%get_nrows() + else + n = a%get_ncols() + m = a%get_nrows() + end if + + if (size(x,1) psb_d_cyy_cssv implicit none class(psb_d_cyy_sparse_mat), intent(in) :: a @@ -59,16 +759,262 @@ subroutine psb_d_cyy_cssv(alpha,a,x,beta,y,info,trans) integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: tmp(:) - logical :: tra + logical :: tra,ctra integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cyy_cssv' logical, parameter :: debug=.false. info = psb_success_ + call psb_erractionsave(err_act) + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + if (.not.a%is_asb()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + tra = (psb_toupper(trans_) == 'T') + ctra = (psb_toupper(trans_) == 'C') + m = a%get_nrows() + + if (.not. (a%is_triangle())) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + end if + + if (size(x) psb_d_cyy_cssm implicit none class(psb_d_cyy_sparse_mat), intent(in) :: a @@ -81,17 +1027,279 @@ subroutine psb_d_cyy_cssm(alpha,a,x,beta,y,info,trans) integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc, nc real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: tmp(:,:) - logical :: tra + logical :: tra, ctra integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cyy_cssm' logical, parameter :: debug=.false. info = psb_success_ + call psb_erractionsave(err_act) + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + if (.not.a%is_asb()) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + endif + + + tra = (psb_toupper(trans_) == 'T') + ctra = (psb_toupper(trans_) == 'C') + + m = a%get_nrows() + nc = min(size(x,2) , size(y,2)) + + if (.not. (a%is_triangle())) then + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + 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_cyysm(tra,ctra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,& + & a%irp,a%ja,a%val,x,size(x,1,kind=psb_ipk_),y,size(y,1,kind=psb_ipk_),info) + do i = 1, m + y(i,1:nc) = alpha*y(i,1:nc) + end do + else + allocate(tmp(m,nc), stat=info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='allocate') + goto 9999 + end if + + call inner_cyysm(tra,ctra,a%is_lower(),a%is_unit(),a%get_nrows(),nc,& + & a%irp,a%ja,a%val,x,size(x,1,kind=psb_ipk_),tmp,size(tmp,1,kind=psb_ipk_),info) + do i = 1, m + y(i,1:nc) = alpha*tmp(i,1:nc) + beta*y(i,1:nc) + end do + end if + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='inner_cyysm') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + +contains + + subroutine inner_cyysm(tra,ctra,lower,unit,nr,nc,& + & irp,ja,val,x,ldx,y,ldy,info) + implicit none + logical, intent(in) :: tra,ctra,lower,unit + integer(psb_ipk_), intent(in) :: nr,nc,ldx,ldy,irp(*),ja(*) + real(psb_dpk_), intent(in) :: val(*), x(ldx,*) + real(psb_dpk_), intent(out) :: y(ldy,*) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: i,j,k,m, ir, jc + real(psb_dpk_), allocatable :: acc(:) + + info = psb_success_ + allocate(acc(nc), stat=info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + return + end if + + + if ((.not.tra).and.(.not.ctra)) then + if (lower) then + if (unit) then + do i=1, nr + acc = dzero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j)*y(ja(j),1:nc) + end do + y(i,1:nc) = x(i,1:nc) - acc + end do + else if (.not.unit) then + do i=1, nr + acc = dzero + do j=irp(i), irp(i+1)-2 + acc = acc + val(j)*y(ja(j),1:nc) + end do + y(i,1:nc) = (x(i,1:nc) - acc)/val(irp(i+1)-1) + end do + end if + else if (.not.lower) then + + if (unit) then + do i=nr, 1, -1 + acc = dzero + do j=irp(i), irp(i+1)-1 + acc = acc + val(j)*y(ja(j),1:nc) + end do + y(i,1:nc) = x(i,1:nc) - acc + end do + else if (.not.unit) then + do i=nr, 1, -1 + acc = dzero + do j=irp(i)+1, irp(i+1)-1 + acc = acc + val(j)*y(ja(j),1:nc) + end do + y(i,1:nc) = (x(i,1:nc) - acc)/val(irp(i)) + end do + end if + + end if + + else if (tra) then + + do i=1, nr + y(i,1:nc) = x(i,1:nc) + end do + + if (lower) then + if (unit) then + do i=nr, 1, -1 + acc = y(i,1:nc) + do j=irp(i), irp(i+1)-1 + jc = ja(j) + y(jc,1:nc) = y(jc,1:nc) - val(j)*acc + end do + end do + else if (.not.unit) then + do i=nr, 1, -1 + y(i,1:nc) = y(i,1:nc)/val(irp(i+1)-1) + acc = y(i,1:nc) + do j=irp(i), irp(i+1)-2 + jc = ja(j) + y(jc,1:nc) = y(jc,1:nc) - val(j)*acc + end do + end do + end if + else if (.not.lower) then + + if (unit) then + do i=1, nr + acc = y(i,1:nc) + do j=irp(i), irp(i+1)-1 + jc = ja(j) + y(jc,1:nc) = y(jc,1:nc) - val(j)*acc + end do + end do + else if (.not.unit) then + do i=1, nr + y(i,1:nc) = y(i,1:nc)/val(irp(i)) + acc = y(i,1:nc) + do j=irp(i)+1, irp(i+1)-1 + jc = ja(j) + y(jc,1:nc) = y(jc,1:nc) - val(j)*acc + end do + end do + end if + + end if + + else if (ctra) then + + do i=1, nr + y(i,1:nc) = x(i,1:nc) + end do + + if (lower) then + if (unit) then + do i=nr, 1, -1 + acc = y(i,1:nc) + do j=irp(i), irp(i+1)-1 + jc = ja(j) + y(jc,1:nc) = y(jc,1:nc) - (val(j))*acc + end do + end do + else if (.not.unit) then + do i=nr, 1, -1 + y(i,1:nc) = y(i,1:nc)/(val(irp(i+1)-1)) + acc = y(i,1:nc) + do j=irp(i), irp(i+1)-2 + jc = ja(j) + y(jc,1:nc) = y(jc,1:nc) - (val(j))*acc + end do + end do + end if + else if (.not.lower) then + + if (unit) then + do i=1, nr + acc = y(i,1:nc) + do j=irp(i), irp(i+1)-1 + jc = ja(j) + y(jc,1:nc) = y(jc,1:nc) - (val(j))*acc + end do + end do + else if (.not.unit) then + do i=1, nr + y(i,1:nc) = y(i,1:nc)/(val(irp(i))) + acc = y(i,1:nc) + do j=irp(i)+1, irp(i+1)-1 + jc = ja(j) + y(jc,1:nc) = y(jc,1:nc) - (val(j))*acc + end do + end do + end if + + end if + end if + end subroutine inner_cyysm end subroutine psb_d_cyy_cssm +function psb_d_cyy_maxval(a) result(res) + use psb_error_mod + use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_maxval + implicit none + class(psb_d_cyy_sparse_mat), intent(in) :: a + real(psb_dpk_) :: res + + integer(psb_ipk_) :: i,j,k,m,n, nnz, ir, jc, nc, info + integer(psb_ipk_) :: ierr(5) + character(len=20) :: name='d_cyy_maxval' + logical, parameter :: debug=.false. + + + res = dzero + nnz = a%get_nzeros() + if (allocated(a%val)) then + nnz = min(nnz,size(a%val)) + res = maxval(abs(a%val(1:nnz))) + end if +end function psb_d_cyy_maxval + function psb_d_cyy_csnmi(a) result(res) - use psb_base_mod + use psb_error_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_csnmi implicit none class(psb_d_cyy_sparse_mat), intent(in) :: a @@ -101,13 +1309,26 @@ function psb_d_cyy_csnmi(a) result(res) real(psb_dpk_) :: acc logical :: tra integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_csnmi' logical, parameter :: debug=.false. - res = dzero + + + res = dzero + + do i = 1, a%get_nrows() + acc = dzero + do j=a%irp(i),a%irp(i+1)-1 + acc = acc + abs(a%val(j)) + end do + res = max(res,acc) + end do + end function psb_d_cyy_csnmi function psb_d_cyy_csnm1(a) result(res) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_csnm1 implicit none @@ -119,14 +1340,34 @@ function psb_d_cyy_csnm1(a) result(res) real(psb_dpk_), allocatable :: vt(:) logical :: tra integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cyy_csnm1' logical, parameter :: debug=.false. - res = -done + + + res = dzero + nnz = a%get_nzeros() + m = a%get_nrows() + n = a%get_ncols() + allocate(vt(n),stat=info) + if (info /= 0) return + vt(:) = dzero + do i=1, m + do j=a%irp(i),a%irp(i+1)-1 + k = a%ja(j) + vt(k) = vt(k) + abs(a%val(j)) + end do + end do + res = maxval(vt(1:n)) + deallocate(vt,stat=info) + return + end function psb_d_cyy_csnm1 subroutine psb_d_cyy_rowsum(d,a) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_rowsum class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) @@ -135,14 +1376,46 @@ subroutine psb_d_cyy_rowsum(d,a) real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: vt(:) logical :: tra - integer(psb_ipk_) :: err_act, info, int_err(5) + integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='rowsum' logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + m = a%get_nrows() + if (size(d) < m) then + info=psb_err_input_asize_small_i_ + ierr(1) = 1; ierr(2) = size(d); ierr(3) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + do i = 1, a%get_nrows() + d(i) = dzero + do j=a%irp(i),a%irp(i+1)-1 + d(i) = d(i) + (a%val(j)) + end do + end do + return + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine psb_d_cyy_rowsum subroutine psb_d_cyy_arwsum(d,a) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_arwsum class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) @@ -151,13 +1424,46 @@ subroutine psb_d_cyy_arwsum(d,a) real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: vt(:) logical :: tra - integer(psb_ipk_) :: err_act, info, int_err(5) + integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='rowsum' logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + m = a%get_nrows() + if (size(d) < m) then + info=psb_err_input_asize_small_i_ + ierr(1) = 1; ierr(2) = size(d); ierr(3) = m + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + + do i = 1, a%get_nrows() + d(i) = dzero + do j=a%irp(i),a%irp(i+1)-1 + d(i) = d(i) + abs(a%val(j)) + end do + end do + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine psb_d_cyy_arwsum subroutine psb_d_cyy_colsum(d,a) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_colsum class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) @@ -166,13 +1472,49 @@ subroutine psb_d_cyy_colsum(d,a) real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: vt(:) logical :: tra - integer(psb_ipk_) :: err_act, info, int_err(5) + integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='colsum' logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + m = a%get_nrows() + n = a%get_ncols() + if (size(d) < n) then + info=psb_err_input_asize_small_i_ + ierr(1) = 1; ierr(2) = size(d); ierr(3) = n + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + d = dzero + + do i=1, m + do j=a%irp(i),a%irp(i+1)-1 + k = a%ja(j) + d(k) = d(k) + (a%val(j)) + end do + end do + + return + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine psb_d_cyy_colsum subroutine psb_d_cyy_aclsum(d,a) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_aclsum class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) @@ -181,13 +1523,49 @@ subroutine psb_d_cyy_aclsum(d,a) real(psb_dpk_) :: acc real(psb_dpk_), allocatable :: vt(:) logical :: tra - integer(psb_ipk_) :: err_act, info, int_err(5) + integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='aclsum' logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + m = a%get_nrows() + n = a%get_ncols() + if (size(d) < n) then + info=psb_err_input_asize_small_i_ + ierr(1) = 1; ierr(2) = size(d); ierr(3) = n + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + d = dzero + + do i=1, m + do j=a%irp(i),a%irp(i+1)-1 + k = a%ja(j) + d(k) = d(k) + abs(a%val(j)) + end do + end do + + return + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine psb_d_cyy_aclsum subroutine psb_d_cyy_get_diag(a,d,info) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_get_diag implicit none class(psb_d_cyy_sparse_mat), intent(in) :: a @@ -195,27 +1573,129 @@ subroutine psb_d_cyy_get_diag(a,d,info) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act, mnm, i, j, k + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='get_diag' logical, parameter :: debug=.false. + info = psb_success_ + call psb_erractionsave(err_act) + + mnm = min(a%get_nrows(),a%get_ncols()) + if (size(d) < mnm) then + info=psb_err_input_asize_invalid_i_ + ierr(1) = 2; ierr(2) = size(d); + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + + if (a%is_triangle().and.a%is_unit()) then + d(1:mnm) = done + else + do i=1, mnm + d(i) = dzero + do k=a%irp(i),a%irp(i+1)-1 + j=a%ja(k) + if ((j == i) .and.(j <= mnm )) then + d(i) = a%val(k) + endif + enddo + end do + end if + do i=mnm+1,size(d) + d(i) = dzero + end do + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine psb_d_cyy_get_diag -subroutine psb_d_cyy_scal(d,a,info) - use psb_base_mod + +subroutine psb_d_cyy_scal(d,a,info,side) + use psb_error_mod + use psb_const_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_scal + use psb_string_mod implicit none class(psb_d_cyy_sparse_mat), intent(inout) :: a real(psb_dpk_), intent(in) :: d(:) integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: side integer(psb_ipk_) :: err_act,mnm, i, j, m + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='scal' + character :: side_ + logical :: left logical, parameter :: debug=.false. + info = psb_success_ + call psb_erractionsave(err_act) + + side_ = 'L' + if (present(side)) then + side_ = psb_toupper(side) + end if + + left = (side_ == 'L') + + if (left) then + m = a%get_nrows() + if (size(d) < m) then + info=psb_err_input_asize_invalid_i_ + ierr(1) = 2; ierr(2) = size(d); + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + do i=1, m + do j = a%irp(i), a%irp(i+1) -1 + a%val(j) = a%val(j) * d(i) + end do + enddo + else + m = a%get_ncols() + if (size(d) < m) then + info=psb_err_input_asize_invalid_i_ + ierr(1) = 2; ierr(2) = size(d); + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + do i=1,a%get_nzeros() + j = a%ja(i) + a%val(i) = a%val(i) * d(j) + enddo + end if + + + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine psb_d_cyy_scal + subroutine psb_d_cyy_scals(d,a,info) - use psb_base_mod + use psb_error_mod + use psb_const_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_scals implicit none class(psb_d_cyy_sparse_mat), intent(inout) :: a @@ -223,50 +1703,188 @@ subroutine psb_d_cyy_scals(d,a,info) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act,mnm, i, j, m + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='scal' logical, parameter :: debug=.false. + info = psb_success_ + call psb_erractionsave(err_act) + + + do i=1,a%get_nzeros() + a%val(i) = a%val(i) * d + enddo + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine psb_d_cyy_scals + + + +! == =================================== +! +! +! +! Data management +! +! +! +! +! +! == =================================== + + subroutine psb_d_cyy_reallocate_nz(nz,a) - use psb_base_mod + use psb_error_mod + use psb_realloc_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_reallocate_nz implicit none integer(psb_ipk_), intent(in) :: nz class(psb_d_cyy_sparse_mat), intent(inout) :: a integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cyy_reallocate_nz' logical, parameter :: debug=.false. + call psb_erractionsave(err_act) + + call psb_realloc(nz,a%ja,info) + if (info == psb_success_) call psb_realloc(nz,a%val,info) + if (info == psb_success_) call psb_realloc(& + & max(nz,a%get_nrows()+1,a%get_ncols()+1),a%irp,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine psb_d_cyy_reallocate_nz subroutine psb_d_cyy_mold(a,b,info) - use psb_base_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_mold + use psb_error_mod implicit none class(psb_d_cyy_sparse_mat), intent(in) :: a class(psb_d_base_sparse_mat), intent(out), allocatable :: b integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='reallocate_nz' logical, parameter :: debug=.false. + + call psb_get_erraction(err_act) + + allocate(psb_d_cyy_sparse_mat :: b, stat=info) + + if (info /= psb_success_) then + info = psb_err_alloc_dealloc_ + call psb_errpush(info, name) + goto 9999 + end if + return +9999 continue + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + end subroutine psb_d_cyy_mold subroutine psb_d_cyy_allocate_mnnz(m,n,a,nz) - use psb_base_mod + use psb_error_mod + use psb_realloc_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_allocate_mnnz implicit none integer(psb_ipk_), intent(in) :: m,n class(psb_d_cyy_sparse_mat), intent(inout) :: a integer(psb_ipk_), intent(in), optional :: nz integer(psb_ipk_) :: err_act, info, nz_ + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='allocate_mnz' logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = psb_success_ + if (m < 0) then + info = psb_err_iarg_neg_ + ierr(1) = ione; ierr(2) = izero; + call psb_errpush(info,name,i_err=ierr) + goto 9999 + endif + if (n < 0) then + info = psb_err_iarg_neg_ + ierr(1) = 2; ierr(2) = izero; + call psb_errpush(info,name,i_err=ierr) + goto 9999 + endif + if (present(nz)) then + nz_ = nz + else + nz_ = max(7*m,7*n,1) + end if + if (nz_ < 0) then + info = psb_err_iarg_neg_ + ierr(1) = 3; ierr(2) = izero; + call psb_errpush(info,name,i_err=ierr) + goto 9999 + endif + + if (info == psb_success_) call psb_realloc(m+1,a%irp,info) + if (info == psb_success_) call psb_realloc(nz_,a%ja,info) + if (info == psb_success_) call psb_realloc(nz_,a%val,info) + if (info == psb_success_) then + a%irp=0 + call a%set_nrows(m) + call a%set_ncols(n) + call a%set_bld() + call a%set_triangle(.false.) + call a%set_unit(.false.) + call a%set_dupl(psb_dupl_def_) + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine psb_d_cyy_allocate_mnnz + subroutine psb_d_cyy_csgetptn(imin,imax,a,nz,ia,ja,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) - use psb_base_mod + ! Output is always in COO format + use psb_error_mod + use psb_const_mod + use psb_error_mod + use psb_d_base_mat_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_csgetptn implicit none @@ -282,14 +1900,167 @@ subroutine psb_d_cyy_csgetptn(imin,imax,a,nz,ia,ja,info,& logical :: append_, rscale_, cscale_ integer(psb_ipk_) :: nzin_, jmin_, jmax_, err_act, i + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='csget' logical, parameter :: debug=.false. + call psb_erractionsave(err_act) + info = psb_success_ + + if (present(jmin)) then + jmin_ = jmin + else + jmin_ = 1 + endif + if (present(jmax)) then + jmax_ = jmax + else + jmax_ = a%get_ncols() + endif + + if ((imax psb_d_cyy_csgetrow implicit none @@ -306,16 +2077,169 @@ subroutine psb_d_cyy_csgetrow(imin,imax,a,nz,ia,ja,val,info,& logical :: append_, rscale_, cscale_ integer(psb_ipk_) :: nzin_, jmin_, jmax_, err_act, i + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='csget' logical, parameter :: debug=.false. call psb_erractionsave(err_act) info = psb_success_ + + if (present(jmin)) then + jmin_ = jmin + else + jmin_ = 1 + endif + if (present(jmax)) then + jmax_ = jmax + else + jmax_ = a%get_ncols() + endif + + if ((imax psb_d_cyy_csgetblk implicit none @@ -328,15 +2252,54 @@ subroutine psb_d_cyy_csgetblk(imin,imax,a,b,info,& integer(psb_ipk_), intent(in), optional :: jmin,jmax logical, intent(in), optional :: rscale,cscale integer(psb_ipk_) :: err_act, nzin, nzout + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='csget' logical :: append_ logical, parameter :: debug=.false. + call psb_erractionsave(err_act) info = psb_success_ + + if (present(append)) then + append_ = append + else + append_ = .false. + endif + if (append_) then + nzin = a%get_nzeros() + else + nzin = 0 + endif + + call a%csget(imin,imax,nzout,b%ia,b%ja,b%val,info,& + & jmin=jmin, jmax=jmax, iren=iren, append=append_, & + & nzin=nzin, rscale=rscale, cscale=cscale) + + if (info /= psb_success_) goto 9999 + + call b%set_nzeros(nzin+nzout) + call b%fix(info) + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine psb_d_cyy_csgetblk + + subroutine psb_d_cyy_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) - use psb_base_mod + use psb_error_mod + use psb_realloc_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_csput implicit none @@ -346,16 +2309,284 @@ subroutine psb_d_cyy_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: gtl(:) + integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cyy_csput' logical, parameter :: debug=.false. - integer(psb_ipk_) :: nza, i,j,k, nzl, isza, int_err(5) + integer(psb_ipk_) :: nza, i,j,k, nzl, isza + + + call psb_erractionsave(err_act) info = psb_success_ + if (nz <= 0) then + info = psb_err_iarg_neg_ + ierr(1)=1 + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (size(ia) < nz) then + info = psb_err_input_asize_invalid_i_ + ierr(1)=2 + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + if (size(ja) < nz) then + info = psb_err_input_asize_invalid_i_ + ierr(1)=3 + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + if (size(val) < nz) then + info = psb_err_input_asize_invalid_i_ + ierr(1)=4 + call psb_errpush(info,name,i_err=ierr) + goto 9999 + end if + + if (nz == 0) return + + nza = a%get_nzeros() + + if (a%is_bld()) then + ! Build phase should only ever be in COO + info = psb_err_invalid_mat_state_ + + else if (a%is_upd()) then + call psb_d_cyy_srch_upd(nz,ia,ja,val,a,& + & imin,imax,jmin,jmax,info,gtl) + + if (info /= psb_success_) then + + info = psb_err_invalid_mat_state_ + end if + + else + ! State is wrong. + info = psb_err_invalid_mat_state_ + end if + if (info /= psb_success_) then + call psb_errpush(info,name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + + +contains + + subroutine psb_d_cyy_srch_upd(nz,ia,ja,val,a,& + & imin,imax,jmin,jmax,info,gtl) + + use psb_const_mod + use psb_realloc_mod + use psb_string_mod + use psb_sort_mod + implicit none + + class(psb_d_cyy_sparse_mat), intent(inout) :: a + integer(psb_ipk_), intent(in) :: nz, imin,imax,jmin,jmax + integer(psb_ipk_), intent(in) :: ia(:),ja(:) + real(psb_dpk_), intent(in) :: val(:) + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: gtl(:) + integer(psb_ipk_) :: i,ir,ic, ilr, ilc, ip, & + & i1,i2,nr,nc,nnz,dupl,ng + integer(psb_ipk_) :: debug_level, debug_unit + character(len=20) :: name='d_cyy_srch_upd' + + info = psb_success_ + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + dupl = a%get_dupl() + + if (.not.a%is_sorted()) then + info = -4 + return + end if + + ilr = -1 + ilc = -1 + nnz = a%get_nzeros() + nr = a%get_nrows() + nc = a%get_ncols() + + if (present(gtl)) then + ng = size(gtl) + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if ((ir > 0).and.(ir <= nr)) then + i1 = a%irp(ir) + i2 = a%irp(ir+1) + nc=i2-i1 + + ip = psb_ibsrch(ic,nc,a%ja(i1:i2-1)) + if (ip>0) then + a%val(i1+ip-1) = val(i) + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Was searching ',ic,' in: ',i1,i2,& + & ' : ',a%ja(i1:i2-1) + info = i + return + end if + + else + + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + end if + end do + + case(psb_dupl_add_) + ! Add + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=1).and.(ir<=ng).and.(ic>=1).and.(ic<=ng)) then + ir = gtl(ir) + ic = gtl(ic) + if ((ir > 0).and.(ir <= nr)) then + i1 = a%irp(ir) + i2 = a%irp(ir+1) + nc = i2-i1 + ip = psb_ibsrch(ic,nc,a%ja(i1:i2-1)) + if (ip>0) then + a%val(i1+ip-1) = a%val(i1+ip-1) + val(i) + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Was searching ',ic,' in: ',i1,i2,& + & ' : ',a%ja(i1:i2-1) + info = i + return + end if + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + + end if + end do + + case default + info = -3 + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Duplicate handling: ',dupl + end select + + else + + select case(dupl) + case(psb_dupl_ovwrt_,psb_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + + if ((ir > 0).and.(ir <= nr)) then + + i1 = a%irp(ir) + i2 = a%irp(ir+1) + nc=i2-i1 + + ip = psb_ibsrch(ic,nc,a%ja(i1:i2-1)) + if (ip>0) then + a%val(i1+ip-1) = val(i) + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Was searching ',ic,' in: ',i1,i2,& + & ' : ',a%ja(i1:i2-1) + info = i + return + end if + + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + + end do + + case(psb_dupl_add_) + ! Add + ilr = -1 + ilc = -1 + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir > 0).and.(ir <= nr)) then + i1 = a%irp(ir) + i2 = a%irp(ir+1) + nc = i2-i1 + ip = psb_ibsrch(ic,nc,a%ja(i1:i2-1)) + if (ip>0) then + a%val(i1+ip-1) = a%val(i1+ip-1) + val(i) + else + info = i + return + end if + else + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Discarding row that does not belong to us.' + end if + end do + + case default + info = -3 + if (debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),& + & ': Duplicate handling: ',dupl + end select + + end if + + end subroutine psb_d_cyy_srch_upd + end subroutine psb_d_cyy_csput + subroutine psb_d_cyy_reinit(a,clear) - use psb_base_mod + use psb_error_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_reinit implicit none @@ -363,27 +2594,84 @@ subroutine psb_d_cyy_reinit(a,clear) logical, intent(in), optional :: clear integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='reinit' logical :: clear_ logical, parameter :: debug=.false. + call psb_erractionsave(err_act) info = psb_success_ + + + if (present(clear)) then + clear_ = clear + else + clear_ = .true. + end if + + if (a%is_bld() .or. a%is_upd()) then + ! do nothing + return + else if (a%is_asb()) then + if (clear_) a%val(:) = dzero + call a%set_upd() + else + info = psb_err_invalid_mat_state_ + call psb_errpush(info,name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine psb_d_cyy_reinit subroutine psb_d_cyy_trim(a) - use psb_base_mod + use psb_realloc_mod + use psb_error_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_trim implicit none class(psb_d_cyy_sparse_mat), intent(inout) :: a integer(psb_ipk_) :: err_act, info, nz, m + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='trim' logical, parameter :: debug=.false. + call psb_erractionsave(err_act) info = psb_success_ + m = a%get_nrows() + nz = a%get_nzeros() + if (info == psb_success_) call psb_realloc(m+1,a%irp,info) + + if (info == psb_success_) call psb_realloc(nz,a%ja,info) + if (info == psb_success_) call psb_realloc(nz,a%val,info) + + if (info /= psb_success_) goto 9999 + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + if (err_act == psb_act_abort_) then + call psb_error() + return + end if + return + end subroutine psb_d_cyy_trim subroutine psb_d_cyy_print(iout,a,iv,head,ivr,ivc) - use psb_base_mod + use psb_string_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_print implicit none @@ -394,14 +2682,73 @@ subroutine psb_d_cyy_print(iout,a,iv,head,ivr,ivc) integer(psb_ipk_), intent(in), optional :: ivr(:), ivc(:) integer(psb_ipk_) :: err_act + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='d_cyy_print' logical, parameter :: debug=.false. - + character(len=*), parameter :: datatype='real' + character(len=80) :: frmtv integer(psb_ipk_) :: irs,ics,i,j, nmx, ni, nr, nc, nz + + if (present(head)) then + write(iout,'(a)') '%%MatrixMarket matrix coordinate real general' + write(iout,'(a,a)') '% ',head + write(iout,'(a)') '%' + write(iout,'(a,a)') '% COO' + endif + + nr = a%get_nrows() + nc = a%get_ncols() + nz = a%get_nzeros() + nmx = max(nr,nc,1) + ni = floor(log10(1.0*nmx)) + 1 + + if (datatype=='real') then + write(frmtv,'(a,i3.3,a,i3.3,a)') '(2(i',ni,',1x),es26.18,1x,2(i',ni,',1x))' + else + write(frmtv,'(a,i3.3,a,i3.3,a)') '(2(i',ni,',1x),2(es26.18,1x),2(i',ni,',1x))' + end if + write(iout,*) nr, nc, nz + if(present(iv)) then + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + write(iout,frmtv) iv(i),iv(a%ja(j)),a%val(j) + end do + enddo + else + if (present(ivr).and..not.present(ivc)) then + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + write(iout,frmtv) ivr(i),(a%ja(j)),a%val(j) + end do + enddo + else if (present(ivr).and.present(ivc)) then + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + write(iout,frmtv) ivr(i),ivc(a%ja(j)),a%val(j) + end do + enddo + else if (.not.present(ivr).and.present(ivc)) then + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + write(iout,frmtv) (i),ivc(a%ja(j)),a%val(j) + end do + enddo + else if (.not.present(ivr).and..not.present(ivc)) then + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + write(iout,frmtv) (i),(a%ja(j)),a%val(j) + end do + enddo + endif + endif + end subroutine psb_d_cyy_print + subroutine psb_d_cp_cyy_from_coo(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_realloc_mod + use psb_d_base_mat_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cp_cyy_from_coo implicit none @@ -413,16 +2760,23 @@ subroutine psb_d_cp_cyy_from_coo(a,b,info) integer(psb_ipk_), allocatable :: itemp(:) !locals logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, i,j,irw, idl,err_act, nc + integer(psb_ipk_) :: nza, nr, i,j,irw, err_act, nc integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name - info = psb_success_ + info = psb_success_ + ! This is to have fix_coo called behind the scenes + call tmp%cp_from_coo(b,info) + if (info == psb_success_) call a%mv_from_coo(tmp,info) + end subroutine psb_d_cp_cyy_from_coo + + subroutine psb_d_cp_cyy_to_coo(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_d_base_mat_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cp_cyy_to_coo implicit none @@ -433,16 +2787,38 @@ subroutine psb_d_cp_cyy_to_coo(a,b,info) integer(psb_ipk_), allocatable :: itemp(:) !locals logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, nc,i,j,irw, idl,err_act + integer(psb_ipk_) :: nza, nr, nc,i,j,irw, err_act integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name info = psb_success_ + + nr = a%get_nrows() + nc = a%get_ncols() + nza = a%get_nzeros() + + call b%allocate(nr,nc,nza) + call b%psb_d_base_sparse_mat%cp_from(a%psb_d_base_sparse_mat) + + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + b%ia(j) = i + b%ja(j) = a%ja(j) + b%val(j) = a%val(j) + end do + end do + call b%set_nzeros(a%get_nzeros()) + call b%fix(info) + + end subroutine psb_d_cp_cyy_to_coo + subroutine psb_d_mv_cyy_to_coo(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_realloc_mod + use psb_d_base_mat_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_mv_cyy_to_coo implicit none @@ -453,16 +2829,41 @@ subroutine psb_d_mv_cyy_to_coo(a,b,info) integer(psb_ipk_), allocatable :: itemp(:) !locals logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, nc,i,j,irw, idl,err_act + integer(psb_ipk_) :: nza, nr, nc,i,j,irw, err_act integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name info = psb_success_ + + nr = a%get_nrows() + nc = a%get_ncols() + nza = a%get_nzeros() + + call b%psb_d_base_sparse_mat%mv_from(a%psb_d_base_sparse_mat) + call b%set_nzeros(a%get_nzeros()) + call move_alloc(a%ja,b%ja) + call move_alloc(a%val,b%val) + call psb_realloc(nza,b%ia,info) + if (info /= psb_success_) return + do i=1, nr + do j=a%irp(i),a%irp(i+1)-1 + b%ia(j) = i + end do + end do + call a%free() + call b%fix(info) + + end subroutine psb_d_mv_cyy_to_coo + + subroutine psb_d_mv_cyy_from_coo(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_realloc_mod + use psb_error_mod + use psb_d_base_mat_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_mv_cyy_from_coo implicit none @@ -473,16 +2874,88 @@ subroutine psb_d_mv_cyy_from_coo(a,b,info) integer(psb_ipk_), allocatable :: itemp(:) !locals logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, i,j,irw, idl,err_act, nc + integer(psb_ipk_) :: nza, nr, i,j,irw, err_act, nc integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit - character(len=20) :: name + character(len=20) :: name='mv_from_coo' info = psb_success_ + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + + call b%fix(info) + if (info /= psb_success_) return + + nr = b%get_nrows() + nc = b%get_ncols() + nza = b%get_nzeros() + + call a%psb_d_base_sparse_mat%mv_from(b%psb_d_base_sparse_mat) + + ! Dirty trick: call move_alloc to have the new data allocated just once. + call move_alloc(b%ia,itemp) + call move_alloc(b%ja,a%ja) + call move_alloc(b%val,a%val) + call psb_realloc(max(nr+1,nc+1),a%irp,info) + call b%free() + + if (nza <= 0) then + a%irp(:) = 1 + else + a%irp(1) = 1 + if (nr < itemp(nza)) then + write(debug_unit,*) trim(name),': RWSHR=.false. : ',& + &nr,itemp(nza),' Expect trouble!' + info = 12 + end if + + j = 1 + i = 1 + irw = itemp(j) + + outer: do + inner: do + if (i >= irw) exit inner + if (i>nr) then + write(debug_unit,*) trim(name),& + & 'Strange situation: i>nr ',i,nr,j,nza,irw + exit outer + end if + a%irp(i+1) = a%irp(i) + i = i + 1 + end do inner + j = j + 1 + if (j > nza) exit + if (itemp(j) /= irw) then + a%irp(i+1) = j + irw = itemp(j) + i = i + 1 + endif + if (i>nr) exit + enddo outer + ! + ! Cleanup empty rows at the end + ! + if (j /= (nza+1)) then + write(debug_unit,*) trim(name),': Problem from loop :',j,nza + info = 13 + endif + do + if (i>nr) exit + a%irp(i+1) = j + i = i + 1 + end do + + endif + + end subroutine psb_d_mv_cyy_from_coo + subroutine psb_d_mv_cyy_to_fmt(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_d_base_mat_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_mv_cyy_to_fmt implicit none @@ -493,17 +2966,36 @@ subroutine psb_d_mv_cyy_to_fmt(a,b,info) !locals type(psb_d_coo_sparse_mat) :: tmp logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, i,j,irw, idl,err_act, nc + integer(psb_ipk_) :: nza, nr, i,j,irw, err_act, nc integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name info = psb_success_ + select type (b) + type is (psb_d_coo_sparse_mat) + call a%mv_to_coo(b,info) + ! Need to fix trivial copies! + type is (psb_d_cyy_sparse_mat) + call b%psb_d_base_sparse_mat%mv_from(a%psb_d_base_sparse_mat) + call move_alloc(a%irp, b%irp) + call move_alloc(a%ja, b%ja) + call move_alloc(a%val, b%val) + call a%free() + + class default + call a%mv_to_coo(tmp,info) + if (info == psb_success_) call b%mv_from_coo(tmp,info) + end select + end subroutine psb_d_mv_cyy_to_fmt + subroutine psb_d_cp_cyy_to_fmt(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_d_base_mat_mod + use psb_realloc_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cp_cyy_to_fmt implicit none @@ -514,16 +3006,35 @@ subroutine psb_d_cp_cyy_to_fmt(a,b,info) !locals type(psb_d_coo_sparse_mat) :: tmp logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, i,j,irw, idl,err_act, nc + integer(psb_ipk_) :: nza, nr, i,j,irw, err_act, nc integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name info = psb_success_ + + + select type (b) + type is (psb_d_coo_sparse_mat) + call a%cp_to_coo(b,info) + + type is (psb_d_cyy_sparse_mat) + call b%psb_d_base_sparse_mat%cp_from(a%psb_d_base_sparse_mat) + call psb_safe_cpy( a%irp, b%irp , info) + call psb_safe_cpy( a%ja , b%ja , info) + call psb_safe_cpy( a%val, b%val , info) + + class default + call a%cp_to_coo(tmp,info) + if (info == psb_success_) call b%mv_from_coo(tmp,info) + end select + end subroutine psb_d_cp_cyy_to_fmt + subroutine psb_d_mv_cyy_from_fmt(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_d_base_mat_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_mv_cyy_from_fmt implicit none @@ -534,17 +3045,37 @@ subroutine psb_d_mv_cyy_from_fmt(a,b,info) !locals type(psb_d_coo_sparse_mat) :: tmp logical :: rwshr_ - integer(psb_ipk_) :: nza, nr, i,j,irw, idl,err_act, nc + integer(psb_ipk_) :: nza, nr, i,j,irw, err_act, nc integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name info = psb_success_ + select type (b) + type is (psb_d_coo_sparse_mat) + call a%mv_from_coo(b,info) + + type is (psb_d_cyy_sparse_mat) + call a%psb_d_base_sparse_mat%mv_from(b%psb_d_base_sparse_mat) + call move_alloc(b%irp, a%irp) + call move_alloc(b%ja, a%ja) + call move_alloc(b%val, a%val) + call b%free() + + class default + call b%mv_to_coo(tmp,info) + if (info == psb_success_) call a%mv_from_coo(tmp,info) + end select + end subroutine psb_d_mv_cyy_from_fmt + + subroutine psb_d_cp_cyy_from_fmt(a,b,info) - use psb_base_mod + use psb_const_mod + use psb_d_base_mat_mod + use psb_realloc_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cp_cyy_from_fmt implicit none @@ -555,41 +3086,106 @@ subroutine psb_d_cp_cyy_from_fmt(a,b,info) !locals type(psb_d_coo_sparse_mat) :: tmp logical :: rwshr_ - integer(psb_ipk_) :: nz, nr, i,j,irw, idl,err_act, nc + integer(psb_ipk_) :: nz, nr, i,j,irw, err_act, nc integer(psb_ipk_), Parameter :: maxtry=8 integer(psb_ipk_) :: debug_level, debug_unit character(len=20) :: name info = psb_success_ + + select type (b) + type is (psb_d_coo_sparse_mat) + call a%cp_from_coo(b,info) + + type is (psb_d_cyy_sparse_mat) + call a%psb_d_base_sparse_mat%cp_from(b%psb_d_base_sparse_mat) + call psb_safe_cpy( b%irp, a%irp , info) + call psb_safe_cpy( b%ja , a%ja , info) + call psb_safe_cpy( b%val, a%val , info) + + class default + call b%cp_to_coo(tmp,info) + if (info == psb_success_) call a%mv_from_coo(tmp,info) + end select end subroutine psb_d_cp_cyy_from_fmt + subroutine psb_d_cyy_cp_from(a,b) - use psb_base_mod + use psb_error_mod + use psb_realloc_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_cp_from implicit none class(psb_d_cyy_sparse_mat), intent(inout) :: a type(psb_d_cyy_sparse_mat), intent(in) :: b + + integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='cp_from' logical, parameter :: debug=.false. + call psb_erractionsave(err_act) + info = psb_success_ + call a%allocate(b%get_nrows(),b%get_ncols(),b%get_nzeros()) + call a%psb_d_base_sparse_mat%cp_from(b%psb_d_base_sparse_mat) + call psb_safe_cpy( b%irp, a%irp , info) + call psb_safe_cpy( b%ja , a%ja , info) + call psb_safe_cpy( b%val, a%val , info) + + if (info /= psb_success_) goto 9999 + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + end subroutine psb_d_cyy_cp_from subroutine psb_d_cyy_mv_from(a,b) - use psb_base_mod + use psb_error_mod use psb_d_cyy_mat_mod, psb_protect_name => psb_d_cyy_mv_from implicit none class(psb_d_cyy_sparse_mat), intent(inout) :: a type(psb_d_cyy_sparse_mat), intent(inout) :: b + + integer(psb_ipk_) :: err_act, info + integer(psb_ipk_) :: ierr(5) character(len=20) :: name='mv_from' logical, parameter :: debug=.false. call psb_erractionsave(err_act) info = psb_success_ + call a%psb_d_base_sparse_mat%mv_from(b%psb_d_base_sparse_mat) + call move_alloc(b%irp, a%irp) + call move_alloc(b%ja, a%ja) + call move_alloc(b%val, a%val) + call b%free() + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + end subroutine psb_d_cyy_mv_from + diff --git a/test/serial/psb_d_cyy_mat_mod.f90 b/test/serial/psb_d_cyy_mat_mod.f90 index ef7231c45..5f7dde619 100644 --- a/test/serial/psb_d_cyy_mat_mod.f90 +++ b/test/serial/psb_d_cyy_mat_mod.f90 @@ -1,10 +1,58 @@ +!!$ +!!$ Parallel Sparse BLAS version 3.0 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010, 2012 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ 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. +!!$ +!!$ +! +! package: psb_d_cyy_mat_mod +! +! This module contains the definition of the psb_d_cyy_sparse_mat type +! which implements an actual storage format. It is a fictitious example, +! identical to CSR except in name, to show how to extend the library. +! module psb_d_cyy_mat_mod use psb_d_base_mat_mod + !> \namespace psb_base_mod \class psb_d_cyy_sparse_mat + !! \extends psb_d_base_mat_mod::psb_d_base_sparse_mat + !! + !! psb_d_cyy_sparse_mat type and the related methods. + !! This is a very common storage type, and is the default for assembled + !! matrices in our library type, extends(psb_d_base_sparse_mat) :: psb_d_cyy_sparse_mat - integer(psb_ipk_), allocatable :: irp(:), ja(:) + !> Pointers to beginning of rows in JA and VAL. + integer(psb_ipk_), allocatable :: irp(:) + !> Column indices. + integer(psb_ipk_), allocatable :: ja(:) + !> Coefficient values. real(psb_dpk_), allocatable :: val(:) contains @@ -18,6 +66,7 @@ module psb_d_cyy_mat_mod procedure, pass(a) :: d_inner_cssv => psb_d_cyy_cssv procedure, pass(a) :: d_scals => psb_d_cyy_scals procedure, pass(a) :: d_scal => psb_d_cyy_scal + procedure, pass(a) :: maxval => psb_d_cyy_maxval procedure, pass(a) :: csnmi => psb_d_cyy_csnmi procedure, pass(a) :: csnm1 => psb_d_cyy_csnm1 procedure, pass(a) :: rowsum => psb_d_cyy_rowsum @@ -54,50 +103,65 @@ module psb_d_cyy_mat_mod private :: d_cyy_get_nzeros, d_cyy_free, d_cyy_get_fmt, & & d_cyy_get_size, d_cyy_sizeof, d_cyy_get_nz_row + !> \memberof psb_d_cyy_sparse_mat + !| \see psb_base_mat_mod::psb_base_reallocate_nz interface subroutine psb_d_cyy_reallocate_nz(nz,a) - import :: psb_ipk_, psb_d_cyy_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat integer(psb_ipk_), intent(in) :: nz class(psb_d_cyy_sparse_mat), intent(inout) :: a end subroutine psb_d_cyy_reallocate_nz end interface + !> \memberof psb_d_cyy_sparse_mat + !| \see psb_base_mat_mod::psb_base_reinit interface subroutine psb_d_cyy_reinit(a,clear) - import :: psb_ipk_, psb_d_cyy_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat class(psb_d_cyy_sparse_mat), intent(inout) :: a logical, intent(in), optional :: clear end subroutine psb_d_cyy_reinit end interface + !> \memberof psb_d_cyy_sparse_mat + !| \see psb_base_mat_mod::psb_base_trim interface subroutine psb_d_cyy_trim(a) - import :: psb_ipk_, psb_d_cyy_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat class(psb_d_cyy_sparse_mat), intent(inout) :: a end subroutine psb_d_cyy_trim end interface + + !> \memberof psb_d_cyy_sparse_mat + !| \see psb_base_mat_mod::psb_base_mold interface subroutine psb_d_cyy_mold(a,b,info) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_base_sparse_mat, psb_long_int_k_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_base_sparse_mat, psb_long_int_k_ class(psb_d_cyy_sparse_mat), intent(in) :: a class(psb_d_base_sparse_mat), intent(out), allocatable :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_cyy_mold end interface + + !> \memberof psb_d_cyy_sparse_mat + !| \see psb_base_mat_mod::psb_base_allocate_mnnz interface subroutine psb_d_cyy_allocate_mnnz(m,n,a,nz) - import :: psb_ipk_, psb_d_cyy_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat integer(psb_ipk_), intent(in) :: m,n class(psb_d_cyy_sparse_mat), intent(inout) :: a integer(psb_ipk_), intent(in), optional :: nz end subroutine psb_d_cyy_allocate_mnnz end interface + + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_print interface subroutine psb_d_cyy_print(iout,a,iv,head,ivr,ivc) - import :: psb_ipk_, psb_d_cyy_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat integer(psb_ipk_), intent(in) :: iout class(psb_d_cyy_sparse_mat), intent(in) :: a integer(psb_ipk_), intent(in), optional :: iv(:) @@ -106,98 +170,120 @@ module psb_d_cyy_mat_mod end subroutine psb_d_cyy_print end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cp_to_coo interface subroutine psb_d_cp_cyy_to_coo(a,b,info) - import :: psb_ipk_, psb_d_coo_sparse_mat, psb_d_cyy_sparse_mat + import :: psb_ipk_, psb_d_coo_sparse_mat, psb_d_cyy_sparse_mat class(psb_d_cyy_sparse_mat), intent(in) :: a class(psb_d_coo_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_cp_cyy_to_coo end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cp_from_coo interface subroutine psb_d_cp_cyy_from_coo(a,b,info) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_coo_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_coo_sparse_mat class(psb_d_cyy_sparse_mat), intent(inout) :: a class(psb_d_coo_sparse_mat), intent(in) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_cp_cyy_from_coo end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cp_to_fmt interface subroutine psb_d_cp_cyy_to_fmt(a,b,info) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_base_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_base_sparse_mat class(psb_d_cyy_sparse_mat), intent(in) :: a class(psb_d_base_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_cp_cyy_to_fmt end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cp_from_fmt interface subroutine psb_d_cp_cyy_from_fmt(a,b,info) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_base_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_base_sparse_mat class(psb_d_cyy_sparse_mat), intent(inout) :: a class(psb_d_base_sparse_mat), intent(in) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_cp_cyy_from_fmt end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_mv_to_coo interface subroutine psb_d_mv_cyy_to_coo(a,b,info) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_coo_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_coo_sparse_mat class(psb_d_cyy_sparse_mat), intent(inout) :: a class(psb_d_coo_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_mv_cyy_to_coo end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_mv_from_coo interface subroutine psb_d_mv_cyy_from_coo(a,b,info) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_coo_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_coo_sparse_mat class(psb_d_cyy_sparse_mat), intent(inout) :: a class(psb_d_coo_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_mv_cyy_from_coo end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_mv_to_fmt interface subroutine psb_d_mv_cyy_to_fmt(a,b,info) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_base_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_base_sparse_mat class(psb_d_cyy_sparse_mat), intent(inout) :: a class(psb_d_base_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_mv_cyy_to_fmt end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_mv_from_fmt interface subroutine psb_d_mv_cyy_from_fmt(a,b,info) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_base_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_d_base_sparse_mat class(psb_d_cyy_sparse_mat), intent(inout) :: a class(psb_d_base_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_mv_cyy_from_fmt end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cp_from interface subroutine psb_d_cyy_cp_from(a,b) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(inout) :: a type(psb_d_cyy_sparse_mat), intent(in) :: b end subroutine psb_d_cyy_cp_from end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_mv_from interface subroutine psb_d_cyy_mv_from(a,b) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(inout) :: a type(psb_d_cyy_sparse_mat), intent(inout) :: b end subroutine psb_d_cyy_mv_from end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csput interface subroutine psb_d_cyy_csput(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(inout) :: a real(psb_dpk_), intent(in) :: val(:) integer(psb_ipk_), intent(in) :: nz,ia(:), ja(:),& @@ -207,10 +293,12 @@ module psb_d_cyy_mat_mod end subroutine psb_d_cyy_csput end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_base_mat_mod::psb_base_csgetptn interface subroutine psb_d_cyy_csgetptn(imin,imax,a,nz,ia,ja,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a integer(psb_ipk_), intent(in) :: imin,imax integer(psb_ipk_), intent(out) :: nz @@ -223,10 +311,12 @@ module psb_d_cyy_mat_mod end subroutine psb_d_cyy_csgetptn end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csgetrow interface subroutine psb_d_cyy_csgetrow(imin,imax,a,nz,ia,ja,val,info,& & jmin,jmax,iren,append,nzin,rscale,cscale) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a integer(psb_ipk_), intent(in) :: imin,imax integer(psb_ipk_), intent(out) :: nz @@ -240,10 +330,12 @@ module psb_d_cyy_mat_mod end subroutine psb_d_cyy_csgetrow end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csgetblk interface subroutine psb_d_cyy_csgetblk(imin,imax,a,b,info,& & jmin,jmax,iren,append,rscale,cscale) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_, psb_d_coo_sparse_mat + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_, psb_d_coo_sparse_mat class(psb_d_cyy_sparse_mat), intent(in) :: a class(psb_d_coo_sparse_mat), intent(inout) :: b integer(psb_ipk_), intent(in) :: imin,imax @@ -255,17 +347,23 @@ module psb_d_cyy_mat_mod end subroutine psb_d_cyy_csgetblk end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cssv interface subroutine psb_d_cyy_cssv(alpha,a,x,beta,y,info,trans) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(in) :: alpha, beta, x(:) real(psb_dpk_), intent(inout) :: y(:) integer(psb_ipk_), intent(out) :: info character, optional, intent(in) :: trans end subroutine psb_d_cyy_cssv + end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_cssm + interface subroutine psb_d_cyy_cssm(alpha,a,x,beta,y,info,trans) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) real(psb_dpk_), intent(inout) :: y(:,:) @@ -274,17 +372,24 @@ module psb_d_cyy_mat_mod end subroutine psb_d_cyy_cssm end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csmv interface subroutine psb_d_cyy_csmv(alpha,a,x,beta,y,info,trans) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(in) :: alpha, beta, x(:) real(psb_dpk_), intent(inout) :: y(:) integer(psb_ipk_), intent(out) :: info character, optional, intent(in) :: trans end subroutine psb_d_cyy_csmv + end interface + + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csmm + interface subroutine psb_d_cyy_csmm(alpha,a,x,beta,y,info,trans) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) real(psb_dpk_), intent(inout) :: y(:,:) @@ -294,75 +399,104 @@ module psb_d_cyy_mat_mod end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_maxval + interface + function psb_d_cyy_maxval(a) result(res) + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + class(psb_d_cyy_sparse_mat), intent(in) :: a + real(psb_dpk_) :: res + end function psb_d_cyy_maxval + end interface + + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csnmi interface function psb_d_cyy_csnmi(a) result(res) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_) :: res end function psb_d_cyy_csnmi end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_csnm1 interface function psb_d_cyy_csnm1(a) result(res) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_) :: res end function psb_d_cyy_csnm1 end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_rowsum interface subroutine psb_d_cyy_rowsum(d,a) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) end subroutine psb_d_cyy_rowsum end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_arwsum interface subroutine psb_d_cyy_arwsum(d,a) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) end subroutine psb_d_cyy_arwsum end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_colsum interface subroutine psb_d_cyy_colsum(d,a) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) end subroutine psb_d_cyy_colsum end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_aclsum interface subroutine psb_d_cyy_aclsum(d,a) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) end subroutine psb_d_cyy_aclsum end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_get_diag interface subroutine psb_d_cyy_get_diag(a,d,info) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) integer(psb_ipk_), intent(out) :: info end subroutine psb_d_cyy_get_diag end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_scal interface - subroutine psb_d_cyy_scal(d,a,info) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + subroutine psb_d_cyy_scal(d,a,info,side) + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(inout) :: a real(psb_dpk_), intent(in) :: d(:) integer(psb_ipk_), intent(out) :: info + character, intent(in), optional :: side end subroutine psb_d_cyy_scal end interface + !> \memberof psb_d_cyy_sparse_mat + !! \see psb_d_base_mat_mod::psb_d_base_scals interface subroutine psb_d_cyy_scals(d,a,info) - import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ + import :: psb_ipk_, psb_d_cyy_sparse_mat, psb_dpk_ class(psb_d_cyy_sparse_mat), intent(inout) :: a real(psb_dpk_), intent(in) :: d integer(psb_ipk_), intent(out) :: info @@ -373,11 +507,28 @@ module psb_d_cyy_mat_mod contains + ! == =================================== + ! + ! + ! + ! Getters + ! + ! + ! + ! + ! + ! == =================================== + + function d_cyy_sizeof(a) result(res) implicit none class(psb_d_cyy_sparse_mat), intent(in) :: a integer(psb_long_int_k_) :: res - res = 0 + res = 8 + res = res + psb_sizeof_dp * size(a%val) + res = res + psb_sizeof_int * size(a%irp) + res = res + psb_sizeof_int * size(a%ja) + end function d_cyy_sizeof function d_cyy_get_fmt() result(res) @@ -390,14 +541,27 @@ contains implicit none class(psb_d_cyy_sparse_mat), intent(in) :: a integer(psb_ipk_) :: res - res = 0 + res = a%irp(a%get_nrows()+1)-1 end function d_cyy_get_nzeros function d_cyy_get_size(a) result(res) implicit none class(psb_d_cyy_sparse_mat), intent(in) :: a integer(psb_ipk_) :: res - res = 0 + + res = -1 + + if (allocated(a%ja)) then + res = size(a%ja) + end if + if (allocated(a%val)) then + if (res >= 0) then + res = min(res,size(a%val)) + else + res = size(a%val) + end if + end if + end function d_cyy_get_size @@ -409,13 +573,44 @@ contains class(psb_d_cyy_sparse_mat), intent(in) :: a integer(psb_ipk_), intent(in) :: idx integer(psb_ipk_) :: res + res = 0 + + if ((1<=idx).and.(idx<=a%get_nrows())) then + res = a%irp(idx+1)-a%irp(idx) + end if + end function d_cyy_get_nz_row + + + ! == =================================== + ! + ! + ! + ! Data management + ! + ! + ! + ! + ! + ! == =================================== + subroutine d_cyy_free(a) implicit none + class(psb_d_cyy_sparse_mat), intent(inout) :: a + + if (allocated(a%irp)) deallocate(a%irp) + if (allocated(a%ja)) deallocate(a%ja) + if (allocated(a%val)) deallocate(a%val) + call a%set_null() + call a%set_nrows(izero) + call a%set_ncols(izero) + return + end subroutine d_cyy_free + end module psb_d_cyy_mat_mod