! Parallel Sparse BLAS GPU plugin ! (C) Copyright 2013 ! ! Salvatore Filippone ! Alessandro Fanfarillo ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: ! 1. Redistributions of source code must retain the above copyright ! notice, this list of conditions and the following disclaimer. ! 2. Redistributions in binary form must reproduce the above copyright ! notice, this list of conditions, and the following disclaimer in the ! documentation and/or other materials provided with the distribution. ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR ! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS ! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR ! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF ! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS ! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. ! module psb_d_mfree_mat_mod use psb_d_base_mat_mod type, extends(psb_d_base_sparse_mat) :: psb_d_mfree_sparse_mat ! ! MFREE: try out a matrix-free approach ! integer(psb_ipk_) :: nnz real(psb_dpk_), allocatable :: val(:,:) contains procedure, pass(a) :: get_size => d_mfree_get_size procedure, pass(a) :: get_nzeros => d_mfree_get_nzeros procedure, nopass :: get_fmt => d_mfree_get_fmt procedure, pass(a) :: sizeof => d_mfree_sizeof procedure, pass(a) :: csmv => psb_d_mfree_csmv procedure, pass(a) :: csmm => psb_d_mfree_csmm procedure, pass(a) :: csnmi => psb_d_mfree_csnmi procedure, pass(a) :: reallocate_nz => psb_d_mfree_reallocate_nz procedure, pass(a) :: allocate_mnnz => psb_d_mfree_allocate_mnnz procedure, pass(a) :: cp_to_coo => psb_d_cp_free_to_coo procedure, pass(a) :: cp_from_coo => psb_d_cp_free_from_coo procedure, pass(a) :: mv_to_coo => psb_d_mv_free_to_coo procedure, pass(a) :: mv_from_coo => psb_d_mv_free_from_coo procedure, pass(a) :: get_diag => psb_d_mfree_get_diag procedure, pass(a) :: csgetrow => psb_d_mfree_csgetrow procedure, pass(a) :: get_nz_row => d_mfree_get_nz_row procedure, pass(a) :: trim => psb_d_mfree_trim procedure, pass(a) :: free => d_mfree_free procedure, pass(a) :: mold => psb_d_mfree_mold end type psb_d_mfree_sparse_mat private :: d_mfree_get_nzeros, d_mfree_free, d_mfree_get_fmt, & & d_mfree_get_size, d_mfree_sizeof, d_mfree_get_nz_row ! ! !> Function reallocate_nz !! \memberof psb_d_mfree_sparse_mat !! \brief One--parameters version of (re)allocate !! !! \param nz number of nonzeros to allocate for !! i.e. makes sure that the internal storage !! allows for NZ coefficients and their indices. ! interface subroutine psb_d_mfree_reallocate_nz(nz,a) import :: psb_d_mfree_sparse_mat, psb_ipk_ integer(psb_ipk_), intent(in) :: nz class(psb_d_mfree_sparse_mat), intent(inout) :: a end subroutine psb_d_mfree_reallocate_nz end interface !> Function trim !! \memberof psb_d_mfree_sparse_mat !! \brief Memory trim !! Make sure the memory allocation of the sparse matrix is as tight as !! possible given the actual number of nonzeros it contains. ! interface subroutine psb_d_mfree_trim(a) import :: psb_d_mfree_sparse_mat class(psb_d_mfree_sparse_mat), intent(inout) :: a end subroutine psb_d_mfree_trim end interface ! !> Function mold: !! \memberof psb_d_mfree_sparse_mat !! \brief Allocate a class(psb_d_mfree_sparse_mat) with the !! same dynamic type as the input. !! This is equivalent to allocate( mold= ) and is provided !! for those compilers not yet supporting mold. !! \param b The output variable !! \param info return code ! interface subroutine psb_d_mfree_mold(a,b,info) import :: psb_d_mfree_sparse_mat, psb_d_base_sparse_mat, psb_epk_, psb_ipk_ class(psb_d_mfree_sparse_mat), intent(in) :: a class(psb_d_base_sparse_mat), intent(inout), allocatable :: b integer(psb_ipk_), intent(out) :: info end subroutine psb_d_mfree_mold end interface ! ! !> Function allocate_mnnz !! \memberof psb_d_mfree_sparse_mat !! \brief Three-parameters version of allocate !! !! \param m number of rows !! \param n number of cols !! \param nz [estimated internally] number of nonzeros to allocate for ! interface subroutine psb_d_mfree_allocate_mnnz(m,n,a,nz) import :: psb_d_mfree_sparse_mat, psb_ipk_ integer(psb_ipk_), intent(in) :: m,n class(psb_d_mfree_sparse_mat), intent(inout) :: a integer(psb_ipk_), intent(in), optional :: nz end subroutine psb_d_mfree_allocate_mnnz end interface ! !> Function cp_to_coo: !! \memberof psb_d_mfree_sparse_mat !! \brief Copy and convert to psb_d_coo_sparse_mat !! Invoked from the source object. !! \param b The output variable !! \param info return code ! interface subroutine psb_d_cp_free_to_coo(a,b,info) import :: psb_d_coo_sparse_mat, psb_d_mfree_sparse_mat, psb_ipk_ class(psb_d_mfree_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_free_to_coo end interface ! !> Function cp_from_coo: !! \memberof psb_d_mfree_sparse_mat !! \brief Copy and convert from psb_d_coo_sparse_mat !! Invoked from the target object. !! \param b The input variable !! \param info return code ! interface subroutine psb_d_cp_free_from_coo(a,b,info) import :: psb_d_mfree_sparse_mat, psb_d_coo_sparse_mat, psb_ipk_ class(psb_d_mfree_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_free_from_coo end interface ! !> Function mv_to_coo: !! \memberof psb_d_mfree_sparse_mat !! \brief Convert to psb_d_coo_sparse_mat, freeing the source. !! Invoked from the source object. !! \param b The output variable !! \param info return code ! interface subroutine psb_d_mv_free_to_coo(a,b,info) import :: psb_d_mfree_sparse_mat, psb_d_coo_sparse_mat, psb_ipk_ class(psb_d_mfree_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_free_to_coo end interface ! !> Function mv_from_coo: !! \memberof psb_d_mfree_sparse_mat !! \brief Convert from psb_d_coo_sparse_mat, freeing the source. !! Invoked from the target object. !! \param b The input variable !! \param info return code ! interface subroutine psb_d_mv_free_from_coo(a,b,info) import :: psb_d_mfree_sparse_mat, psb_d_coo_sparse_mat, psb_ipk_ class(psb_d_mfree_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_free_from_coo end interface ! ! !> Function csgetrow: !! \memberof psb_d_mfree_sparse_mat !! \brief Get a (subset of) row(s) !! !! getrow is the basic method by which the other (getblk, clip) can !! be implemented. !! !! Returns the set !! NZ, IA(1:nz), JA(1:nz), VAL(1:NZ) !! each identifying the position of a nonzero in A !! i.e. !! VAL(1:NZ) = A(IA(1:NZ),JA(1:NZ)) !! with IMIN<=IA(:)<=IMAX !! with JMIN<=JA(:)<=JMAX !! IA,JA are reallocated as necessary. !! !! \param imin the minimum row index we are interested in !! \param imax the minimum row index we are interested in !! \param nz the number of output coefficients !! \param ia(:) the output row indices !! \param ja(:) the output col indices !! \param val(:) the output coefficients !! \param info return code !! \param jmin [1] minimum col index !! \param jmax [a\%get_ncols()] maximum col index !! \param iren(:) [none] an array to return renumbered indices (iren(ia(:)),iren(ja(:)) !! \param rscale [false] map [min(ia(:)):max(ia(:))] onto [1:max(ia(:))-min(ia(:))+1] !! \param cscale [false] map [min(ja(:)):max(ja(:))] onto [1:max(ja(:))-min(ja(:))+1] !! ( iren cannot be specified with rscale/cscale) !! \param append [false] append to ia,ja !! \param nzin [none] if append, then first new entry should go in entry nzin+1 !! ! interface subroutine psb_d_mfree_csgetrow(imin,imax,a,nz,ia,ja,val,info,& & jmin,jmax,iren,append,nzin,rscale,cscale,chksz) import :: psb_d_mfree_sparse_mat, psb_dpk_, psb_ipk_ class(psb_d_mfree_sparse_mat), intent(in) :: a integer(psb_ipk_), intent(in) :: imin,imax integer(psb_ipk_), intent(out) :: nz integer(psb_ipk_), allocatable, intent(inout) :: ia(:), ja(:) real(psb_dpk_), allocatable, intent(inout) :: val(:) integer(psb_ipk_),intent(out) :: info logical, intent(in), optional :: append integer(psb_ipk_), intent(in), optional :: iren(:) integer(psb_ipk_), intent(in), optional :: jmin,jmax, nzin logical, intent(in), optional :: rscale,cscale,chksz end subroutine psb_d_mfree_csgetrow end interface !> Function csmv: !! \memberof psb_d_mfree_sparse_mat !! \brief Product by a dense rank 1 array. !! !! Compute !! Y = alpha*op(A)*X + beta*Y !! !! \param alpha Scaling factor for Ax !! \param A the input sparse matrix !! \param x(:) the input dense X !! \param beta Scaling factor for y !! \param y(:) the input/output dense Y !! \param info return code !! \param trans [N] Whether to use A (N), its transpose (T) !! or its conjugate transpose (C) !! ! interface subroutine psb_d_mfree_csmv(alpha,a,x,beta,y,info,trans) import :: psb_d_mfree_sparse_mat, psb_dpk_, psb_ipk_ class(psb_d_mfree_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_mfree_csmv end interface !> Function csmm: !! \memberof psb_d_mfree_sparse_mat !! \brief Product by a dense rank 2 array. !! !! Compute !! Y = alpha*op(A)*X + beta*Y !! !! \param alpha Scaling factor for Ax !! \param A the input sparse matrix !! \param x(:,:) the input dense X !! \param beta Scaling factor for y !! \param y(:,:) the input/output dense Y !! \param info return code !! \param trans [N] Whether to use A (N), its transpose (T) !! or its conjugate transpose (C) !! ! interface subroutine psb_d_mfree_csmm(alpha,a,x,beta,y,info,trans) import :: psb_d_mfree_sparse_mat, psb_dpk_, psb_ipk_ class(psb_d_mfree_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_mfree_csmm end interface ! ! !> Function csnmi: !! \memberof psb_d_mfree_sparse_mat !! \brief Operator infinity norm !! CSNMI = MAXVAL(SUM(ABS(A(:,:)),dim=2)) !! ! interface function psb_d_mfree_csnmi(a) result(res) import :: psb_d_mfree_sparse_mat, psb_dpk_ class(psb_d_mfree_sparse_mat), intent(in) :: a real(psb_dpk_) :: res end function psb_d_mfree_csnmi end interface ! !> Function get_diag: !! \memberof psb_d_mfree_sparse_mat !! \brief Extract the diagonal of A. !! !! D(i) = A(i:i), i=1:min(nrows,ncols) !! !! \param d(:) The output diagonal !! \param info return code. ! interface subroutine psb_d_mfree_get_diag(a,d,info) import :: psb_d_mfree_sparse_mat, psb_dpk_, psb_ipk_ class(psb_d_mfree_sparse_mat), intent(in) :: a real(psb_dpk_), intent(out) :: d(:) integer(psb_ipk_), intent(out) :: info end subroutine psb_d_mfree_get_diag end interface contains ! !> Function sizeof !! \memberof psb_d_mfree_sparse_mat !! \brief Memory occupation in bytes ! function d_mfree_sizeof(a) result(res) implicit none class(psb_d_mfree_sparse_mat), intent(in) :: a integer(psb_epk_) :: res res = psb_sizeof_dp * size(a%val) res = res + psb_sizeof_ip end function d_mfree_sizeof ! !> Function get_fmt !! \memberof psb_d_mfree_sparse_mat !! \brief return a short descriptive name (e.g. COO CSR etc.) ! function d_mfree_get_fmt() result(res) implicit none character(len=5) :: res res = 'MFREE' end function d_mfree_get_fmt ! !> Function get_nzeros !! \memberof psb_d_mfree_sparse_mat !! \brief Current number of nonzero entries ! function d_mfree_get_nzeros(a) result(res) implicit none class(psb_d_mfree_sparse_mat), intent(in) :: a integer(psb_ipk_) :: res res = a%nnz end function d_mfree_get_nzeros ! !> Function get_size !! \memberof psb_d_mfree_sparse_mat !! \brief Maximum number of nonzeros the current structure can hold ! this is fixed once you initialize the matrix, with dense storage ! you can hold up to MxN entries function d_mfree_get_size(a) result(res) implicit none class(psb_d_mfree_sparse_mat), intent(in) :: a integer(psb_ipk_) :: res res = size(a%val) end function d_mfree_get_size ! !> Function get_nz_row. !! \memberof psb_d_coo_sparse_mat !! \brief How many nonzeros in a row? !! !! \param idx The row to search. !! ! function d_mfree_get_nz_row(idx,a) result(res) implicit none class(psb_d_mfree_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 = count(a%val(idx,:) /= dzero) end if end function d_mfree_get_nz_row ! !> Function free !! \memberof psb_d_mfree_sparse_mat !! Name says all subroutine d_mfree_free(a) implicit none class(psb_d_mfree_sparse_mat), intent(inout) :: a if (allocated(a%val)) deallocate(a%val) a%nnz = 0 ! ! Mark the object as empty just in case ! call a%set_null() call a%set_nrows(izero) call a%set_ncols(izero) return end subroutine d_mfree_free end module psb_d_mfree_mat_mod