From 0ce5545cea0bef6b17d27446d877a78ccc08004e Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Thu, 27 Aug 2009 16:16:59 +0000 Subject: [PATCH] psblas3: Changelog base/modules/Makefile base/modules/psb_ip_reord_mod.f90 base/newserial/Makefile base/newserial/psbn_base_mat_mod.f03 base/newserial/psbn_coo_mat.f03 base/newserial/psbn_d_base_mat_mod.f03 base/newserial/psbn_d_coo_impl.f03 base/newserial/psbn_mat_mod.f03 base/serial/Makefile base/serial/psb_ip_reord_mod.f90 config/pac.m4 Reworked base_mat to include COO (everybody must know about COO to define A%TO_COO() and A%FROM_COO) --- Changelog | 376 +--- base/modules/Makefile | 3 +- base/{serial => modules}/psb_ip_reord_mod.f90 | 0 base/newserial/Makefile | 4 +- base/newserial/psbn_base_mat_mod.f03 | 124 ++ base/newserial/psbn_coo_mat.f03 | 1577 -------------- base/newserial/psbn_d_base_mat_mod.f03 | 561 ++++- base/newserial/psbn_d_coo_impl.f03 | 1845 +++++++++++++++++ base/newserial/psbn_mat_mod.f03 | 369 ++++ base/serial/Makefile | 4 +- config/pac.m4 | 17 +- 11 files changed, 2923 insertions(+), 1957 deletions(-) rename base/{serial => modules}/psb_ip_reord_mod.f90 (100%) delete mode 100644 base/newserial/psbn_coo_mat.f03 create mode 100644 base/newserial/psbn_d_coo_impl.f03 diff --git a/Changelog b/Changelog index 382dbe54..041653a3 100644 --- a/Changelog +++ b/Changelog @@ -1,372 +1,15 @@ -2008/09/18: Defined psb_sizeof to be integer(8). Added support - into psb_sum, psb_amx and other reductions for long int - scalars. - -2008/09/16: Implemented new scheme for index conversion. - Changed cdall with an option to suppress global checks. - -2008/09/02: Improved psi_fnd_owner performace. - -2008/09/01: Better timings in the pargen test cases. - -2008/08/28: Changed CDALL in case of VL to handle overlapped indices. - -2008/07/28: New sorting/reordering modules. - -2008/07/24: Addded HTML version of user's guide. - -2008/07/22: Fixed I/O for Harwell-Boeing and Matrix Market examples - -2008/05/27: Merged single precision branch. - -2008/04/28: Fixed trimming space in sparse matrix conversion. - Fixed performance issue in cdins. - -2008/03/25: Fix performance bug in psi_idx_ins_cnv. Changed names of - some internal components of preconditioner data structure. - -======= -2009/01/27: Renamed psb_transfer into psb_move_alloc. -2009/01/08: Require GNU Fortran 4.3 or later. -2008/11/04: Repackaged and streamlined linear maps. - -2008/10/16: Fixed internal structure of psb_inter_desc. - -2008/09/23: Fix borderline cases where one process does not own any - indices from the global space. - -2008/09/18: Defined psb_sizeof to be integer(8). Added support - into psb_sum, psb_amx and other reductions for long int - scalars. - -2008/09/16: Implemented new scheme for index conversion. - Changed cdall with an option to suppress global checks. - -2008/09/02: Improved psi_fnd_owner performace. - -2008/09/01: Better timings in the pargen test cases. - -2008/08/28: Changed CDALL in case of VL to handle overlapped indices. - -2008/07/28: New sorting/reordering modules. - -2008/07/24: Addded HTML version of user's guide. - -2008/07/22: Fixed I/O for Harwell-Boeing and Matrix Market examples - -2008/05/27: Merged single precision branch. - -2008/04/28: Fixed trimming space in sparse matrix conversion. - Fixed performance issue in cdins. - -2008/03/25: Fix performance bug in psi_idx_ins_cnv. Changed names of - some internal components of preconditioner data structure. - ->>>>>>> .merge-dx.r3587 -2008/03/27: Merged the experimental branch for implementing the AVL tree - data structure in Fortran instead of relying on C and passing - functions around to perform comparisons. There seems to be - some performance advantage, although not very large. - -2008/03/25: Merged in changes from the 2.2-maint branch re: error - messages, performance bug in psi_idx_ins_cnv. - -2008/02/26: New psb_linmap_init, psb_linmap_ins, psb_linmap_asb for a - general linear operator mapping among index spaces. - -2008/02/18: Branched off for Version 2.2 - -2008/02/08: Merged changes from intermesh branch: we now have an - inter_desc_type object. Currently we only implement the - version needed for aggregation algorithms in the algebraic - multigrid preconditioners, but we'll define more general - (linear) maps soon enough. - -2008/01/25: Various changes to variables controlling conditional - compilation on the Fortran side: removed NETLIB_BLACS, now - HAVE_METIS HAVE_ESSL_BLACS HAVE_KSENDID. - Files impacted: Make.inc.XXX, base/modules/psb_penv_mod, - util/psb_metispart_mod - - -2008/01/18: Centralized convergence checks. Still partial for RGMRES. - -2008/01/14: Merged changes for handling of transpose vs. overlap. - -2008/01/10: Changed name of GMRESR into RGMRES for consistency. - -2007/12/21: Merged in debug infrastructure, internal and html docs. - -2007/11/14: Fix INTENT(IN) on X vector in preconditioner routines. - -2007/10/15: Repackaged the sorting routines in a submodule of their - own, adding some heap management and heapsort utilities for the - benefit of the multilevel preconditioners. - -2007/09/28: Moved gelp and csrp to serial. Changed interface to - sphalo: the new one makes more sense. - Updated documentation. - -2007/09/14: Second round of serial changes: merged into trunk, fixed - JAD regeneration and srch_upd now works. - -2007/09/10: First round of serial changes: implemented serial - psb_spcnv unifying multiple functionalities. - -2007/09/04: Implemented RGMRES for complex data. - -2007/06/04: Fixed implementation of fctint and coins: assume size - arrays caused troubles on some compilers. Documentation of - set_large_threshold. - -2007/05/22: Defined psb_precinit. - -2007/05/15: Defined psb_sizeof. - -2007/05/15: Merged in various fixes coming from tests on SP5 and - HP-Itanium. - -2007/04/08: Changed the implementation of psb_sp_getrow & friends. - -2007/03/27: Merged in changes for enabling compilation on SUN. - -2007/02/22: Fixed various misalignments between real and complex. - Defined new psb_sp_clip routines. - Fixed psb_rwextd. - Changed the USE statements, minimizing size of modules and - maximizing consistency checks. - -2007/02/01: Merged serial version: we provide a minimal fake mpi to - allow compiling and running without mpi and blacs. Only - tested with gnu42 so far. - -2007/01/23: Defined new field ext_index in desc_type, and - fixed long standing inconsistency in usage of overlap for - AS preconditioners. Modified halo to accept selector for - halo_index vs. ext_index. - -2007/01/11: Migrated repository to SVN. - -2007/01/11: MLD2P4 has been moved to the new org. Now tackling the - test dirs. - -2007/01/09: First try at reorganizing directories. Subdir MLD2P4 still - to be fixed. Documentation still to be updated. - -2006/12/11: Documented options in glob_to_loc. - -2006/12/06: Fixed raw aggregation. - -2006/12/05: Taken out extra interfaces; inserted use modules with ONLY - clauses where appropriate. - -2006/11/30: Fixed a bug in raw aggregation. Note: raw aggregation - gives different results from smoothed with omega=0.0, - because in the latter we have explicitly stored zero - coefficients that would be absent in the first, thus - generating different ILU factorizations. - -2006/11/28: Merged the mods for descriptors of large index spaces to - avoid having the GLOB_TO_LOC array. Took the chance to - reorganize the descriptor build routines and define some - access functions for descriptor features and entries, so - as not to use the descriptor components directly. Tested - with AS, 2- and 3- level Post smoothers. - -2006/11/09: The allocatable version works, but under gcc42 there is a - compiler bug when using -fbounds-check. - -2006/11/08: Merged the allocatable version; hope everything works! - -2006/11/08: Branched version psblas2-2-0-maint, and defined tag - 2.0.2.6 - -2006/11/02: Done in the allocatable branch: repackaging of cdasb and - friends, taking out AVL trees where they were not - absolutely needed, and new dcsrmv routine. - -2006/11/01: Merged changes in the handling of data exchange. - -2006/10/03: Merged in the multilevel preconditioner stuff. This is - still experimental, especially the interfaces are not - stable yet. - -2006/10/03: Declared version 2.0.2.5 for reference purposes. - -2006/10/03: Fixed a bunch of minor bugs, incuding the sorting routines - imsr and imsrx. Added a default call to blacs_exit inside - psb_exit fixed a bad termination in test/pargen/ppde90.f90 - -2006/09/02: Declared version 2.0.2, after having fixed a lot of - details in the environment routines. - -2006/07/25: Defined a new psb_wtime function. Modified precset to - have a non-optional INFO dummy argument. - -2006/07/06: Fixed bug in swaptran. Added psb_krylov generic interface. - -2006/07/04: Ooops, the GetRow mod in SMMP is a performance hit. - Need to investigate further. - -2006/06/21: Bug fix in hb_read when dealing with symmetric matrices. - -2006/06/20: Rewritten symbmm and numbmm from SMMP to be intependent of - CSR storage by using GetRow. Still need to test for - performance. - -2006/06/16: Defined GetRow. This way we may close the mat objects. - Next we will rewrite SMMP to only make use of GetRow, - not to rely on CSR storage format. - -2006/05/29: Added BLACS-like routines for data communication, - broadcasts, reductions, send/receive. - -2006/05/25: Added environment management routines. - -2006/05/03: Bug fixes, plus some change in the internals for SPINS, - preparing hooks for insertion with local numbering. - -2006/04/24: Minor changes to the interface of dense tools routines, - trying to achieve a uniform look & feel. - Rewritten documentation; it is now reasonable, though not - perfect, except for the preconditioner routines. - We can now declare RC3. - -2006/04/21: A bunch of fixes related to various matrix initialization - problems that were revealed while testing on SP5. - -2006/04/18: Changed interface to spasb and csdp: better handling of - regeneration. To be tested further for sophisticated uses. - -2006/03/31: We declare RC2 now. Improved I/O routines in test/Fileread. - -2006/03/24: We have a complex version now, working (not necessarily bug free). - -2006/03/15: Started move to complex version. - -2006/03/01: Complete restructure of PREC section. - -2006/02/01: New naming scheme. - -2006/01/01: New multilevel preconditioning wih smoothed aggregation. - -2005/09 : Now enabled UMFPACK complete factorization as basis for AS. - -2005/05/04: Now enabled SuperLU complete factorization as basis for AS. - -2005/04/29: First version with decoupled 2-level. - -2005/04/06: Started work on decoupling the preconditioner aggregation - for 2-level from the main factorization. - -2005/03/30: First version of new DSC/SP allocate/insert/assembly - routines. -2005/03/17: First version of RGMRES. To be refined. - -2005/03/08: dSwapTran aligned with dSwapData. Taken out SwapOverlap. - also moved onto iSwapX. - -2005/03/07: dSwapData rewritten to achieve: 1. better performance; - 2. more flexible functionality. It is now possible to - avoid SwapOvrlap entirely, relying on just SwapData. - SwapTran is still alive, since it reads the descriptors in - "transpose" mode. Also, added work areas to preconditioner - routine, to avoid excessive allocation in the halo/overlap - exchange. - -2005/03/04: Had to put in a workaround for a gfortran bug: - tolower/toupper cannot be functions. - -2005/02/09: Explicit storage choice for the smoother. This seems - to be changing a little bit the actual preconditioner. - To be evaluated further. - -2005/02/08: Renamed F90_PSPREC to PSB_PRCAPLY and Preconditioner to - PSB_PRCBLD. Changed the way PRCAPLY decides what to do. - Still needs a PSB_PRCSET to be called before PRCBLD. - -2005/01/28: Started moving functionalities to a SERIAL F90 layer. Also - defined a new COMM layer, to enable implementing SPMM - directly in F90. - -2005/01/20: Finally taken out a direct call to the F77 DCSDP from - SPASB. - -2005/01/18: After much work, we now have 2-level Additive Schwarz - prototype implemented and working. We now start a major - code cleanup that will take some time. Mainly we want to - move a lot of the serial F77 functionality into a new F95 - serial layer, to simplify the parallel F95 code. - -2004/11/25: Following the introduction of Additive Shwarz and - variants, we have now renamed DECOMP_ and friends as - DESC_; this makes things more readable. Sooner or later - we're going to merge this into mainline, but this version - is still very much in a state of flux. - -2004/07/18: For use with gfortran we need to declare the pointer - components with NULL() initialization. This rules out - VAST and PGI. - -2004/07/15: First development version with gfortran from the current - snapshot of gcc 3.5.0. - It is now possible in PSI_dSwapData to opt for - SEND|RECEIVE|SYNC data exchange; plan is to extend to all - data exchange functions, plus making it available as an - option from the F90 level. - -2004/07/06: Merged in a lot of stuff coming mainly from the ASM - development; full merge will have to wait a little more. - Among other things: - use of psimod - new choice parms for overlap - new data exchange for swapdata, to be extended. - multicolumn CSMM. - use psrealloc - new format for marking a matrix as suitable for update. - - -2003/12/09: Changed DSALLOC and DSASB to make sure whenever a dense - matrix is allocated it is also zeroed out. - -2003/10/13: Added call to BLACS_SET in the solvers to ensure global - heterogeneous coherence in the combine operations. - -2003/09/30: Added LOC_TO_GLOB and GLOB_TO_LOC support routines. - -2003/09/30: Changed interface for smart update capabilities: choose - with optional parameters in ASB routines. - -2003/09/16: IFC 7.0 had a strange behaviour in the test programs: - sometimes the declaration of PARTS dummy argument with an - INTERFACE would not work, requiring an EXTERNAL - declaration. The proper INTERFACE works now with 7.1. - -2003/03/10: Halo data exchange in F90_PSHALO can now be applied to - integer data; create appropriate support routines. - -2002/12/05: Initial version of Fileread sample programs. - -2002/11/19: Fixes for JAD preconditioner. - -2002/11/19: Methods for patterns: create a descriptor without a - matrix. - -2001/11/16: Reviewed the interfaces: in the tools section we really - need the POINTER attribute for dense vectors, but not in - the computational routines; taking it out allows more - flexibility. - -2001/09/16: Smart update capabilities. - -2001/03/16: Renumbering routines. - -2001/01/14: Added extensions to compute multiple DOTs and AMAXs at once; - -======= Changelog. A lot less detailed than usual, at least for past history. +2009/08/25: New configure flag + --enable-serial + for serial-only compilation. +2009/06/24: Changed order of arguments in sp_scal to make it uniform with + rest of library. +2009/05/15: Changed interface to matdist. +2009/05/12: Added support for NAG Fortran compiler +2009/03/16: Release 2.3.1 + 2009/01/27: Renamed psb_transfer into psb_move_alloc. 2009/01/08: Require GNU Fortran 4.3 or later. 2008/11/04: Repackaged and streamlined linear maps. @@ -703,4 +346,3 @@ history. 2001/01/14: Added extensions to compute multiple DOTs and AMAXs at once; ->>>>>>> .merge-dx.r3592 diff --git a/base/modules/Makefile b/base/modules/Makefile index 494d594c..4d7a5d9a 100644 --- a/base/modules/Makefile +++ b/base/modules/Makefile @@ -5,7 +5,7 @@ UTIL_MODS = psb_string_mod.o psb_spmat_type.o \ psb_desc_type.o psb_sort_mod.o psb_penv_mod.o \ psb_serial_mod.o psb_tools_mod.o psb_blacs_mod.o \ psb_linmap_type_mod.o psb_comm_mod.o psb_psblas_mod.o \ - psi_serial_mod.o psi_mod.o \ + psi_serial_mod.o psi_mod.o psb_ip_reord_mod.o\ psb_check_mod.o psb_gps_mod.o psb_linmap_mod.o psb_hash_mod.o MODULES=$(BASIC_MODS) $(UTIL_MODS) @@ -26,6 +26,7 @@ lib: $(BASIC_MODS) blacsmod $(UTIL_MODS) $(OBJS) $(LIBMOD) psb_realloc_mod.o : psb_error_mod.o psb_spmat_type.o : psb_realloc_mod.o psb_error_mod.o psb_const_mod.o psb_string_mod.o psb_sort_mod.o psb_error_mod.o: psb_const_mod.o +psb_ip_reord_mod.o: psb_const_mod.o psb_penv_mod.o: psb_const_mod.o psb_error_mod.o psb_realloc_mod.o psb_blacs_mod.o psb_blacs_mod.o: psb_const_mod.o psi_serial_mod.o: psb_const_mod.o psb_realloc_mod.o diff --git a/base/serial/psb_ip_reord_mod.f90 b/base/modules/psb_ip_reord_mod.f90 similarity index 100% rename from base/serial/psb_ip_reord_mod.f90 rename to base/modules/psb_ip_reord_mod.f90 diff --git a/base/newserial/Makefile b/base/newserial/Makefile index d8138d75..5d0398c9 100644 --- a/base/newserial/Makefile +++ b/base/newserial/Makefile @@ -1,6 +1,7 @@ include ../../Make.inc -MODULES = psbn_base_mat_mod.o psbn_d_base_mat_mod.o psbn_coo_mat.o psbn_csr_mat.o psbn_mat_mod.o +MODULES = psbn_base_mat_mod.o psbn_d_base_mat_mod.o psbn_d_coo_impl.o +# psbn_csr_mat.o LIBMOD= @@ -18,6 +19,7 @@ lib: $(MODULES) $(OBJS) $(LIBMOD) psbn_mat_mod.o: psbn_base_mat_mod.o psbn_coo_mat.o psbn_csr_mat.o: psbn_d_base_mat_mod.o +psbn_d_mat_impl.o: psbn_d_base_mat_mod.o clean: /bin/rm -f $(MODULES) $(OBJS) $(MPFOBJS) *$(.mod) diff --git a/base/newserial/psbn_base_mat_mod.f03 b/base/newserial/psbn_base_mat_mod.f03 index 444a1983..159a4d95 100644 --- a/base/newserial/psbn_base_mat_mod.f03 +++ b/base/newserial/psbn_base_mat_mod.f03 @@ -2,6 +2,7 @@ module psbn_base_mat_mod use psb_const_mod + integer, parameter :: psbn_invalid_ = -1 integer, parameter :: psbn_spmat_null_=0, psbn_spmat_bld_=1 integer, parameter :: psbn_spmat_asb_=2, psbn_spmat_upd_=4 @@ -29,12 +30,29 @@ module psbn_base_mat_mod integer, private :: state, duplicate logical, private :: triangle, unitd, upper, sorted contains + procedure, pass(a) :: set_nrows + procedure, pass(a) :: set_ncols + procedure, pass(a) :: set_dupl + procedure, pass(a) :: set_state + procedure, pass(a) :: set_null + procedure, pass(a) :: set_bld + procedure, pass(a) :: set_upd + procedure, pass(a) :: set_asb + procedure, pass(a) :: set_sorted + procedure, pass(a) :: set_upper + procedure, pass(a) :: set_lower + procedure, pass(a) :: set_triangle + procedure, pass(a) :: set_unit + + procedure, pass(a) :: get_nrows procedure, pass(a) :: get_ncols procedure, pass(a) :: get_nzeros procedure, pass(a) :: get_size procedure, pass(a) :: get_state procedure, pass(a) :: get_dupl + + procedure, pass(a) :: is_null procedure, pass(a) :: is_bld procedure, pass(a) :: is_upd @@ -81,6 +99,112 @@ contains res = a%n end function get_ncols + + subroutine set_nrows(m,a) + class(psbn_base_sparse_mat), intent(inout) :: a + integer, intent(in) :: m + a%m = m + end subroutine set_nrows + + subroutine set_ncols(n,a) + class(psbn_base_sparse_mat), intent(inout) :: a + integer, intent(in) :: n + a%n = n + end subroutine set_ncols + + + subroutine set_state(n,a) + class(psbn_base_sparse_mat), intent(inout) :: a + integer, intent(in) :: n + a%state = n + end subroutine set_state + + + subroutine set_dupl(n,a) + class(psbn_base_sparse_mat), intent(inout) :: a + integer, intent(in) :: n + a%duplicate = n + end subroutine set_dupl + + subroutine set_null(a) + class(psbn_base_sparse_mat), intent(inout) :: a + + a%state = psbn_spmat_null_ + end subroutine set_null + + subroutine set_bld(a) + class(psbn_base_sparse_mat), intent(inout) :: a + + a%state = psbn_spmat_bld_ + end subroutine set_bld + + subroutine set_upd(a) + class(psbn_base_sparse_mat), intent(inout) :: a + + a%state = psbn_spmat_upd_ + end subroutine set_upd + + subroutine set_asb(a) + class(psbn_base_sparse_mat), intent(inout) :: a + + a%state = psbn_spmat_asb_ + end subroutine set_asb + + subroutine set_sorted(a,val) + class(psbn_base_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + + if (present(val)) then + a%sorted = val + else + a%sorted = .true. + end if + end subroutine set_sorted + + subroutine set_triangle(a,val) + class(psbn_base_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + + if (present(val)) then + a%triangle = val + else + a%triangle = .true. + end if + end subroutine set_triangle + + subroutine set_unit(a,val) + class(psbn_base_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + + if (present(val)) then + a%unitd = val + else + a%unitd = .true. + end if + end subroutine set_unit + + subroutine set_lower(a,val) + class(psbn_base_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + + if (present(val)) then + a%upper = .not.val + else + a%upper = .false. + end if + end subroutine set_lower + + subroutine set_upper(a,val) + class(psbn_base_sparse_mat), intent(inout) :: a + logical, intent(in), optional :: val + + if (present(val)) then + a%upper = val + else + a%upper = .true. + end if + end subroutine set_upper + function is_triangle(a) result(res) class(psbn_base_sparse_mat), intent(in) :: a logical :: res diff --git a/base/newserial/psbn_coo_mat.f03 b/base/newserial/psbn_coo_mat.f03 deleted file mode 100644 index c3da01b5..00000000 --- a/base/newserial/psbn_coo_mat.f03 +++ /dev/null @@ -1,1577 +0,0 @@ -module psbn_d_coo_sparse_mat_mod - - use psbn_d_base_mat_mod - - type, extends(psbn_d_base_sparse_mat) :: psbn_d_coo_sparse_mat - - integer :: nnz - integer, allocatable :: ia(:), ja(:) - real(psb_dpk_), allocatable :: val(:) - - contains - - procedure, pass(a) :: get_nzeros => d_coo_get_nzeros - procedure, pass(a) :: set_nzeros => d_coo_set_nzeros - procedure, pass(a) :: d_base_csmm => d_coo_csmm - procedure, pass(a) :: d_base_csmv => d_coo_csmv - procedure, pass(a) :: d_base_cssm => d_coo_cssm - procedure, pass(a) :: d_base_cssv => d_coo_cssv - procedure, pass(a) :: csins => d_coo_csins - procedure, pass(a) :: reallocate_nz => d_coo_reallocate_nz - - end type psbn_d_coo_sparse_mat - -contains - - - subroutine d_coo_reallocate_nz(nz,a) - use psb_error_mod - use psb_realloc_mod - integer, intent(in) :: nz - class(psbn_d_coo_sparse_mat), intent(inout) :: a - Integer :: err_act - character(len=20) :: name='d_coo_reallocate_nz' - logical, parameter :: debug=.false. - - call psb_erractionsave(err_act) - - call psb_realloc(nx,a%ia,a%ja,a%val,info) - - if (info /= 0) then - call psb_errpush(4000,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 d_coo_reallocate_nz - - - function d_coo_get_nzeros(a) result(res) - class(psbn_d_coo_sparse_mat), intent(in) :: a - integer :: res - res = a%nnz - end function d_coo_get_nzeros - - - subroutine d_coo_set_nzeros(nz,a) - integer, intent(in) :: nz - class(psbn_d_coo_sparse_mat), intent(inout) :: a - - a%nnz = nz - - end subroutine d_coo_set_nzeros - - - subroutine d_coo_csins(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) - use psb_error_mod - use psb_realloc_mod - class(psbn_d_coo_sparse_mat), intent(inout) :: a - real(psb_dpk_), intent(in) :: val(:) - integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax - integer, intent(out) :: info - integer, intent(in), optional :: gtl(:) - - - Integer :: err_act - character(len=20) :: name='d_coo_csins' - logical, parameter :: debug=.false. - integer :: nza, i,j,k, nzl, isza, int_err(5) - - call psb_erractionsave(err_act) - info = 0 - - if (nz <= 0) then - info = 10 - int_err(1)=1 - call psb_errpush(info,name,i_err=int_err) - goto 9999 - end if - if (size(ia) < nz) then - info = 35 - int_err(1)=2 - call psb_errpush(info,name,i_err=int_err) - goto 9999 - end if - - if (size(ja) < nz) then - info = 35 - int_err(1)=3 - call psb_errpush(info,name,i_err=int_err) - goto 9999 - end if - if (size(val) < nz) then - info = 35 - int_err(1)=4 - call psb_errpush(info,name,i_err=int_err) - goto 9999 - end if - - if (nz == 0) return - - - nza = a%get_nzeros() - isza = a%get_size() - if (a%is_bld()) then - ! Build phase. Must handle reallocations in a sensible way. - if (isza < (nza+nz)) then - call a%reallocate(max(nza+nz,int(1.5*isza))) - isza = a%get_size() - endif - - call psb_inner_ins(nz,ia,ja,val,nza,a%ia,a%ja,a%val,isza,& - & imin,imax,jmin,jmax,info,gtl) - call a%set_nzeros(nz+nza) - - else if (a%is_upd()) then - if (a%is_sorted()) then - - -!!$#ifdef FIXED_NAG_SEGV -!!$ call d_coo_srch_upd(nz,ia,ja,val,a,& -!!$ & imin,imax,jmin,jmax,info,gtl) -!!$#else - call d_coo_srch_upd(nz,ia,ja,val,& - & a%ia,a%ja,a%val,& - & a%get_dupl(),a%get_nzeros(),a%get_nrows(),& - & info,gtl) -!!$#endif - - else - info = 1121 - end if - - else - ! State is wrong. - info = 1121 - end if - if (info /= 0) 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_inner_upd(nz,ia,ja,val,nza,aspk,maxsz,& -!!$ & imin,imax,jmin,jmax,info,gtl) -!!$ implicit none -!!$ -!!$ integer, intent(in) :: nz, imin,imax,jmin,jmax,maxsz -!!$ integer, intent(in) :: ia(:),ja(:) -!!$ integer, intent(inout) :: nza -!!$ real(psb_dpk_), intent(in) :: val(:) -!!$ real(psb_dpk_), intent(inout) :: aspk(:) -!!$ integer, intent(out) :: info -!!$ integer, intent(in), optional :: gtl(:) -!!$ integer :: i,ir,ic, ng,nzl -!!$ character(len=20) :: name, ch_err -!!$ -!!$ -!!$ name='psb_inner_upd' -!!$ nzl = 0 -!!$ if (present(gtl)) then -!!$ ng = size(gtl) -!!$ if ((nza > nzl)) then -!!$ do i=1, nz -!!$ nza = nza + 1 -!!$ if (nza>maxsz) then -!!$ call psb_errpush(50,name,i_err=(/7,maxsz,5,0,nza /)) -!!$ info = -71 -!!$ return -!!$ endif -!!$ aspk(nza) = val(i) -!!$ end do -!!$ else -!!$ 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 >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then -!!$ nza = nza + 1 -!!$ if (nza>maxsz) then -!!$ info = -72 -!!$ return -!!$ endif -!!$ aspk(nza) = val(i) -!!$ end if -!!$ end if -!!$ end do -!!$ end if -!!$ else -!!$ if ((nza >= nzl)) then -!!$ do i=1, nz -!!$ nza = nza + 1 -!!$ if (nza>maxsz) then -!!$ info = -73 -!!$ return -!!$ endif -!!$ aspk(nza) = val(i) -!!$ end do -!!$ else -!!$ do i=1, nz -!!$ ir = ia(i) -!!$ ic = ja(i) -!!$ if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then -!!$ nza = nza + 1 -!!$ if (nza>maxsz) then -!!$ info = -74 -!!$ return -!!$ endif -!!$ aspk(nza) = val(i) -!!$ end if -!!$ end do -!!$ end if -!!$ end if -!!$ end subroutine psb_inner_upd - - subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,maxsz,& - & imin,imax,jmin,jmax,info,gtl) - implicit none - - integer, intent(in) :: nz, imin,imax,jmin,jmax,maxsz - integer, intent(in) :: ia(:),ja(:) - integer, intent(inout) :: nza,ia1(:),ia2(:) - real(psb_dpk_), intent(in) :: val(:) - real(psb_dpk_), intent(inout) :: aspk(:) - integer, intent(out) :: info - integer, intent(in), optional :: gtl(:) - integer :: i,ir,ic,ng - - info = 0 - if (present(gtl)) then - ng = size(gtl) - - 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 >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then - nza = nza + 1 - if (nza > maxsz) then - info = -91 - return - endif - ia1(nza) = ir - ia2(nza) = ic - aspk(nza) = val(i) - end if - end if - end do - else - - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then - nza = nza + 1 - if (nza > maxsz) then - info = -92 - return - endif - ia1(nza) = ir - ia2(nza) = ic - aspk(nza) = val(i) - end if - end do - end if - - end subroutine psb_inner_ins - - -!!$#ifdef FIXED_NAG_SEGV -!!$ subroutine d_coo_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_serial_mod -!!$ implicit none -!!$ -!!$ class(psbn_d_coo_sparse_mat), intent(inout) :: a -!!$ integer, intent(in) :: nz, imin,imax,jmin,jmax -!!$ integer, intent(in) :: ia(:),ja(:) -!!$ real(psb_dpk_), intent(in) :: val(:) -!!$ integer, intent(out) :: info -!!$ integer, intent(in), optional :: gtl(:) -!!$ integer :: i,ir,ic, ilr, ilc, ip, & -!!$ & i1,i2,nc,nnz,dupl,ng -!!$ integer :: debug_level, debug_unit -!!$ character(len=20) :: name='d_coo_srch_upd' -!!$ -!!$ info = 0 -!!$ 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() -!!$ -!!$ -!!$ if (present(gtl)) then -!!$ ng = size(gtl) -!!$ -!!$ select case(dupl) -!!$ case(psbn_dupl_ovwrt_,psbn_dupl_err_) -!!$ ! Overwrite. -!!$ ! Cannot test for error, should have been caught earlier. -!!$ 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) -!!$ if ((ir > 0).and.(ir <= a%m)) then -!!$ ic = gtl(ic) -!!$ if (ir /= ilr) then -!!$ i1 = psb_ibsrch(ir,nnz,a%ia) -!!$ i2 = i1 -!!$ do -!!$ if (i2+1 > nnz) exit -!!$ if (a%ia(i2+1) /= a%ia(i2)) exit -!!$ i2 = i2 + 1 -!!$ end do -!!$ do -!!$ if (i1-1 < 1) exit -!!$ if (a%ia(i1-1) /= a%ia(i1)) exit -!!$ i1 = i1 - 1 -!!$ end do -!!$ ilr = ir -!!$ else -!!$ i1 = 1 -!!$ i2 = 1 -!!$ end if -!!$ nc = i2-i1+1 -!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) -!!$ if (ip>0) then -!!$ 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.' -!!$ endif -!!$ end if -!!$ end do -!!$ case(psbn_dupl_add_) -!!$ ! Add -!!$ 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 <= a%m)) then -!!$ -!!$ if (ir /= ilr) then -!!$ i1 = psb_ibsrch(ir,nnz,a%ia) -!!$ i2 = i1 -!!$ do -!!$ if (i2+1 > nnz) exit -!!$ if (a%ia(i2+1) /= a%ia(i2)) exit -!!$ i2 = i2 + 1 -!!$ end do -!!$ do -!!$ if (i1-1 < 1) exit -!!$ if (a%ia(i1-1) /= a%ia(i1)) exit -!!$ i1 = i1 - 1 -!!$ end do -!!$ ilr = ir -!!$ else -!!$ i1 = 1 -!!$ i2 = 1 -!!$ end if -!!$ nc = i2-i1+1 -!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) -!!$ 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 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(psbn_dupl_ovwrt_,psbn_dupl_err_) -!!$ ! Overwrite. -!!$ ! Cannot test for error, should have been caught earlier. -!!$ do i=1, nz -!!$ ir = ia(i) -!!$ ic = ja(i) -!!$ if ((ir > 0).and.(ir <= a%m)) then -!!$ -!!$ if (ir /= ilr) then -!!$ i1 = psb_ibsrch(ir,nnz,a%ia) -!!$ i2 = i1 -!!$ do -!!$ if (i2+1 > nnz) exit -!!$ if (a%ia(i2+1) /= a%ia(i2)) exit -!!$ i2 = i2 + 1 -!!$ end do -!!$ do -!!$ if (i1-1 < 1) exit -!!$ if (a%ia(i1-1) /= a%ia(i1)) exit -!!$ i1 = i1 - 1 -!!$ end do -!!$ ilr = ir -!!$ else -!!$ i1 = 1 -!!$ i2 = 1 -!!$ end if -!!$ nc = i2-i1+1 -!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) -!!$ if (ip>0) then -!!$ a%val(i1+ip-1) = val(i) -!!$ else -!!$ info = i -!!$ return -!!$ end if -!!$ end if -!!$ end do -!!$ -!!$ case(psbn_dupl_add_) -!!$ ! Add -!!$ do i=1, nz -!!$ ir = ia(i) -!!$ ic = ja(i) -!!$ if ((ir > 0).and.(ir <= a%m)) then -!!$ -!!$ if (ir /= ilr) then -!!$ i1 = psb_ibsrch(ir,nnz,a%ia) -!!$ i2 = i1 -!!$ do -!!$ if (i2+1 > nnz) exit -!!$ if (a%ia(i2+1) /= a%ia(i2)) exit -!!$ i2 = i2 + 1 -!!$ end do -!!$ do -!!$ if (i1-1 < 1) exit -!!$ if (a%ia(i1-1) /= a%ia(i1)) exit -!!$ i1 = i1 - 1 -!!$ end do -!!$ ilr = ir -!!$ else -!!$ i1 = 1 -!!$ i2 = 1 -!!$ end if -!!$ nc = i2-i1+1 -!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) -!!$ if (ip>0) then -!!$ a%val(i1+ip-1) = a%val(i1+ip-1) + val(i) -!!$ else -!!$ info = i -!!$ return -!!$ 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 -!!$ -!!$ end if -!!$ -!!$ end subroutine d_coo_srch_upd -!!$ -!!$#else - subroutine d_coo_srch_upd(nz,ia,ja,val,& - & aia,aja,aval,dupl,nza,nra,& - & info,gtl) - - use psb_error_mod - use psb_sort_mod - implicit none - - integer, intent(inout) :: aia(:),aja(:) - real(psb_dpk_), intent(inout) :: aval(:) - integer, intent(in) :: nz, dupl,nza, nra - integer, intent(in) :: ia(:),ja(:) - real(psb_dpk_), intent(in) :: val(:) - integer, intent(out) :: info - integer, intent(in), optional :: gtl(:) - integer :: i,ir,ic, ilr, ilc, ip, & - & i1,i2,nc,ng - integer :: debug_level, debug_unit - character(len=20) :: name='d_coo_srch_upd' - - info = 0 - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - - - ilr = -1 - ilc = -1 - - - if (present(gtl)) then - ng = size(gtl) - - select case(dupl) - - case(psbn_dupl_ovwrt_,psbn_dupl_err_) - ! Overwrite. - ! Cannot test for error, should have been caught earlier. - 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) - if ((ir > 0).and.(ir <= nra)) then - ic = gtl(ic) - if (ir /= ilr) then - i1 = psb_ibsrch(ir,nza,aia) - i2 = i1 - do - if (i2+1 > nza) exit - if (aia(i2+1) /= aia(i2)) exit - i2 = i2 + 1 - end do - do - if (i1-1 < 1) exit - if (aia(i1-1) /= aia(i1)) exit - i1 = i1 - 1 - end do - ilr = ir - else - i1 = 1 - i2 = 1 - end if - nc = i2-i1+1 - ip = psb_issrch(ic,nc,aja(i1:i2)) - if (ip>0) then - aval(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.' - endif - end if - end do - case(psbn_dupl_add_) - ! Add - 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 <= nra)) then - - if (ir /= ilr) then - i1 = psb_ibsrch(ir,nza,aia) - i2 = i1 - do - if (i2+1 > nza) exit - if (aia(i2+1) /= aia(i2)) exit - i2 = i2 + 1 - end do - do - if (i1-1 < 1) exit - if (aia(i1-1) /= aia(i1)) exit - i1 = i1 - 1 - end do - ilr = ir - else - i1 = 1 - i2 = 1 - end if - nc = i2-i1+1 - ip = psb_issrch(ic,nc,aja(i1:i2)) - if (ip>0) then - aval(i1+ip-1) = aval(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 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(psbn_dupl_ovwrt_,psbn_dupl_err_) - ! Overwrite. - ! Cannot test for error, should have been caught earlier. - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir > 0).and.(ir <= nra)) then - - if (ir /= ilr) then - i1 = psb_ibsrch(ir,nza,aia) - i2 = i1 - do - if (i2+1 > nza) exit - if (aia(i2+1) /= aia(i2)) exit - i2 = i2 + 1 - end do - do - if (i1-1 < 1) exit - if (aia(i1-1) /= aia(i1)) exit - i1 = i1 - 1 - end do - ilr = ir - else - i1 = 1 - i2 = 1 - end if - nc = i2-i1+1 - ip = psb_issrch(ic,nc,aja(i1:i2)) - if (ip>0) then - aval(i1+ip-1) = val(i) - else - info = i - return - end if - end if - end do - - case(psbn_dupl_add_) - ! Add - do i=1, nz - ir = ia(i) - ic = ja(i) - if ((ir > 0).and.(ir <= nra)) then - - if (ir /= ilr) then - i1 = psb_ibsrch(ir,nza,aia) - i2 = i1 - do - if (i2+1 > nza) exit - if (aia(i2+1) /= aia(i2)) exit - i2 = i2 + 1 - end do - do - if (i1-1 < 1) exit - if (aia(i1-1) /= aia(i1)) exit - i1 = i1 - 1 - end do - ilr = ir - else - i1 = 1 - i2 = 1 - end if - nc = i2-i1+1 - ip = psb_issrch(ic,nc,aja(i1:i2)) - if (ip>0) then - aval(i1+ip-1) = aval(i1+ip-1) + val(i) - else - info = i - return - 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 - - end if - - end subroutine d_coo_srch_upd -!!$#endif - end subroutine d_coo_csins - - - subroutine d_coo_csmv(alpha,a,x,beta,y,info,trans) - use psb_error_mod - class(psbn_d_coo_sparse_mat), intent(in) :: a - real(psb_dpk_), intent(in) :: alpha, beta, x(:) - real(psb_dpk_), intent(inout) :: y(:) - integer, intent(out) :: info - character, optional, intent(in) :: trans - - character :: trans_ - integer :: i,j,k,m,n, nnz, ir, jc - real(psb_dpk_) :: acc - logical :: tra - Integer :: err_act - character(len=20) :: name='d_co_csmv' - logical, parameter :: debug=.false. - - call psb_erractionsave(err_act) - - if (present(trans)) then - trans_ = trans - else - trans_ = 'N' - end if - - tra = ((trans_=='T').or.(trans_=='t')) - - if (.not.a%is_asb()) then - info = 1121 - call psb_errpush(info,name) - goto 9999 - endif - - - if (tra) then - m = a%get_ncols() - n = a%get_nrows() - else - n = a%get_ncols() - m = a%get_nrows() - end if - nnz = a%get_nzeros() - - if (alpha == dzero) then - if (beta == dzero) then - do i = 1, m - y(i) = dzero - enddo - else - do i = 1, m - y(i) = beta*y(i) - end do - endif - return - else - if (.not.a%is_unit()) then - if (beta == dzero) then - do i = 1, m - y(i) = dzero - enddo - else - do i = 1, m - y(i) = beta*y(i) - end do - endif - - else if (a%is_unit()) then - if (beta == dzero) then - do i = 1, min(m,n) - y(i) = alpha*x(i) - enddo - do i = min(m,n)+1, m - y(i) = dzero - enddo - else - do i = 1, min(m,n) - y(i) = beta*y(i) + alpha*x(i) - end do - do i = min(m,n)+1, m - y(i) = beta*y(i) - enddo - endif - - endif - - end if - - if (.not.tra) then - i = 1 - j = i - if (nnz > 0) then - ir = a%ia(1) - acc = zero - do - if (i>nnz) then - y(ir) = y(ir) + alpha * acc - exit - endif - if (ia(i) /= ir) then - y(ir) = y(ir) + alpha * acc - ir = ia(i) - acc = zero - endif - acc = acc + a%val(i) * x(a%ja(i)) - i = i + 1 - enddo - end if - - else if (tra) then - if (alpha.eq.done) then - i = 1 - do i=1,nnz - ir = a%ja(i) - jc = a%ia(i) - y(ir) = y(ir) + a%val(i)*x(jc) - enddo - - else if (alpha.eq.-done) then - - do i=1,nnz - ir = a%ja(i) - jc = a%ia(i) - y(ir) = y(ir) - a%val(i)*x(jc) - enddo - - else - - do i=1,nnz - ir = ja(i) - jc = ia(i) - y(ir) = y(ir) + alpha*a%val(i)*x(jc) - enddo - - end if !.....end testing on alpha - - endif - - 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 d_coo_csmv - - subroutine d_coo_csmm(alpha,a,x,beta,y,info,trans) - use psb_error_mod - class(psbn_d_coo_sparse_mat), intent(in) :: a - real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) - real(psb_dpk_), intent(inout) :: y(:,:) - integer, intent(out) :: info - character, optional, intent(in) :: trans - - character :: trans_ - integer :: i,j,k,m,n, nnz, ir, jc, nc - real(psb_dpk_), allocatable :: acc(:) - logical :: tra - Integer :: err_act - character(len=20) :: name='d_coo_csmm' - logical, parameter :: debug=.false. - - call psb_erractionsave(err_act) - - if (present(trans)) then - trans_ = trans - else - trans_ = 'N' - end if - - if (.not.a%is_asb()) then - info = 1121 - call psb_errpush(info,name) - goto 9999 - endif - - - tra = ((trans_=='T').or.(trans_=='t')) - - if (tra) then - m = a%get_ncols() - n = a%get_nrows() - else - n = a%get_ncols() - m = a%get_nrows() - end if - nnz = a%get_nzeros() - - nc = min(size(x,2), size(y,2)) - allocate(acc(nc),stat=info) - if(info /= 0) then - info=4010 - call psb_errpush(info,name,a_err='allocate') - 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 - else - if (.not.a%is_unit()) then - if (beta == dzero) then - do i = 1, m - y(i,:) = dzero - enddo - else - do i = 1, m - y(i,:) = beta*y(i,:) - end do - endif - - else if (a%is_unit()) then - if (beta == dzero) then - do i = 1, min(m,n) - y(i,:) = alpha*x(i,:) - enddo - do i = min(m,n)+1, m - y(i,:) = dzero - enddo - else - do i = 1, min(m,n) - y(i,:) = beta*y(i,:) + alpha*x(i,:) - end do - do i = min(m,n)+1, m - y(i,:) = beta*y(i,:) - enddo - endif - - endif - - end if - - if (.not.tra) then - i = 1 - j = i - if (nnz > 0) then - ir = a%ia(1) - acc = zero - do - if (i>nnz) then - y(ir,:) = y(ir,:) + alpha * acc - exit - endif - if (ia(i) /= ir) then - y(ir,:) = y(ir,:) + alpha * acc - ir = ia(i) - acc = zero - endif - acc = acc + a%val(i) * x(a%ja(i),:) - i = i + 1 - enddo - end if - - else if (tra) then - if (alpha.eq.done) then - i = 1 - do i=1,nnz - ir = a%ja(i) - jc = a%ia(i) - y(ir,:) = y(ir,:) + a%val(i)*x(jc,:) - enddo - - else if (alpha.eq.-done) then - - do i=1,nnz - ir = a%ja(i) - jc = a%ia(i) - y(ir,:) = y(ir,:) - a%val(i)*x(jc,:) - enddo - - else - - do i=1,nnz - ir = ja(i) - jc = ia(i) - y(ir,:) = y(ir,:) + alpha*a%val(i)*x(jc,:) - enddo - - end if !.....end testing on alpha - - endif - - 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 d_coo_csmm - - subroutine d_coo_cssv(alpha,a,x,beta,y,info,trans) - use psb_error_mod - class(psbn_d_coo_sparse_mat), intent(in) :: a - real(psb_dpk_), intent(in) :: alpha, beta, x(:) - real(psb_dpk_), intent(inout) :: y(:) - integer, intent(out) :: info - character, optional, intent(in) :: trans - - character :: trans_ - integer :: i,j,k,m,n, nnz, ir, jc - real(psb_dpk_) :: acc - real(psb_dpk_), allocatable :: tmp(:) - logical :: tra - Integer :: err_act - character(len=20) :: name='d_coo_cssv' - logical, parameter :: debug=.false. - - call psb_erractionsave(err_act) - - if (present(trans)) then - trans_ = trans - else - trans_ = 'N' - end if - if (.not.a%is_asb()) then - info = 1121 - call psb_errpush(info,name) - goto 9999 - endif - - tra = ((trans_=='T').or.(trans_=='t')) - m = a%get_nrows() - - if (.not. (a%is_triangle())) then - info = 1121 - 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_coosv(tra,a,x,y,info) - if (info /= 0) then - call psb_errpush(info,name) - goto 9999 - end if - do i = 1, m - y(i) = alpha*y(i) - end do - else - allocate(tmp(m), stat=info) - if (info /= 0) then - info=4010 - call psb_errpush(info,name,a_err='allocate') - goto 9999 - end if - - tmp(1:m) = x(1:m) - call inner_coosv(tra,a,tmp,y,info) - if (info /= 0) then - call psb_errpush(info,name) - goto 9999 - end if - do i = 1, m - y(i) = alpha*tmp(i) + beta*y(i) - end do - 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_coosv(tra,a,x,y,info) - - logical, intent(in) :: tra - class(psbn_d_coo_sparse_mat), intent(in) :: a - real(psb_dpk_), intent(in) :: x(:) - real(psb_dpk_), intent(out) :: y(:) - - integer :: i,j,k,m, ir, jc, nnz - real(psb_dpk_) :: acc - - if (.not.a%is_sorted()) then - info = 1121 - return - end if - - nnz = a%get_nzeros() - - if (.not.tra) then - - if (a%is_lower()) then - if (a%is_unit()) then - j = 1 - do i=1, a%get_nrows() - acc = dzero - do - if (j > nnz) exit - if (ia(j) > i) exit - acc = acc + a%val(j)*x(a%ja(j)) - j = j + 1 - end do - y(i) = x(i) - acc - end do - else if (.not.a%is_unit()) then - j = 1 - do i=1, a%get_nrows() - acc = dzero - do - if (j > nnz) exit - if (ia(j) > i) exit - if (ja(j) == i) then - y(i) = (x(i) - acc)/a%val(j) - j = j + 1 - exit - end if - acc = acc + a%val(j)*x(a%ja(j)) - j = j + 1 - end do - end do - end if - - else if (a%is_upper()) then - if (a%is_unit()) then - j = nnz - do i=a%get_nrows(), 1, -1 - acc = dzero - do - if (j < 1) exit - if (ia(j) < i) exit - acc = acc + a%val(j)*x(a%ja(j)) - j = j - 1 - end do - y(i) = x(i) - acc - end do - - else if (.not.a%is_unit()) then - - j = nnz - do i=a%get_nrows(), 1, -1 - acc = dzero - do - if (j < 1) exit - if (ia(j) < i) exit - if (ja(j) == i) then - y(i) = (x(i) - acc)/a%val(j) - j = j - 1 - exit - end if - acc = acc + a%val(j)*x(a%ja(j)) - j = j - 1 - end do - end do - end if - - end if - - else if (tra) then - - do i=1, a%get_nrows() - y(i) = x(i) - end do - - if (a%is_lower()) then - if (a%is_unit()) then - j = nnz - do i=a%get_nrows(), 1, -1 - acc = y(i) - do - if (j < 1) exit - if (ia(j) < i) exit - jc = a%ja(j) - y(jc) = y(jc) - a%val(j)*acc - j = j - 1 - end do - end do - else if (.not.a%is_unit()) then - j = nnz - do i=a%get_nrows(), 1, -1 - if (ja(j) == i) then - y(i) = y(i) /a%val(j) - j = j - 1 - end if - acc = y(i) - do - if (j < 1) exit - if (ia(j) < i) exit - jc = a%ja(j) - y(jc) = y(jc) - a%val(j)*acc - j = j - 1 - end do - end do - - else if (a%is_upper()) then - if (a%is_unit()) then - j = 1 - do i=1, a%get_nrows() - acc = y(i) - do - if (j > nnz) exit - if (ia(j) > i) exit - jc = a%ja(j) - y(jc) = y(jc) - a%val(j)*acc - j = j + 1 - end do - end do - else if (.not.a%is_unit()) then - j = 1 - do i=1, a%get_nrows() - if (ja(j) == i) then - y(i) = y(i) /a%val(j) - j = j + 1 - end if - acc = y(i) - do - if (j > nnz) exit - if (ia(j) > i) exit - jc = a%ja(j) - y(jc) = y(jc) - a%val(j)*acc - j = j + 1 - end do - end do - end if - end if - end if - end if - - end subroutine inner_coosv - - end subroutine d_coo_cssv - - - - subroutine d_coo_cssm(alpha,a,x,beta,y,info,trans) - use psb_error_mod - class(psbn_d_coo_sparse_mat), intent(in) :: a - real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) - real(psb_dpk_), intent(inout) :: y(:,:) - integer, intent(out) :: info - character, optional, intent(in) :: trans - - character :: trans_ - integer :: i,j,k,m,n, nnz, ir, jc, nc - real(psb_dpk_) :: acc - real(psb_dpk_), allocatable :: tmp(:,:) - logical :: tra - Integer :: err_act - character(len=20) :: name='d_base_csmm' - logical, parameter :: debug=.false. - - call psb_erractionsave(err_act) - - if (present(trans)) then - trans_ = trans - else - trans_ = 'N' - end if - if (.not.a%is_asb()) then - info = 1121 - call psb_errpush(info,name) - goto 9999 - endif - - - tra = ((trans_=='T').or.(trans_=='t')) - m = a%get_nrows() - nc = min(size(x,2) , size(y,2)) - - if (.not. (a%is_triangle())) then - write(0,*) 'Called SM on a non-triangular mat!' - info = 1121 - 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_coosm(tra,a,x,y,info) - do i = 1, m - y(i,:) = alpha*y(i,:) - end do - else - allocate(tmp(m,nc), stat=info) - if(info /= 0) then - info=4010 - call psb_errpush(info,name,a_err='allocate') - goto 9999 - end if - - tmp(1:m,:) = x(1:m,:) - call inner_coosm(tra,a,tmp,y,info) - do i = 1, m - y(i,:) = alpha*tmp(i,:) + beta*y(i,:) - end do - end if - - if(info /= 0) then - info=4010 - call psb_errpush(info,name,a_err='inner_coosm') - 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_coosm(tra,a,x,y,info) - logical, intent(in) :: tra - class(psbn_d_coo_sparse_mat), intent(in) :: a - real(psb_dpk_), intent(in) :: x(:,:) - real(psb_dpk_), intent(out) :: y(:,:) - integer, intent(out) :: info - integer :: i,j,k,m, ir, jc - real(psb_dpk_), allocatable :: acc(:) - - allocate(acc(size(x,2)), stat=info) - if(info /= 0) then - info=4010 - return - end if - - - if (.not.a%is_sorted()) then - info = 1121 - return - end if - - nnz = a%get_nzeros() - - if (.not.tra) then - - if (a%is_lower()) then - if (a%is_unit()) then - j = 1 - do i=1, a%get_nrows() - acc = dzero - do - if (j > nnz) exit - if (ia(j) > i) exit - acc = acc + a%val(j)*x(a%ja(j),:) - j = j + 1 - end do - y(i,:) = x(i,:) - acc - end do - else if (.not.a%is_unit()) then - j = 1 - do i=1, a%get_nrows() - acc = dzero - do - if (j > nnz) exit - if (ia(j) > i) exit - if (ja(j) == i) then - y(i,:) = (x(i,:) - acc)/a%val(j) - j = j + 1 - exit - end if - acc = acc + a%val(j)*x(a%ja(j),:) - j = j + 1 - end do - end do - end if - - else if (a%is_upper()) then - if (a%is_unit()) then - j = nnz - do i=a%get_nrows(), 1, -1 - acc = dzero - do - if (j < 1) exit - if (ia(j) < i) exit - acc = acc + a%val(j)*x(a%ja(j),:) - j = j - 1 - end do - y(i,:) = x(i,:) - acc - end do - - else if (.not.a%is_unit()) then - - j = nnz - do i=a%get_nrows(), 1, -1 - acc = dzero - do - if (j < 1) exit - if (ia(j) < i) exit - if (ja(j) == i) then - y(i,:) = (x(i,:) - acc)/a%val(j) - j = j - 1 - exit - end if - acc = acc + a%val(j)*x(a%ja(j),:) - j = j - 1 - end do - end do - end if - - end if - - else if (tra) then - - do i=1, a%get_nrows() - y(i,:) = x(i,:) - end do - - if (a%is_lower()) then - if (a%is_unit()) then - j = nnz - do i=a%get_nrows(), 1, -1 - acc = y(i,:) - do - if (j < 1) exit - if (ia(j) < i) exit - jc = a%ja(j) - y(jc,:) = y(jc,:) - a%val(j)*acc - j = j - 1 - end do - end do - else if (.not.a%is_unit()) then - j = nnz - do i=a%get_nrows(), 1, -1 - if (ja(j) == i) then - y(i,:) = y(i,:) /a%val(j) - j = j - 1 - end if - acc = y(i,:) - do - if (j < 1) exit - if (ia(j) < i) exit - jc = a%ja(j) - y(jc,:) = y(jc,:) - a%val(j)*acc - j = j - 1 - end do - end do - - else if (a%is_upper()) then - if (a%is_unit()) then - j = 1 - do i=1, a%get_nrows() - acc = y(i,:) - do - if (j > nnz) exit - if (ia(j) > i) exit - jc = a%ja(j) - y(jc,:) = y(jc,:) - a%val(j)*acc - j = j + 1 - end do - end do - else if (.not.a%is_unit()) then - j = 1 - do i=1, a%get_nrows() - if (ja(j) == i) then - y(i,:) = y(i,:) /a%val(j) - j = j + 1 - end if - acc = y(i,:) - do - if (j > nnz) exit - if (ia(j) > i) exit - jc = a%ja(j) - y(jc,:) = y(jc,:) - a%val(j)*acc - j = j + 1 - end do - end do - end if - end if - end if - end if - end subroutine inner_coosm - - end subroutine d_coo_cssm - -end module psbn_d_coo_sparse_mat_mod - diff --git a/base/newserial/psbn_d_base_mat_mod.f03 b/base/newserial/psbn_d_base_mat_mod.f03 index 6e1fab84..f58c3046 100644 --- a/base/newserial/psbn_d_base_mat_mod.f03 +++ b/base/newserial/psbn_d_base_mat_mod.f03 @@ -12,11 +12,269 @@ module psbn_d_base_mat_mod procedure, pass(a) :: d_base_cssm generic, public :: psbn_cssm => d_base_cssm, d_base_cssv procedure, pass(a) :: csins - + procedure, pass(a) :: to_coo + procedure, pass(a) :: from_coo + end type psbn_d_base_sparse_mat + type, extends(psbn_d_base_sparse_mat) :: psbn_d_coo_sparse_mat + + integer :: nnz + integer, allocatable :: ia(:), ja(:) + real(psb_dpk_), allocatable :: val(:) + + contains + + procedure, pass(a) :: get_nzeros => d_coo_get_nzeros + procedure, pass(a) :: set_nzeros => d_coo_set_nzeros + procedure, pass(a) :: d_base_csmm => d_coo_csmm + procedure, pass(a) :: d_base_csmv => d_coo_csmv + procedure, pass(a) :: d_base_cssm => d_coo_cssm + procedure, pass(a) :: d_base_cssv => d_coo_cssv + procedure, pass(a) :: csins => d_coo_csins + procedure, pass(a) :: reallocate_nz => d_coo_reallocate_nz + procedure, pass(a) :: to_coo => d_coo_to_coo + procedure, pass(a) :: from_coo => d_coo_from_coo + procedure, pass(a) :: fix => d_fix_coo + + end type psbn_d_coo_sparse_mat + + + interface + subroutine d_fix_coo_impl(a,info,idir) + use psb_const_mod + import psbn_d_coo_sparse_mat + class(psbn_d_coo_sparse_mat), intent(inout) :: a + integer, intent(out) :: info + integer, intent(in), optional :: idir + end subroutine d_fix_coo_impl + end interface + + interface + subroutine d_coo_to_coo_impl(a,b,info) + use psb_const_mod + import psbn_d_coo_sparse_mat + class(psbn_d_coo_sparse_mat), intent(in) :: a + class(psbn_d_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + end subroutine d_coo_to_coo_impl + end interface + + interface + subroutine d_coo_from_coo_impl(a,b,info) + use psb_const_mod + import psbn_d_coo_sparse_mat + class(psbn_d_coo_sparse_mat), intent(inout) :: a + class(psbn_d_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info + end subroutine d_coo_from_coo_impl + end interface + + interface + subroutine d_coo_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) + use psb_const_mod + import psbn_d_coo_sparse_mat + class(psbn_d_coo_sparse_mat), intent(inout) :: a + real(psb_dpk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + end subroutine d_coo_csins_impl + end interface + + interface d_coo_cssm_impl + subroutine d_coo_cssv_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psbn_d_coo_sparse_mat + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine d_coo_cssv_impl + subroutine d_coo_cssm_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psbn_d_coo_sparse_mat + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine d_coo_cssm_impl + end interface + + interface d_coo_csmm_impl + subroutine d_coo_csmv_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psbn_d_coo_sparse_mat + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine d_coo_csmv_impl + subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + import psbn_d_coo_sparse_mat + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + end subroutine d_coo_csmm_impl + end interface + contains + + + subroutine to_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + class(psbn_d_base_sparse_mat), intent(in) :: a + class(psbn_d_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine to_coo + + + subroutine d_fix_coo(a,info,idir) + use psb_error_mod + use psb_const_mod + class(psbn_d_coo_sparse_mat), intent(inout) :: a + integer, intent(out) :: info + integer, intent(in), optional :: idir + Integer :: err_act + character(len=20) :: name='fix_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call d_fix_coo_impl(a,info,idir) + + if (info /= 0) 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 d_fix_coo + + subroutine d_coo_to_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + class(psbn_d_coo_sparse_mat), intent(in) :: a + class(psbn_d_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call d_coo_to_coo_impl(a,b,info) + if (info /= 0) 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 d_coo_to_coo + + subroutine d_coo_from_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + class(psbn_d_coo_sparse_mat), intent(inout) :: a + class(psbn_d_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + info = 0 + call d_coo_from_coo_impl(a,b,info) + if (info /= 0) 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 d_coo_from_coo + + subroutine from_coo(a,b,info) + use psb_error_mod + use psb_realloc_mod + class(psbn_d_base_sparse_mat), intent(inout) :: a + class(psbn_d_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + info = 700 + call psb_errpush(info,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine from_coo + + + subroutine csins(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) use psb_error_mod use psb_realloc_mod @@ -150,4 +408,305 @@ contains end subroutine d_base_cssv + subroutine d_coo_reallocate_nz(nz,a) + use psb_error_mod + use psb_realloc_mod + integer, intent(in) :: nz + class(psbn_d_coo_sparse_mat), intent(inout) :: a + Integer :: err_act + character(len=20) :: name='d_coo_reallocate_nz' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + call psb_realloc(nx,a%ia,a%ja,a%val,info) + + if (info /= 0) then + call psb_errpush(4000,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 d_coo_reallocate_nz + + + function d_coo_get_nzeros(a) result(res) + class(psbn_d_coo_sparse_mat), intent(in) :: a + integer :: res + res = a%nnz + end function d_coo_get_nzeros + + + subroutine d_coo_set_nzeros(nz,a) + integer, intent(in) :: nz + class(psbn_d_coo_sparse_mat), intent(inout) :: a + + a%nnz = nz + + end subroutine d_coo_set_nzeros + + + subroutine d_coo_csins(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) + use psb_error_mod + use psb_realloc_mod + class(psbn_d_coo_sparse_mat), intent(inout) :: a + real(psb_dpk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + + + Integer :: err_act + character(len=20) :: name='d_coo_csins' + logical, parameter :: debug=.false. + integer :: nza, i,j,k, nzl, isza, int_err(5) + + call psb_erractionsave(err_act) + info = 0 + + if (nz <= 0) then + info = 10 + int_err(1)=1 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + if (size(ia) < nz) then + info = 35 + int_err(1)=2 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + if (size(ja) < nz) then + info = 35 + int_err(1)=3 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + if (size(val) < nz) then + info = 35 + int_err(1)=4 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + if (nz == 0) return + + call d_coo_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) + if (info /= 0) 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 d_coo_csins + + + subroutine d_coo_csmv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc + real(psb_dpk_) :: acc + logical :: tra + Integer :: err_act + character(len=20) :: name='d_coo_csmv' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + call d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) + + if (info /= 0) 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 d_coo_csmv + + subroutine d_coo_csmm(alpha,a,x,beta,y,info,trans) + use psb_error_mod + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_dpk_), allocatable :: acc(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='d_coo_csmm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + + + call d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) + + if (info /= 0) 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 d_coo_csmm + + + subroutine d_coo_cssv(alpha,a,x,beta,y,info,trans) + use psb_error_mod + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc + real(psb_dpk_) :: acc + real(psb_dpk_), allocatable :: tmp(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='d_coo_cssv' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + if (.not. (a%is_triangle())) then + write(0,*) 'Called SM on a non-triangular mat!' + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + call d_coo_cssm_impl(alpha,a,x,beta,y,info,trans) + + 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 d_coo_cssv + + + + subroutine d_coo_cssm(alpha,a,x,beta,y,info,trans) + use psb_error_mod + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_dpk_) :: acc + real(psb_dpk_), allocatable :: tmp(:,:) + logical :: tra + Integer :: err_act + character(len=20) :: name='d_coo_csmm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + if (.not. (a%is_triangle())) then + write(0,*) 'Called SM on a non-triangular mat!' + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + call d_coo_cssm_impl(alpha,a,x,beta,y,info,trans) + 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 d_coo_cssm + + end module psbn_d_base_mat_mod + + + diff --git a/base/newserial/psbn_d_coo_impl.f03 b/base/newserial/psbn_d_coo_impl.f03 new file mode 100644 index 00000000..3c08ab5e --- /dev/null +++ b/base/newserial/psbn_d_coo_impl.f03 @@ -0,0 +1,1845 @@ + +subroutine d_coo_cssm_impl(alpha,a,x,beta,y,info,trans) + use psb_error_mod + use psbn_d_base_mat_mod, psb_protect_name => d_coo_cssm_impl + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_dpk_) :: acc + real(psb_dpk_), allocatable :: tmp(:,:) + logical :: tra + Integer :: err_act + character(len=20) :: name='d_base_csmm' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + if (.not. (a%is_triangle())) then + write(0,*) 'Called SM on a non-triangular mat!' + info = 1121 + call psb_errpush(info,name) + goto 9999 + end if + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + tra = ((trans_=='T').or.(trans_=='t')) + m = a%get_nrows() + nc = min(size(x,2) , size(y,2)) + + + 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_coosm(tra,a,x,y,info) + do i = 1, m + y(i,:) = alpha*y(i,:) + end do + else + allocate(tmp(m,nc), stat=info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='allocate') + goto 9999 + end if + + tmp(1:m,:) = x(1:m,:) + call inner_coosm(tra,a,tmp,y,info) + do i = 1, m + y(i,:) = alpha*tmp(i,:) + beta*y(i,:) + end do + end if + + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='inner_coosm') + 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_coosm(tra,a,x,y,info) + logical, intent(in) :: tra + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: x(:,:) + real(psb_dpk_), intent(out) :: y(:,:) + integer, intent(out) :: info + integer :: i,j,k,m, ir, jc + real(psb_dpk_), allocatable :: acc(:) + + allocate(acc(size(x,2)), stat=info) + if(info /= 0) then + info=4010 + return + end if + + + if (.not.a%is_sorted()) then + info = 1121 + return + end if + + nnz = a%get_nzeros() + + if (.not.tra) then + + if (a%is_lower()) then + if (a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + acc = dzero + do + if (j > nnz) exit + if (ia(j) > i) exit + acc = acc + a%val(j)*x(a%ja(j),:) + j = j + 1 + end do + y(i,:) = x(i,:) - acc + end do + else if (.not.a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + acc = dzero + do + if (j > nnz) exit + if (ia(j) > i) exit + if (ja(j) == i) then + y(i,:) = (x(i,:) - acc)/a%val(j) + j = j + 1 + exit + end if + acc = acc + a%val(j)*x(a%ja(j),:) + j = j + 1 + end do + end do + end if + + else if (a%is_upper()) then + if (a%is_unit()) then + j = nnz + do i=a%get_nrows(), 1, -1 + acc = dzero + do + if (j < 1) exit + if (ia(j) < i) exit + acc = acc + a%val(j)*x(a%ja(j),:) + j = j - 1 + end do + y(i,:) = x(i,:) - acc + end do + + else if (.not.a%is_unit()) then + + j = nnz + do i=a%get_nrows(), 1, -1 + acc = dzero + do + if (j < 1) exit + if (ia(j) < i) exit + if (ja(j) == i) then + y(i,:) = (x(i,:) - acc)/a%val(j) + j = j - 1 + exit + end if + acc = acc + a%val(j)*x(a%ja(j),:) + j = j - 1 + end do + end do + end if + + end if + + else if (tra) then + + do i=1, a%get_nrows() + y(i,:) = x(i,:) + end do + + if (a%is_lower()) then + if (a%is_unit()) then + j = nnz + do i=a%get_nrows(), 1, -1 + acc = y(i,:) + do + if (j < 1) exit + if (ia(j) < i) exit + jc = a%ja(j) + y(jc,:) = y(jc,:) - a%val(j)*acc + j = j - 1 + end do + end do + else if (.not.a%is_unit()) then + j = nnz + do i=a%get_nrows(), 1, -1 + if (ja(j) == i) then + y(i,:) = y(i,:) /a%val(j) + j = j - 1 + end if + acc = y(i,:) + do + if (j < 1) exit + if (ia(j) < i) exit + jc = a%ja(j) + y(jc,:) = y(jc,:) - a%val(j)*acc + j = j - 1 + end do + end do + + else if (a%is_upper()) then + if (a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + acc = y(i,:) + do + if (j > nnz) exit + if (ia(j) > i) exit + jc = a%ja(j) + y(jc,:) = y(jc,:) - a%val(j)*acc + j = j + 1 + end do + end do + else if (.not.a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + if (ja(j) == i) then + y(i,:) = y(i,:) /a%val(j) + j = j + 1 + end if + acc = y(i,:) + do + if (j > nnz) exit + if (ia(j) > i) exit + jc = a%ja(j) + y(jc,:) = y(jc,:) - a%val(j)*acc + j = j + 1 + end do + end do + end if + end if + end if + end if + end subroutine inner_coosm + +end subroutine d_coo_cssm_impl + + + +subroutine d_coo_cssv_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + use psb_error_mod + use psbn_d_base_mat_mod, psb_protect_name => d_coo_cssv_impl + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc + real(psb_dpk_) :: acc + real(psb_dpk_), allocatable :: tmp(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='d_coo_cssv_impl' + logical, parameter :: debug=.false. + + + call psb_erractionsave(err_act) + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + tra = ((trans_=='T').or.(trans_=='t')) + m = a%get_nrows() + + if (.not. (a%is_triangle())) then + info = 1121 + 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_coosv(tra,a,x,y,info) + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + do i = 1, m + y(i) = alpha*y(i) + end do + else + allocate(tmp(m), stat=info) + if (info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='allocate') + goto 9999 + end if + + tmp(1:m) = x(1:m) + call inner_coosv(tra,a,tmp,y,info) + if (info /= 0) then + call psb_errpush(info,name) + goto 9999 + end if + do i = 1, m + y(i) = alpha*tmp(i) + beta*y(i) + end do + 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_coosv(tra,a,x,y,info) + + logical, intent(in) :: tra + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: x(:) + real(psb_dpk_), intent(out) :: y(:) + + integer :: i,j,k,m, ir, jc, nnz + real(psb_dpk_) :: acc + + if (.not.a%is_sorted()) then + info = 1121 + return + end if + + nnz = a%get_nzeros() + + if (.not.tra) then + + if (a%is_lower()) then + if (a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + acc = dzero + do + if (j > nnz) exit + if (ia(j) > i) exit + acc = acc + a%val(j)*x(a%ja(j)) + j = j + 1 + end do + y(i) = x(i) - acc + end do + else if (.not.a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + acc = dzero + do + if (j > nnz) exit + if (ia(j) > i) exit + if (ja(j) == i) then + y(i) = (x(i) - acc)/a%val(j) + j = j + 1 + exit + end if + acc = acc + a%val(j)*x(a%ja(j)) + j = j + 1 + end do + end do + end if + + else if (a%is_upper()) then + if (a%is_unit()) then + j = nnz + do i=a%get_nrows(), 1, -1 + acc = dzero + do + if (j < 1) exit + if (ia(j) < i) exit + acc = acc + a%val(j)*x(a%ja(j)) + j = j - 1 + end do + y(i) = x(i) - acc + end do + + else if (.not.a%is_unit()) then + + j = nnz + do i=a%get_nrows(), 1, -1 + acc = dzero + do + if (j < 1) exit + if (ia(j) < i) exit + if (ja(j) == i) then + y(i) = (x(i) - acc)/a%val(j) + j = j - 1 + exit + end if + acc = acc + a%val(j)*x(a%ja(j)) + j = j - 1 + end do + end do + end if + + end if + + else if (tra) then + + do i=1, a%get_nrows() + y(i) = x(i) + end do + + if (a%is_lower()) then + if (a%is_unit()) then + j = nnz + do i=a%get_nrows(), 1, -1 + acc = y(i) + do + if (j < 1) exit + if (ia(j) < i) exit + jc = a%ja(j) + y(jc) = y(jc) - a%val(j)*acc + j = j - 1 + end do + end do + else if (.not.a%is_unit()) then + j = nnz + do i=a%get_nrows(), 1, -1 + if (ja(j) == i) then + y(i) = y(i) /a%val(j) + j = j - 1 + end if + acc = y(i) + do + if (j < 1) exit + if (ia(j) < i) exit + jc = a%ja(j) + y(jc) = y(jc) - a%val(j)*acc + j = j - 1 + end do + end do + + else if (a%is_upper()) then + if (a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + acc = y(i) + do + if (j > nnz) exit + if (ia(j) > i) exit + jc = a%ja(j) + y(jc) = y(jc) - a%val(j)*acc + j = j + 1 + end do + end do + else if (.not.a%is_unit()) then + j = 1 + do i=1, a%get_nrows() + if (ja(j) == i) then + y(i) = y(i) /a%val(j) + j = j + 1 + end if + acc = y(i) + do + if (j > nnz) exit + if (ia(j) > i) exit + jc = a%ja(j) + y(jc) = y(jc) - a%val(j)*acc + j = j + 1 + end do + end do + end if + end if + end if + end if + + end subroutine inner_coosv + + +end subroutine d_coo_cssv_impl + +subroutine d_coo_csmv_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + use psb_error_mod + use psbn_d_base_mat_mod, psb_protect_name => d_coo_csMv_impl + + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:) + real(psb_dpk_), intent(inout) :: y(:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc + real(psb_dpk_) :: acc + logical :: tra + Integer :: err_act + character(len=20) :: name='d_coo_csmv_impl' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + + tra = ((trans_=='T').or.(trans_=='t')) + + + + if (tra) then + m = a%get_ncols() + n = a%get_nrows() + else + n = a%get_ncols() + m = a%get_nrows() + end if + nnz = a%get_nzeros() + + if (alpha == dzero) then + if (beta == dzero) then + do i = 1, m + y(i) = dzero + enddo + else + do i = 1, m + y(i) = beta*y(i) + end do + endif + return + else + if (.not.a%is_unit()) then + if (beta == dzero) then + do i = 1, m + y(i) = dzero + enddo + else + do i = 1, m + y(i) = beta*y(i) + end do + endif + + else if (a%is_unit()) then + if (beta == dzero) then + do i = 1, min(m,n) + y(i) = alpha*x(i) + enddo + do i = min(m,n)+1, m + y(i) = dzero + enddo + else + do i = 1, min(m,n) + y(i) = beta*y(i) + alpha*x(i) + end do + do i = min(m,n)+1, m + y(i) = beta*y(i) + enddo + endif + + endif + + end if + + if (.not.tra) then + i = 1 + j = i + if (nnz > 0) then + ir = a%ia(1) + acc = zero + do + if (i>nnz) then + y(ir) = y(ir) + alpha * acc + exit + endif + if (ia(i) /= ir) then + y(ir) = y(ir) + alpha * acc + ir = ia(i) + acc = zero + endif + acc = acc + a%val(i) * x(a%ja(i)) + i = i + 1 + enddo + end if + + else if (tra) then + if (alpha.eq.done) then + i = 1 + do i=1,nnz + ir = a%ja(i) + jc = a%ia(i) + y(ir) = y(ir) + a%val(i)*x(jc) + enddo + + else if (alpha.eq.-done) then + + do i=1,nnz + ir = a%ja(i) + jc = a%ia(i) + y(ir) = y(ir) - a%val(i)*x(jc) + enddo + + else + + do i=1,nnz + ir = ja(i) + jc = ia(i) + y(ir) = y(ir) + alpha*a%val(i)*x(jc) + enddo + + end if !.....end testing on alpha + + endif + + 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 d_coo_csmv_impl + + +subroutine d_coo_csmm_impl(alpha,a,x,beta,y,info,trans) + use psb_const_mod + use psb_error_mod + use psbn_d_base_mat_mod, psb_protect_name => d_coo_csmm_impl + class(psbn_d_coo_sparse_mat), intent(in) :: a + real(psb_dpk_), intent(in) :: alpha, beta, x(:,:) + real(psb_dpk_), intent(inout) :: y(:,:) + integer, intent(out) :: info + character, optional, intent(in) :: trans + + character :: trans_ + integer :: i,j,k,m,n, nnz, ir, jc, nc + real(psb_dpk_), allocatable :: acc(:) + logical :: tra + Integer :: err_act + character(len=20) :: name='d_coo_csmm_impl' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + + + if (.not.a%is_asb()) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + + if (present(trans)) then + trans_ = trans + else + trans_ = 'N' + end if + + + tra = ((trans_=='T').or.(trans_=='t')) + + if (tra) then + m = a%get_ncols() + n = a%get_nrows() + else + n = a%get_ncols() + m = a%get_nrows() + end if + nnz = a%get_nzeros() + + nc = min(size(x,2), size(y,2)) + allocate(acc(nc),stat=info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='allocate') + 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 + else + if (.not.a%is_unit()) then + if (beta == dzero) then + do i = 1, m + y(i,:) = dzero + enddo + else + do i = 1, m + y(i,:) = beta*y(i,:) + end do + endif + + else if (a%is_unit()) then + if (beta == dzero) then + do i = 1, min(m,n) + y(i,:) = alpha*x(i,:) + enddo + do i = min(m,n)+1, m + y(i,:) = dzero + enddo + else + do i = 1, min(m,n) + y(i,:) = beta*y(i,:) + alpha*x(i,:) + end do + do i = min(m,n)+1, m + y(i,:) = beta*y(i,:) + enddo + endif + + endif + + end if + + if (.not.tra) then + i = 1 + j = i + if (nnz > 0) then + ir = a%ia(1) + acc = zero + do + if (i>nnz) then + y(ir,:) = y(ir,:) + alpha * acc + exit + endif + if (ia(i) /= ir) then + y(ir,:) = y(ir,:) + alpha * acc + ir = ia(i) + acc = zero + endif + acc = acc + a%val(i) * x(a%ja(i),:) + i = i + 1 + enddo + end if + + else if (tra) then + if (alpha.eq.done) then + i = 1 + do i=1,nnz + ir = a%ja(i) + jc = a%ia(i) + y(ir,:) = y(ir,:) + a%val(i)*x(jc,:) + enddo + + else if (alpha.eq.-done) then + + do i=1,nnz + ir = a%ja(i) + jc = a%ia(i) + y(ir,:) = y(ir,:) - a%val(i)*x(jc,:) + enddo + + else + + do i=1,nnz + ir = ja(i) + jc = ia(i) + y(ir,:) = y(ir,:) + alpha*a%val(i)*x(jc,:) + enddo + + end if !.....end testing on alpha + + endif + + 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 d_coo_csmm_impl + + +subroutine d_coo_csins_impl(nz,val,ia,ja,a,imin,imax,jmin,jmax,info,gtl) + use psb_error_mod + use psb_realloc_mod + use psbn_d_base_mat_mod, psb_protect_name => d_coo_csins_impl + + class(psbn_d_coo_sparse_mat), intent(inout) :: a + real(psb_dpk_), intent(in) :: val(:) + integer, intent(in) :: nz, ia(:), ja(:), imin,imax,jmin,jmax + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + + + Integer :: err_act + character(len=20) :: name='d_coo_csins' + logical, parameter :: debug=.false. + integer :: nza, i,j,k, nzl, isza, int_err(5) + + call psb_erractionsave(err_act) + info = 0 + + if (nz <= 0) then + info = 10 + int_err(1)=1 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + if (size(ia) < nz) then + info = 35 + int_err(1)=2 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + if (size(ja) < nz) then + info = 35 + int_err(1)=3 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + if (size(val) < nz) then + info = 35 + int_err(1)=4 + call psb_errpush(info,name,i_err=int_err) + goto 9999 + end if + + if (nz == 0) return + + + nza = a%get_nzeros() + isza = a%get_size() + if (a%is_bld()) then + ! Build phase. Must handle reallocations in a sensible way. + if (isza < (nza+nz)) then + call a%reallocate(max(nza+nz,int(1.5*isza))) + isza = a%get_size() + endif + + call psb_inner_ins(nz,ia,ja,val,nza,a%ia,a%ja,a%val,isza,& + & imin,imax,jmin,jmax,info,gtl) + call a%set_nzeros(nz+nza) + + else if (a%is_upd()) then + if (a%is_sorted()) then + + +!!$#ifdef FIXED_NAG_SEGV +!!$ call d_coo_srch_upd(nz,ia,ja,val,a,& +!!$ & imin,imax,jmin,jmax,info,gtl) +!!$#else + call d_coo_srch_upd(nz,ia,ja,val,& + & a%ia,a%ja,a%val,& + & a%get_dupl(),a%get_nzeros(),a%get_nrows(),& + & info,gtl) +!!$#endif + + else + info = 1121 + end if + + else + ! State is wrong. + info = 1121 + end if + if (info /= 0) 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_inner_upd(nz,ia,ja,val,nza,aspk,maxsz,& +!!$ & imin,imax,jmin,jmax,info,gtl) +!!$ implicit none +!!$ +!!$ integer, intent(in) :: nz, imin,imax,jmin,jmax,maxsz +!!$ integer, intent(in) :: ia(:),ja(:) +!!$ integer, intent(inout) :: nza +!!$ real(psb_dpk_), intent(in) :: val(:) +!!$ real(psb_dpk_), intent(inout) :: aspk(:) +!!$ integer, intent(out) :: info +!!$ integer, intent(in), optional :: gtl(:) +!!$ integer :: i,ir,ic, ng,nzl +!!$ character(len=20) :: name, ch_err +!!$ +!!$ +!!$ name='psb_inner_upd' +!!$ nzl = 0 +!!$ if (present(gtl)) then +!!$ ng = size(gtl) +!!$ if ((nza > nzl)) then +!!$ do i=1, nz +!!$ nza = nza + 1 +!!$ if (nza>maxsz) then +!!$ call psb_errpush(50,name,i_err=(/7,maxsz,5,0,nza /)) +!!$ info = -71 +!!$ return +!!$ endif +!!$ aspk(nza) = val(i) +!!$ end do +!!$ else +!!$ 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 >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then +!!$ nza = nza + 1 +!!$ if (nza>maxsz) then +!!$ info = -72 +!!$ return +!!$ endif +!!$ aspk(nza) = val(i) +!!$ end if +!!$ end if +!!$ end do +!!$ end if +!!$ else +!!$ if ((nza >= nzl)) then +!!$ do i=1, nz +!!$ nza = nza + 1 +!!$ if (nza>maxsz) then +!!$ info = -73 +!!$ return +!!$ endif +!!$ aspk(nza) = val(i) +!!$ end do +!!$ else +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then +!!$ nza = nza + 1 +!!$ if (nza>maxsz) then +!!$ info = -74 +!!$ return +!!$ endif +!!$ aspk(nza) = val(i) +!!$ end if +!!$ end do +!!$ end if +!!$ end if +!!$ end subroutine psb_inner_upd + + subroutine psb_inner_ins(nz,ia,ja,val,nza,ia1,ia2,aspk,maxsz,& + & imin,imax,jmin,jmax,info,gtl) + implicit none + + integer, intent(in) :: nz, imin,imax,jmin,jmax,maxsz + integer, intent(in) :: ia(:),ja(:) + integer, intent(inout) :: nza,ia1(:),ia2(:) + real(psb_dpk_), intent(in) :: val(:) + real(psb_dpk_), intent(inout) :: aspk(:) + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + integer :: i,ir,ic,ng + + info = 0 + if (present(gtl)) then + ng = size(gtl) + + 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 >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then + nza = nza + 1 + if (nza > maxsz) then + info = -91 + return + endif + ia1(nza) = ir + ia2(nza) = ic + aspk(nza) = val(i) + end if + end if + end do + else + + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir >=imin).and.(ir<=imax).and.(ic>=jmin).and.(ic<=jmax)) then + nza = nza + 1 + if (nza > maxsz) then + info = -92 + return + endif + ia1(nza) = ir + ia2(nza) = ic + aspk(nza) = val(i) + end if + end do + end if + + end subroutine psb_inner_ins + + +!!$#ifdef FIXED_NAG_SEGV +!!$ subroutine d_coo_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_serial_mod +!!$ implicit none +!!$ +!!$ class(psbn_d_coo_sparse_mat), intent(inout) :: a +!!$ integer, intent(in) :: nz, imin,imax,jmin,jmax +!!$ integer, intent(in) :: ia(:),ja(:) +!!$ real(psb_dpk_), intent(in) :: val(:) +!!$ integer, intent(out) :: info +!!$ integer, intent(in), optional :: gtl(:) +!!$ integer :: i,ir,ic, ilr, ilc, ip, & +!!$ & i1,i2,nc,nnz,dupl,ng +!!$ integer :: debug_level, debug_unit +!!$ character(len=20) :: name='d_coo_srch_upd' +!!$ +!!$ info = 0 +!!$ 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() +!!$ +!!$ +!!$ if (present(gtl)) then +!!$ ng = size(gtl) +!!$ +!!$ select case(dupl) +!!$ case(psbn_dupl_ovwrt_,psbn_dupl_err_) +!!$ ! Overwrite. +!!$ ! Cannot test for error, should have been caught earlier. +!!$ 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) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ ic = gtl(ic) +!!$ if (ir /= ilr) then +!!$ i1 = psb_ibsrch(ir,nnz,a%ia) +!!$ i2 = i1 +!!$ do +!!$ if (i2+1 > nnz) exit +!!$ if (a%ia(i2+1) /= a%ia(i2)) exit +!!$ i2 = i2 + 1 +!!$ end do +!!$ do +!!$ if (i1-1 < 1) exit +!!$ if (a%ia(i1-1) /= a%ia(i1)) exit +!!$ i1 = i1 - 1 +!!$ end do +!!$ ilr = ir +!!$ else +!!$ i1 = 1 +!!$ i2 = 1 +!!$ end if +!!$ nc = i2-i1+1 +!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) +!!$ if (ip>0) then +!!$ 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.' +!!$ endif +!!$ end if +!!$ end do +!!$ case(psbn_dupl_add_) +!!$ ! Add +!!$ 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 <= a%m)) then +!!$ +!!$ if (ir /= ilr) then +!!$ i1 = psb_ibsrch(ir,nnz,a%ia) +!!$ i2 = i1 +!!$ do +!!$ if (i2+1 > nnz) exit +!!$ if (a%ia(i2+1) /= a%ia(i2)) exit +!!$ i2 = i2 + 1 +!!$ end do +!!$ do +!!$ if (i1-1 < 1) exit +!!$ if (a%ia(i1-1) /= a%ia(i1)) exit +!!$ i1 = i1 - 1 +!!$ end do +!!$ ilr = ir +!!$ else +!!$ i1 = 1 +!!$ i2 = 1 +!!$ end if +!!$ nc = i2-i1+1 +!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) +!!$ 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 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(psbn_dupl_ovwrt_,psbn_dupl_err_) +!!$ ! Overwrite. +!!$ ! Cannot test for error, should have been caught earlier. +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ +!!$ if (ir /= ilr) then +!!$ i1 = psb_ibsrch(ir,nnz,a%ia) +!!$ i2 = i1 +!!$ do +!!$ if (i2+1 > nnz) exit +!!$ if (a%ia(i2+1) /= a%ia(i2)) exit +!!$ i2 = i2 + 1 +!!$ end do +!!$ do +!!$ if (i1-1 < 1) exit +!!$ if (a%ia(i1-1) /= a%ia(i1)) exit +!!$ i1 = i1 - 1 +!!$ end do +!!$ ilr = ir +!!$ else +!!$ i1 = 1 +!!$ i2 = 1 +!!$ end if +!!$ nc = i2-i1+1 +!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) +!!$ if (ip>0) then +!!$ a%val(i1+ip-1) = val(i) +!!$ else +!!$ info = i +!!$ return +!!$ end if +!!$ end if +!!$ end do +!!$ +!!$ case(psbn_dupl_add_) +!!$ ! Add +!!$ do i=1, nz +!!$ ir = ia(i) +!!$ ic = ja(i) +!!$ if ((ir > 0).and.(ir <= a%m)) then +!!$ +!!$ if (ir /= ilr) then +!!$ i1 = psb_ibsrch(ir,nnz,a%ia) +!!$ i2 = i1 +!!$ do +!!$ if (i2+1 > nnz) exit +!!$ if (a%ia(i2+1) /= a%ia(i2)) exit +!!$ i2 = i2 + 1 +!!$ end do +!!$ do +!!$ if (i1-1 < 1) exit +!!$ if (a%ia(i1-1) /= a%ia(i1)) exit +!!$ i1 = i1 - 1 +!!$ end do +!!$ ilr = ir +!!$ else +!!$ i1 = 1 +!!$ i2 = 1 +!!$ end if +!!$ nc = i2-i1+1 +!!$ ip = psb_issrch(ic,nc,a%ja(i1:i2)) +!!$ if (ip>0) then +!!$ a%val(i1+ip-1) = a%val(i1+ip-1) + val(i) +!!$ else +!!$ info = i +!!$ return +!!$ 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 +!!$ +!!$ end if +!!$ +!!$ end subroutine d_coo_srch_upd +!!$ +!!$#else + subroutine d_coo_srch_upd(nz,ia,ja,val,& + & aia,aja,aval,dupl,nza,nra,& + & info,gtl) + + use psb_error_mod + use psb_sort_mod + implicit none + + integer, intent(inout) :: aia(:),aja(:) + real(psb_dpk_), intent(inout) :: aval(:) + integer, intent(in) :: nz, dupl,nza, nra + integer, intent(in) :: ia(:),ja(:) + real(psb_dpk_), intent(in) :: val(:) + integer, intent(out) :: info + integer, intent(in), optional :: gtl(:) + integer :: i,ir,ic, ilr, ilc, ip, & + & i1,i2,nc,ng + integer :: debug_level, debug_unit + character(len=20) :: name='d_coo_srch_upd' + + info = 0 + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + + ilr = -1 + ilc = -1 + + + if (present(gtl)) then + ng = size(gtl) + + select case(dupl) + + case(psbn_dupl_ovwrt_,psbn_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + 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) + if ((ir > 0).and.(ir <= nra)) then + ic = gtl(ic) + if (ir /= ilr) then + i1 = psb_ibsrch(ir,nza,aia) + i2 = i1 + do + if (i2+1 > nza) exit + if (aia(i2+1) /= aia(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (aia(i1-1) /= aia(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + else + i1 = 1 + i2 = 1 + end if + nc = i2-i1+1 + ip = psb_issrch(ic,nc,aja(i1:i2)) + if (ip>0) then + aval(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.' + endif + end if + end do + case(psbn_dupl_add_) + ! Add + 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 <= nra)) then + + if (ir /= ilr) then + i1 = psb_ibsrch(ir,nza,aia) + i2 = i1 + do + if (i2+1 > nza) exit + if (aia(i2+1) /= aia(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (aia(i1-1) /= aia(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + else + i1 = 1 + i2 = 1 + end if + nc = i2-i1+1 + ip = psb_issrch(ic,nc,aja(i1:i2)) + if (ip>0) then + aval(i1+ip-1) = aval(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 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(psbn_dupl_ovwrt_,psbn_dupl_err_) + ! Overwrite. + ! Cannot test for error, should have been caught earlier. + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir > 0).and.(ir <= nra)) then + + if (ir /= ilr) then + i1 = psb_ibsrch(ir,nza,aia) + i2 = i1 + do + if (i2+1 > nza) exit + if (aia(i2+1) /= aia(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (aia(i1-1) /= aia(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + else + i1 = 1 + i2 = 1 + end if + nc = i2-i1+1 + ip = psb_issrch(ic,nc,aja(i1:i2)) + if (ip>0) then + aval(i1+ip-1) = val(i) + else + info = i + return + end if + end if + end do + + case(psbn_dupl_add_) + ! Add + do i=1, nz + ir = ia(i) + ic = ja(i) + if ((ir > 0).and.(ir <= nra)) then + + if (ir /= ilr) then + i1 = psb_ibsrch(ir,nza,aia) + i2 = i1 + do + if (i2+1 > nza) exit + if (aia(i2+1) /= aia(i2)) exit + i2 = i2 + 1 + end do + do + if (i1-1 < 1) exit + if (aia(i1-1) /= aia(i1)) exit + i1 = i1 - 1 + end do + ilr = ir + else + i1 = 1 + i2 = 1 + end if + nc = i2-i1+1 + ip = psb_issrch(ic,nc,aja(i1:i2)) + if (ip>0) then + aval(i1+ip-1) = aval(i1+ip-1) + val(i) + else + info = i + return + 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 + + end if + + end subroutine d_coo_srch_upd +!!$#endif +end subroutine d_coo_csins_impl + + +subroutine d_coo_to_coo_impl(a,b,info) + use psb_error_mod + use psb_realloc_mod + use psbn_d_base_mat_mod, psb_protect_name => d_coo_to_coo_impl + class(psbn_d_coo_sparse_mat), intent(in) :: a + class(psbn_d_coo_sparse_mat), intent(out) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='to_coo' + logical, parameter :: debug=.false. + + + call psb_erractionsave(err_act) + info = 0 + call b%set_nzeros(a%get_nzeros()) + call b%set_nrows(a%get_nrows()) + call b%set_ncols(a%get_ncols()) + call b%set_dupl(a%get_dupl()) + call b%set_state(a%get_state()) + call b%set_triangle(a%is_triangle()) + call b%set_upper(a%is_upper()) + + call b%reallocate(a%get_nzeros()) + + b%ia(:) = a%ia(:) + b%ja(:) = a%ja(:) + b%val(:) = a%val(:) + + call b%fix(info) + + if (info /= 0) 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 d_coo_to_coo_impl + +subroutine d_coo_from_coo_impl(a,b,info) + use psb_error_mod + use psb_realloc_mod + use psbn_d_base_mat_mod, psb_protect_name => d_coo_from_coo_impl + class(psbn_d_coo_sparse_mat), intent(inout) :: a + class(psbn_d_coo_sparse_mat), intent(in) :: b + integer, intent(out) :: info + + Integer :: err_act + character(len=20) :: name='from_coo' + logical, parameter :: debug=.false. + integer :: m,n,nz + + + call psb_erractionsave(err_act) + info = 0 + call a%set_nzeros(b%get_nzeros()) + call a%set_nrows(b%get_nrows()) + call a%set_ncols(b%get_ncols()) + call a%set_dupl(b%get_dupl()) + call a%set_state(b%get_state()) + call a%set_triangle(b%is_triangle()) + call a%set_upper(b%is_upper()) + + call a%reallocate(b%get_nzeros()) + + a%ia(:) = b%ia(:) + a%ja(:) = b%ja(:) + a%val(:) = b%val(:) + + call a%fix(info) + + if (info /= 0) 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 d_coo_from_coo_impl + + +subroutine d_fix_coo_impl(a,info,idir) + use psb_const_mod + use psb_error_mod + use psb_realloc_mod + use psbn_d_base_mat_mod, psb_protect_name => d_fix_coo_impl + use psb_string_mod + use psb_ip_reord_mod + + class(psbn_d_coo_sparse_mat), intent(inout) :: a + integer, intent(out) :: info + integer, intent(in), optional :: idir + integer, allocatable :: iaux(:) + !locals + Integer :: nza, nzl,iret,idir_, dupl_ + integer :: i,j, irw, icl, err_act + integer :: debug_level, debug_unit + character(len=20) :: name = 'psb_fixcoo' + + info = 0 + + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + if(debug_level >= psb_debug_serial_) & + & write(debug_unit,*) trim(name),': start ',& + & size(a%ia),size(a%ja) + if (present(idir)) then + idir_ = idir + else + idir_ = 0 + endif + + nza = a%get_nzeros() + if (nza < 2) return + dupl_ = a%get_dupl() + + allocate(iaux(nza+2),stat=info) + if (info /= 0) return + + + select case(idir_) + + case(0) ! Row major order + + call msort_up(nza,a%ia(1),iaux(1),iret) + if (iret == 0) & + & call psb_ip_reord(nza,a%val,a%ia,a%ja,iaux) + i = 1 + j = i + do while (i <= nza) + do while ((a%ia(j) == a%ia(i))) + j = j+1 + if (j > nza) exit + enddo + nzl = j - i + call msort_up(nzl,a%ja(i),iaux(1),iret) + if (iret == 0) & + & call psb_ip_reord(nzl,a%val(i:i+nzl-1),& + & a%ia(i:i+nzl-1),a%ja(i:i+nzl-1),iaux) + i = j + enddo + + i = 1 + irw = a%ia(i) + icl = a%ja(i) + j = 1 + + select case(dupl_) + case(psbn_dupl_ovwrt_) + + do + j = j + 1 + if (j > nza) exit + if ((a%ia(j) == irw).and.(a%ja(j) == icl)) then + a%val(i) = a%val(j) + else + i = i+1 + a%val(i) = a%val(j) + a%ia(i) = a%ia(j) + a%ja(i) = a%ja(j) + irw = a%ia(i) + icl = a%ja(i) + endif + enddo + + case(psbn_dupl_add_) + + do + j = j + 1 + if (j > nza) exit + if ((a%ia(j) == irw).and.(a%ja(j) == icl)) then + a%val(i) = a%val(i) + a%val(j) + else + i = i+1 + a%val(i) = a%val(j) + a%ia(i) = a%ia(j) + a%ja(i) = a%ja(j) + irw = a%ia(i) + icl = a%ja(i) + endif + enddo + + case(psbn_dupl_err_) + do + j = j + 1 + if (j > nza) exit + if ((a%ia(j) == irw).and.(a%ja(j) == icl)) then + call psb_errpush(130,name) + goto 9999 + else + i = i+1 + a%val(i) = a%val(j) + a%ia(i) = a%ia(j) + a%ja(i) = a%ja(j) + irw = a%ia(i) + icl = a%ja(i) + endif + enddo + + end select + + + if(debug_level >= psb_debug_serial_)& + & write(debug_unit,*) trim(name),': end second loop' + + case(1) ! Col major order + + call msort_up(nza,a%ja(1),iaux(1),iret) + if (iret == 0) & + & call psb_ip_reord(nza,a%val,a%ia,a%ja,iaux) + i = 1 + j = i + do while (i <= nza) + do while ((a%ja(j) == a%ja(i))) + j = j+1 + if (j > nza) exit + enddo + nzl = j - i + call msort_up(nzl,a%ia(i),iaux(1),iret) + if (iret == 0) & + & call psb_ip_reord(nzl,a%val(i:i+nzl-1),& + & a%ia(i:i+nzl-1),a%ja(i:i+nzl-1),iaux) + i = j + enddo + + i = 1 + irw = a%ia(i) + icl = a%ja(i) + j = 1 + + + select case(dupl_) + case(psbn_dupl_ovwrt_) + do + j = j + 1 + if (j > nza) exit + if ((a%ia(j) == irw).and.(a%ja(j) == icl)) then + a%val(i) = a%val(j) + else + i = i+1 + a%val(i) = a%val(j) + a%ia(i) = a%ia(j) + a%ja(i) = a%ja(j) + irw = a%ia(i) + icl = a%ja(i) + endif + enddo + + case(psbn_dupl_add_) + do + j = j + 1 + if (j > nza) exit + if ((a%ia(j) == irw).and.(a%ja(j) == icl)) then + a%val(i) = a%val(i) + a%val(j) + else + i = i+1 + a%val(i) = a%val(j) + a%ia(i) = a%ia(j) + a%ja(i) = a%ja(j) + irw = a%ia(i) + icl = a%ja(i) + endif + enddo + + case(psbn_dupl_err_) + do + j = j + 1 + if (j > nza) exit + if ((a%ia(j) == irw).and.(a%ja(j) == icl)) then + call psb_errpush(130,name) + goto 9999 + else + i = i+1 + a%val(i) = a%val(j) + a%ia(i) = a%ia(j) + a%ja(i) = a%ja(j) + irw = a%ia(i) + icl = a%ja(i) + endif + enddo + end select + if (debug_level >= psb_debug_serial_)& + & write(debug_unit,*) trim(name),': end second loop' + case default + write(debug_unit,*) trim(name),': unknown direction ',idir_ + end select + + call a%set_sorted() + call a%set_nzeros(i) + call a%set_asb() + + deallocate(iaux) + + 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 d_fix_coo_impl diff --git a/base/newserial/psbn_mat_mod.f03 b/base/newserial/psbn_mat_mod.f03 index 93657305..67e2f401 100644 --- a/base/newserial/psbn_mat_mod.f03 +++ b/base/newserial/psbn_mat_mod.f03 @@ -9,6 +9,30 @@ module psbn_d_mat_mod contains + procedure, pass(a) :: get_nrows + procedure, pass(a) :: get_ncols + procedure, pass(a) :: get_nzeros + procedure, pass(a) :: get_size + procedure, pass(a) :: get_state + + procedure, pass(a) :: get_dupl + procedure, pass(a) :: is_null + procedure, pass(a) :: is_bld + procedure, pass(a) :: is_upd + procedure, pass(a) :: is_asb + procedure, pass(a) :: is_sorted + procedure, pass(a) :: is_upper + procedure, pass(a) :: is_lower + procedure, pass(a) :: is_triangle + procedure, pass(a) :: is_unit + procedure, pass(a) :: get_neigh + procedure, pass(a) :: allocate_mn + procedure, pass(a) :: allocate_mnnz + procedure, pass(a) :: reallocate_nz + procedure, pass(a) :: free + generic, public :: allocate => allocate_mn, allocate_mnnz + generic, public :: reallocate => reallocate_nz + procedure, pass(a) :: d_csmv procedure, pass(a) :: d_csmm generic, public :: psbn_csmm => d_csmm, d_csmv @@ -21,6 +45,351 @@ module psbn_d_mat_mod contains + function get_dupl(a) result(res) + use psb_error_mod + class(psbn_d_sparse_mat), intent(in) :: a + integer :: res + + if (allocated(a%a)) then + res = a%a%get_dupl + else + res = psbn_invalid_ + end if + end function get_dupl + + + function get_state(a) result(res) + class(psbn_d_sparse_mat), intent(in) :: a + integer :: res + + if (allocated(a%a)) then + res = a%a%get_state() + else + res = psbn_spmat_null_ + end if + end function get_state + + function get_nrows(a) result(res) + class(psbn_d_sparse_mat), intent(in) :: a + integer :: res + + if (allocated(a%a)) then + res = a%a%get_nrows() + else + res = 0 + end if + + end function get_nrows + + function get_ncols(a) result(res) + class(psbn_d_sparse_mat), intent(in) :: a + integer :: res + + if (allocated(a%a)) then + res = a%a%get_ncols() + else + res = 0 + end if + + end function get_ncols + + function is_triangle(a) result(res) + class(psbn_d_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_triangle() + else + res = .false. + end if + + end function is_triangle + + function is_unit(a) result(res) + class(psbn_d_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_unit() + else + res = .false. + end if + + end function is_unit + + function is_upper(a) result(res) + class(psbn_d_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_upper() + else + res = .false. + end if + + end function is_upper + + function is_lower(a) result(res) + class(psbn_d_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = .not. a%a%is_upper() + else + res = .false. + end if + + end function is_lower + + function is_null(a) result(res) + class(psbn_d_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_null() + else + res = .true. + end if + + end function is_null + + function is_bld(a) result(res) + class(psbn_d_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_bld() + else + res = .false. + end if + + end function is_bld + + function is_upd(a) result(res) + class(psbn_d_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_upd() + else + res = .false. + end if + + end function is_upd + + function is_asb(a) result(res) + class(psbn_d_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_asb() + else + res = .false. + end if + + end function is_asb + + function is_sorted(a) result(res) + class(psbn_d_sparse_mat), intent(in) :: a + logical :: res + + if (allocated(a%a)) then + res = a%a%is_sorted() + else + res = .false. + end if + + end function is_sorted + + + function get_nzeros(a) result(res) + use psb_error_mod + class(psbn_d_sparse_mat), intent(in) :: a + integer :: res + + Integer :: err_act + character(len=20) :: name='get_nzeros' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + res = a%a%get_nzeros() + + 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 + + end function get_nzeros + + function get_size(a) result(res) + use psb_error_mod + class(psbn_d_sparse_mat), intent(in) :: a + integer :: res + + Integer :: err_act + character(len=20) :: name='get_size' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + res = a%a%get_size() + + 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 function get_size + + + subroutine get_neigh(a,idx,neigh,n,info,lev) + use psb_error_mod + class(psbn_d_sparse_mat), intent(in) :: a + integer, intent(in) :: idx + integer, intent(out) :: n + integer, allocatable, intent(out) :: neigh(:) + integer, intent(out) :: info + integer, optional, intent(in) :: lev + + Integer :: err_act + character(len=20) :: name='get_neigh' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%get_neigh(idx,neigh,n,info,lev) + + if (info /= 0) 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 get_neigh + + subroutine allocate_mn(m,n,a) + use psb_error_mod + integer, intent(in) :: m,n + class(psbn_d_sparse_mat), intent(inout) :: a + + Integer :: err_act + character(len=20) :: name='allocate_mn' + logical, parameter :: debug=.false. + + + call psb_erractionsave(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + call psb_errpush(700,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine allocate_mn + + subroutine allocate_mnnz(m,n,nz,a) + use psb_error_mod + integer, intent(in) :: m,n,nz + class(psbn_d_sparse_mat), intent(inout) :: a + Integer :: err_act + character(len=20) :: name='allocate_mnz' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + call psb_errpush(700,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine allocate_mnnz + + subroutine reallocate_nz(nz,a) + use psb_error_mod + integer, intent(in) :: nz + class(psbn_d_sparse_mat), intent(inout) :: a + Integer :: err_act + character(len=20) :: name='reallocate_nz' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + ! This is the base version. If we get here + ! it means the derived class is incomplete, + ! so we throw an error. + call psb_errpush(700,name) + + if (err_act /= psb_act_ret_) then + call psb_error() + end if + return + + end subroutine reallocate_nz + + subroutine free(a) + use psb_error_mod + class(psbn_d_sparse_mat), intent(inout) :: a + Integer :: err_act + character(len=20) :: name='free' + logical, parameter :: debug=.false. + + call psb_erractionsave(err_act) + if (.not.allocated(a%a)) then + info = 1121 + call psb_errpush(info,name) + goto 9999 + endif + + call a%a%free() + + return + + end subroutine free + + subroutine d_csmm(alpha,a,x,beta,y,info,trans) use psb_error_mod class(psbn_d_sparse_mat), intent(in) :: a diff --git a/base/serial/Makefile b/base/serial/Makefile index 7cf0b11a..470472c3 100644 --- a/base/serial/Makefile +++ b/base/serial/Makefile @@ -26,7 +26,7 @@ FOBJS = psb_cest.o psb_dcoins.o psb_dcsmm.o psb_dcsmv.o \ psb_cgelp.o psb_cspgtdiag.o psb_cspgtblk.o psb_cspgetrow.o\ psb_ccssm.o psb_ccssv.o psb_ccsmm.o psb_ccsmv.o psb_ctransp.o psb_ctransc.o\ psb_cspclip.o psb_crwextd.o psb_cspscal.o\ - psb_cnumbmm.o psb_csymbmm.o psb_cneigh.o psb_ip_reord_mod.o + psb_cnumbmm.o psb_csymbmm.o psb_cneigh.o # psb_dcsrp.o psb_zcsrp.o\ LIBDIR=.. @@ -43,7 +43,7 @@ lib1: $(FOBJS) psb_scoins.o psb_dcoins.o psb_zcoins.o: psb_update_mod.o psb_sspgetrow.o psb_dspgetrow.o psb_zspgetrow.o: psb_getrow_mod.o psb_sspcnv.o psb_dspcnv.o pzb_zspcnv.o: psb_regen_mod.o -psb_sfixcoo.o psb_dfixcoo.o psb_cfixcoo.o psb_zfixcoo.o: psb_ip_reord_mod.o +#psb_sfixcoo.o psb_dfixcoo.o psb_cfixcoo.o psb_zfixcoo.o: psb_ip_reord_mod.o auxd: psb_ip_reord_mod.o (cd aux; make lib) diff --git a/config/pac.m4 b/config/pac.m4 index 7a34c67e..146258cb 100644 --- a/config/pac.m4 +++ b/config/pac.m4 @@ -300,19 +300,20 @@ dnl AC_DEFUN([PAC_ARG_SERIAL_MPI], [ AC_MSG_CHECKING([whether we want serial (fake) mpi]) -AC_ARG_WITH(serial-mpi, -AC_HELP_STRING([--with-serial-mpi], -[Specify whether to enable a fake mpi library to run in serial mode. - --with-serial-mpi={yes|no}]), +AC_ARG_ENABLE(serial, +AC_HELP_STRING([--enable-serial], +[Specify whether to enable a fake mpi library to run in serial mode. ]), [ -pac_cv_serial_mpi="${withval}"; -], -[pac_cv_serial_mpi="no";] +pac_cv_serial_mpi="yes"; +] +dnl , +dnl [pac_cv_serial_mpi="no";] ) if test x"$pac_cv_serial_mpi" == x"yes" ; then AC_MSG_RESULT([yes.]) else - AC_MSG_RESULT([no.]) + pac_cv_serial_mpi="no"; + AC_MSG_RESULT([no.]) fi ] )