From f6c145f982ef4256ec0a7768487d50834ed68c1e Mon Sep 17 00:00:00 2001 From: Cirdans-Home Date: Mon, 16 Nov 2020 19:16:21 +0100 Subject: [PATCH] csall changed in allocate with CSR format --- prec/impl/psb_c_ilu0_fact.f90 | 154 ++++++++++---------- prec/impl/psb_c_iluk_fact.f90 | 214 +++++++++++++-------------- prec/impl/psb_c_ilut_fact.f90 | 264 +++++++++++++++++----------------- prec/impl/psb_d_ilu0_fact.f90 | 154 ++++++++++---------- prec/impl/psb_d_iluk_fact.f90 | 214 +++++++++++++-------------- prec/impl/psb_d_ilut_fact.f90 | 264 +++++++++++++++++----------------- prec/impl/psb_s_ilu0_fact.f90 | 154 ++++++++++---------- prec/impl/psb_s_iluk_fact.f90 | 214 +++++++++++++-------------- prec/impl/psb_s_ilut_fact.f90 | 264 +++++++++++++++++----------------- prec/impl/psb_z_ilu0_fact.f90 | 154 ++++++++++---------- prec/impl/psb_z_iluk_fact.f90 | 214 +++++++++++++-------------- prec/impl/psb_z_ilut_fact.f90 | 264 +++++++++++++++++----------------- 12 files changed, 1264 insertions(+), 1264 deletions(-) diff --git a/prec/impl/psb_c_ilu0_fact.f90 b/prec/impl/psb_c_ilu0_fact.f90 index c4097dea..1a3e1046 100644 --- a/prec/impl/psb_c_ilu0_fact.f90 +++ b/prec/impl/psb_c_ilu0_fact.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,21 +27,21 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! +! ! Moved here from MLD2P4, original copyright below. -! -! -! +! +! +! ! MLD2P4 version 2.2 ! MultiLevel Domain Decomposition Parallel Preconditioners Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2008-2018 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Daniela di Serafino -! +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -53,7 +53,7 @@ ! 3. The name of the MLD2P4 group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -65,8 +65,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: psb_cilu0_fact.f90 ! ! Subroutine: psb_cilu0_fact @@ -75,7 +75,7 @@ ! ! This routine computes either the ILU(0) or the MILU(0) factorization of ! the diagonal blocks of a distributed matrix. These factorizations are used -! to build the 'base preconditioner' (block-Jacobi preconditioner/solver, +! to build the 'base preconditioner' (block-Jacobi preconditioner/solver, ! Additive Schwarz preconditioner) corresponding to a given level of a ! multilevel preconditioner. ! @@ -83,10 +83,10 @@ ! Y. Saad, Iterative Methods for Sparse Linear Systems, Second Edition, ! SIAM, 2003, Chapter 10. ! -! The local matrix is stored into a and blck, as specified in the description -! of the arguments below. The storage format for both the L and U factors is CSR. -! The diagonal of the U factor is stored separately (actually, the inverse of the -! diagonal entries is stored; this is then managed in the solve stage associated +! The local matrix is stored into a and blck, as specified in the description +! of the arguments below. The storage format for both the L and U factors is CSR. +! The diagonal of the U factor is stored separately (actually, the inverse of the +! diagonal entries is stored; this is then managed in the solve stage associated ! to the ILU(0)/MILU(0) factorization). ! ! The routine copies and factors "on the fly" from a and blck into l (L factor), @@ -94,7 +94,7 @@ ! ! This implementation of ILU(0)/MILU(0) is faster than the implementation in ! psb_ziluk_fct (the latter routine performs the more general ILU(k)/MILU(k)). -! +! ! ! Arguments: ! ialg - integer, input. @@ -121,7 +121,7 @@ ! factorization. ! Note: its allocation is managed by the calling routine psb_ilu_bld, ! hence it cannot be only intent(out). -! info - integer, output. +! info - integer, output. ! Error code. ! blck - type(psb_cspmat_type), input, optional, target. ! The sparse matrix structure containing the remote rows of the @@ -129,7 +129,7 @@ ! to build an Additive Schwarz base preconditioner with overlap ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (see psb_fact_bld), then blck is empty. -! +! subroutine psb_cilu0_fact(ialg,a,l,u,d,info,blck, upd) use psb_base_mod @@ -157,30 +157,30 @@ subroutine psb_cilu0_fact(ialg,a,l,u,d,info,blck, upd) info = psb_success_ call psb_erractionsave(err_act) - ! + ! ! Point to / allocate memory for the incomplete factorization ! - if (present(blck)) then + if (present(blck)) then blck_ => blck else - allocate(blck_,stat=info) - if (info == psb_success_) call blck_%csall(izero,izero,info,ione) + allocate(blck_,stat=info) + if (info == psb_success_) call blck_%allocate(izero,izero,info,ione,type='CSR') if (info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='csall' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if endif - if (present(upd)) then - upd_ = psb_toupper(upd) + if (present(upd)) then + upd_ = psb_toupper(upd) else upd_ = 'F' end if - + m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& - & (m > size(d)) ) then + & (m > size(d)) ) then write(0,*) 'Wrong allocation status for L,D,U? ',& & l%get_nrows(),size(d),u%get_nrows() info = -1 @@ -212,12 +212,12 @@ subroutine psb_cilu0_fact(ialg,a,l,u,d,info,blck, upd) call u%set_triangle() call u%set_unit() call u%set_upper() - + ! ! Nullify pointer / deallocate memory ! - if (present(blck)) then - blck_ => null() + if (present(blck)) then + blck_ => null() else call blck_%free() if(info.ne.0) then @@ -226,7 +226,7 @@ subroutine psb_cilu0_fact(ialg,a,l,u,d,info,blck, upd) call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - deallocate(blck_) + deallocate(blck_) endif call psb_erractionrestore(err_act) @@ -243,9 +243,9 @@ contains ! Version: complex ! Note: internal subroutine of psb_cilu0_fact. ! - ! This routine computes either the ILU(0) or the MILU(0) factorization of the - ! diagonal blocks of a distributed matrix. - ! These factorizations are used to build the 'base preconditioner' + ! This routine computes either the ILU(0) or the MILU(0) factorization of the + ! diagonal blocks of a distributed matrix. + ! These factorizations are used to build the 'base preconditioner' ! (block-Jacobi preconditioner/solver, Additive Schwarz ! preconditioner) corresponding to a given level of a multilevel preconditioner. ! @@ -257,7 +257,7 @@ contains ! ! The routine copies and factors "on the fly" from the sparse matrix structures a ! and b into the arrays lval, uval, d (L, U without its diagonal, diagonal of U). - ! + ! ! ! Arguments: ! ialg - integer, input. @@ -296,7 +296,7 @@ contains ! according to the CSR storage format. ! lirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the L factor in lval, according to the CSR storage format. + ! of the L factor in lval, according to the CSR storage format. ! uval - complex(psb_spk_), dimension(:), input/output. ! The U factor in the incomplete factorization. ! The entries of U are stored according to the CSR format. @@ -305,18 +305,18 @@ contains ! according to the CSR storage format. ! uirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the U factor in uval, according to the CSR storage format. + ! of the U factor in uval, according to the CSR storage format. ! l1 - integer, output. ! The number of nonzero entries in lval. ! l2 - integer, output. ! The number of nonzero entries in uval. - ! info - integer, output. + ! info - integer, output. ! Error code. ! subroutine psb_cilu0_factint(ialg,a,b,& & d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,info) - implicit none + implicit none ! Arguments integer(psb_ipk_), intent(in) :: ialg @@ -332,7 +332,7 @@ contains complex(psb_spk_) :: dia,temp integer(psb_ipk_), parameter :: nrb=16 type(psb_c_coo_sparse_mat) :: trw - integer(psb_ipk_) :: int_err(5) + integer(psb_ipk_) :: int_err(5) character(len=20) :: name, ch_err name='psb_cilu0_factint' @@ -346,10 +346,10 @@ contains select case(ialg) case(psb_ilu_n_,psb_milu_n_) - ! Ok + ! Ok case default info=psb_err_input_asize_invalid_i_ - call psb_errpush(info,name,& + call psb_errpush(info,name,& & i_err=(/ione,ialg,izero,izero,izero/)) goto 9999 end select @@ -363,7 +363,7 @@ contains end if m = ma+mb - if (psb_toupper(upd) == 'F' ) then + if (psb_toupper(upd) == 'F' ) then lirp(1) = 1 uirp(1) = 1 l1 = 0 @@ -379,14 +379,14 @@ contains if (i <= ma) then ! ! Copy the i-th local row of the matrix, stored in a, - ! into lval/d(i)/uval + ! into lval/d(i)/uval ! call ilu_copyin(i,ma,a,i,ione,m,l1,lja,lval,& & d(i),l2,uja,uval,ktrw,trw,upd) else ! ! Copy the i-th local row of the matrix, stored in b - ! (as (i-ma)-th row), into lval/d(i)/uval + ! (as (i-ma)-th row), into lval/d(i)/uval ! call ilu_copyin(i-ma,mb,b,i,ione,m,l1,lja,lval,& & d(i),l2,uja,uval,ktrw,trw,upd) @@ -437,7 +437,7 @@ contains ! dia = dia - temp*uval(jj) cycle updateloop - ! + ! else if (j > i) then ! ! search u(i,*) (i-th row of U) for a matching index j @@ -454,23 +454,23 @@ contains end if enddo end if - ! - ! If we get here we missed the cycle updateloop, which means + ! + ! If we get here we missed the cycle updateloop, which means ! that this entry does not match; thus we accumulate on the ! diagonal for MILU(0). ! - if (ialg == psb_milu_n_) then + if (ialg == psb_milu_n_) then dia = dia - temp*uval(jj) end if enddo updateloop enddo - ! + ! ! Check the pivot size - ! + ! if (abs(dia) < s_epstol) then ! ! Too small pivot: unstable factorization - ! + ! info = psb_err_pivot_too_small_ int_err(1) = i write(ch_err,'(g20.10)') abs(dia) @@ -493,7 +493,7 @@ contains else write(0,*) 'Update not implemented ' info = 31 - call psb_errpush(info,name,& + call psb_errpush(info,name,& & i_err=(/ione*13,izero,izero,izero,izero/),a_err=upd) goto 9999 @@ -515,7 +515,7 @@ contains ! Version: complex ! Note: internal subroutine of psb_cilu0_fact ! - ! This routine copies a row of a sparse matrix A, stored in the psb_cspmat_type + ! This routine copies a row of a sparse matrix A, stored in the psb_cspmat_type ! data structure a, into the arrays lval and uval and into the scalar variable ! dia, corresponding to the lower and upper triangles of A and to the diagonal ! entry of the row, respectively. The entries in lval and uval are stored @@ -529,14 +529,14 @@ contains ! ! The routine is used by psb_cilu0_factint in the computation of the ILU(0)/MILU(0) ! factorization of a local sparse matrix. - ! + ! ! TODO: modify the routine to allow copying into output L and U that are ! already filled with indices; this would allow computing an ILU(k) pattern, - ! then use the ILU(0) internal for subsequent calls with the same pattern. + ! then use the ILU(0) internal for subsequent calls with the same pattern. ! ! Arguments: ! i - integer, input. - ! The local index of the row to be extracted from the + ! The local index of the row to be extracted from the ! sparse matrix structure a. ! m - integer, input. ! The number of rows of the local matrix stored into a. @@ -570,13 +570,13 @@ contains ! to the CSR storage format. ! uval - complex(psb_spk_), dimension(:), input/output. ! The array where the entries of the row corresponding to the - ! upper triangle are copied. + ! upper triangle are copied. ! ktrw - integer, input/output. ! The index identifying the last entry taken from the ! staging buffer trw. See below. ! trw - type(psb_cspmat_type), input/output. ! A staging buffer. If the matrix A is not in CSR format, we use - ! the psb_sp_getblk routine and store its output in trw; when we + ! the psb_sp_getblk routine and store its output in trw; when we ! need to call psb_sp_getblk we do it for a block of rows, and then ! we consume them from trw in successive calls to this routine, ! until we empty the buffer. Thus we will make a call to psb_sp_getblk @@ -608,10 +608,10 @@ contains if (psb_errstatus_fatal()) then info = psb_err_internal_error_; goto 9999 end if - if (psb_toupper(upd) == 'F') then + if (psb_toupper(upd) == 'F') then - select type(aa => a%a) - type is (psb_c_csr_sparse_mat) + select type(aa => a%a) + type is (psb_c_csr_sparse_mat) ! ! Take a fast shortcut if the matrix is stored in CSR format @@ -633,19 +633,19 @@ contains end if enddo - class default + class default ! - ! Otherwise use psb_sp_getblk, slower but able (in principle) of + ! Otherwise use psb_sp_getblk, slower but able (in principle) of ! handling any format. In this case, a block of rows is extracted ! instead of a single row, for performance reasons, and these ! rows are copied one by one into lval, dia, uval, through ! successive calls to ilu_copyin. ! - if ((mod(i,nrb) == 1).or.(nrb == 1)) then + if ((mod(i,nrb) == 1).or.(nrb == 1)) then irb = min(m-i+1,nrb) - call aa%csget(i,i+irb-1,trw,info) + call aa%csget(i,i+irb-1,trw,info) if(info /= psb_success_) then info=psb_err_from_subroutine_ ch_err='csget' @@ -656,7 +656,7 @@ contains end if nz = trw%get_nzeros() - do + do if (ktrw > nz) exit if (trw%ia(ktrw) > i) exit k = trw%ja(ktrw) @@ -680,7 +680,7 @@ contains write(0,*) 'Update not implemented ' info = 31 - call psb_errpush(info,name,& + call psb_errpush(info,name,& & i_err=(/ione*13,izero,izero,izero,izero/),a_err=upd) goto 9999 diff --git a/prec/impl/psb_c_iluk_fact.f90 b/prec/impl/psb_c_iluk_fact.f90 index 8748816d..c4ebc678 100644 --- a/prec/impl/psb_c_iluk_fact.f90 +++ b/prec/impl/psb_c_iluk_fact.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,21 +27,21 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! +! ! Moved here from MLD2P4, original copyright below. -! -! -! +! +! +! ! MLD2P4 version 2.2 ! MultiLevel Domain Decomposition Parallel Preconditioners Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2008-2018 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Daniela di Serafino -! +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -53,7 +53,7 @@ ! 3. The name of the MLD2P4 group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -65,8 +65,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: psb_ciluk_fact.f90 ! ! Subroutine: psb_ciluk_fact @@ -74,7 +74,7 @@ ! Contains: psb_ciluk_factint, iluk_copyin, iluk_fact, iluk_copyout. ! ! This routine computes either the ILU(k) or the MILU(k) factorization of the -! diagonal blocks of a distributed matrix. These factorizations are used to +! diagonal blocks of a distributed matrix. These factorizations are used to ! build the 'base preconditioner' (block-Jacobi preconditioner/solver, ! Additive Schwarz preconditioner) corresponding to a certain level of a ! multilevel preconditioner. @@ -88,7 +88,7 @@ ! U factors is CSR. The diagonal of the U factor is stored separately (actually, ! the inverse of the diagonal entries is stored; this is then managed in the solve ! stage associated to the ILU(k)/MILU(k) factorization). -! +! ! ! Arguments: ! fill_in - integer, input. @@ -118,7 +118,7 @@ ! factorization. ! Note: its allocation is managed by the calling routine psb_ilu_bld, ! hence it cannot be only intent(out). -! info - integer, output. +! info - integer, output. ! Error code. ! blck - type(psb_cspmat_type), input, optional, target. ! The sparse matrix structure containing the remote rows of the @@ -126,7 +126,7 @@ ! to build an Additive Schwarz base preconditioner with overlap ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (see psb_fact_bld), then blck does not contain any row. -! +! subroutine psb_ciluk_fact(fill_in,ialg,a,l,u,d,info,blck) use psb_base_mod @@ -143,7 +143,7 @@ subroutine psb_ciluk_fact(fill_in,ialg,a,l,u,d,info,blck) complex(psb_spk_), intent(inout) :: d(:) ! Local Variables integer(psb_ipk_) :: l1, l2, m, err_act - + type(psb_cspmat_type), pointer :: blck_ type(psb_c_csr_sparse_mat) :: ll, uu character(len=20) :: name, ch_err @@ -152,17 +152,17 @@ subroutine psb_ciluk_fact(fill_in,ialg,a,l,u,d,info,blck) info = psb_success_ call psb_erractionsave(err_act) - ! + ! ! Point to / allocate memory for the incomplete factorization ! - if (present(blck)) then + if (present(blck)) then blck_ => blck else - allocate(blck_,stat=info) - if (info == psb_success_) call blck_%csall(izero,izero,info,ione) + allocate(blck_,stat=info) + if (info == psb_success_) call blck_%allocate(izero,izero,info,ione,type='CSR') if (info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='csall' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if @@ -170,7 +170,7 @@ subroutine psb_ciluk_fact(fill_in,ialg,a,l,u,d,info,blck) m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& - & (m > size(d)) ) then + & (m > size(d)) ) then write(0,*) 'Wrong allocation status for L,D,U? ',& & l%get_nrows(),size(d),u%get_nrows() info = -1 @@ -203,15 +203,15 @@ subroutine psb_ciluk_fact(fill_in,ialg,a,l,u,d,info,blck) call u%set_triangle() call u%set_unit() call u%set_upper() - + ! ! Nullify pointer / deallocate memory ! - if (present(blck)) then - blck_ => null() + if (present(blck)) then + blck_ => null() else call blck_%free() - deallocate(blck_,stat=info) + deallocate(blck_,stat=info) if(info.ne.0) then info=psb_err_from_subroutine_ ch_err='psb_sp_free' @@ -280,7 +280,7 @@ contains ! according to the CSR storage format. ! lia2 - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the L factor in laspk, according to the CSR storage format. + ! of the L factor in laspk, according to the CSR storage format. ! uval - complex(psb_spk_), dimension(:), input/output. ! The U factor in the incomplete factorization. ! The entries of U are stored according to the CSR format. @@ -289,12 +289,12 @@ contains ! according to the CSR storage format. ! uirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the U factor in uval, according to the CSR storage format. + ! of the U factor in uval, according to the CSR storage format. ! l1 - integer, output ! The number of nonzero entries in laspk. ! l2 - integer, output ! The number of nonzero entries in uval. - ! info - integer, output. + ! info - integer, output. ! Error code. ! subroutine psb_ciluk_factint(fill_in,ialg,a,b,& @@ -304,7 +304,7 @@ contains implicit none - ! Arguments + ! Arguments integer(psb_ipk_), intent(in) :: fill_in, ialg type(psb_cspmat_type),intent(in) :: a,b integer(psb_ipk_),intent(inout) :: l1,l2,info @@ -330,14 +330,14 @@ contains select case(ialg) case(psb_ilu_n_,psb_milu_n_) - ! Ok + ! Ok case default info=psb_err_input_asize_invalid_i_ call psb_errpush(info,name,& & i_err=(/itwo,ialg,izero,izero,izero/)) goto 9999 end select - if (fill_in < 0) then + if (fill_in < 0) then info=psb_err_input_asize_invalid_i_ call psb_errpush(info,name, & & i_err=(/ione,fill_in,izero,izero,izero/)) @@ -349,7 +349,7 @@ contains m = ma+mb ! - ! Allocate a temporary buffer for the iluk_copyin function + ! Allocate a temporary buffer for the iluk_copyin function ! call trw%allocate(izero,izero,ione) @@ -361,7 +361,7 @@ contains call psb_errpush(info,name,a_err='psb_sp_all') goto 9999 end if - + l1=0 l2=0 lirp(1) = 1 @@ -381,34 +381,34 @@ contains uplevs(:) = m+1 row(:) = czero rowlevs(:) = -(m+1) - + ! ! Cycle over the matrix rows ! do i = 1, m - + ! ! At each iteration of the loop we keep in a heap the column indices ! affected by the factorization. The heap is initialized and filled ! in the iluk_copyin routine, and updated during the elimination, in ! the iluk_fact routine. The heap is ideal because at each step we need ! the lowest index, but we also need to insert new items, and the heap - ! allows to do both in log time. + ! allows to do both in log time. ! d(i) = czero if (i<=ma) then ! - ! Copy into trw the i-th local row of the matrix, stored in a - ! + ! Copy into trw the i-th local row of the matrix, stored in a + ! call iluk_copyin(i,ma,a,ione,m,row,rowlevs,heap,ktrw,trw,info) else ! ! Copy into trw the i-th local row of the matrix, stored in b - ! (as (i-ma)-th row) - ! + ! (as (i-ma)-th row) + ! call iluk_copyin(i-ma,mb,b,ione,m,row,rowlevs,heap,ktrw,trw,info) endif - + ! Do an elimination step on the current row. It turns out we only ! need to keep track of fill levels for the upper triangle, hence we ! do not have a lowlevs variable. @@ -417,7 +417,7 @@ contains & d,uja,uirp,uval,uplevs,nidx,idxs,info) ! ! Copy the row into lval/d(i)/uval - ! + ! if (info == psb_success_) call iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,& & l1,l2,lja,lirp,lval,d,uja,uirp,uval,uplevs,info) if (info /= psb_success_) then @@ -474,11 +474,11 @@ contains ! ! This routine is used by psb_ciluk_factint in the computation of the ! ILU(k)/MILU(k) factorization of a local sparse matrix. - ! + ! ! ! Arguments: ! i - integer, input. - ! The local index of the row to be extracted from the + ! The local index of the row to be extracted from the ! sparse matrix structure a. ! m - integer, input. ! The number of rows of the local matrix stored into a. @@ -510,7 +510,7 @@ contains ! staging buffer trw. See below. ! trw - type(psb_cspmat_type), input/output. ! A staging buffer. If the matrix A is not in CSR format, we use - ! the psb_sp_getblk routine and store its output in trw; when we + ! the psb_sp_getblk routine and store its output in trw; when we ! need to call psb_sp_getblk we do it for a block of rows, and then ! we consume them from trw in successive calls to this routine, ! until we empty the buffer. Thus we will make a call to psb_sp_getblk @@ -521,8 +521,8 @@ contains use psb_base_mod implicit none - - ! Arguments + + ! Arguments type(psb_cspmat_type), intent(in) :: a type(psb_c_coo_sparse_mat), intent(inout) :: trw integer(psb_ipk_), intent(in) :: i,m,jmin,jmax @@ -542,17 +542,17 @@ contains if (psb_errstatus_fatal()) then info = psb_err_internal_error_; goto 9999 end if - call heap%init(info) + call heap%init(info) - select type (aa=> a%a) - type is (psb_c_csr_sparse_mat) + select type (aa=> a%a) + type is (psb_c_csr_sparse_mat) ! ! Take a fast shortcut if the matrix is stored in CSR format ! - + do j = aa%irp(i), aa%irp(i+1) - 1 k = aa%ja(j) - if ((jmin<=k).and.(k<=jmax)) then + if ((jmin<=k).and.(k<=jmax)) then row(k) = aa%val(j) rowlevs(k) = 0 call heap%insert(k,info) @@ -562,14 +562,14 @@ contains class default ! - ! Otherwise use psb_sp_getblk, slower but able (in principle) of + ! Otherwise use psb_sp_getblk, slower but able (in principle) of ! handling any format. In this case, a block of rows is extracted ! instead of a single row, for performance reasons, and these ! rows are copied one by one into the array row, through successive ! calls to iluk_copyin. ! - if ((mod(i,nrb) == 1).or.(nrb == 1)) then + if ((mod(i,nrb) == 1).or.(nrb == 1)) then irb = min(m-i+1,nrb) call aa%csget(i,i+irb-1,trw,info) if (info /= psb_success_) then @@ -581,11 +581,11 @@ contains ktrw=1 end if nz = trw%get_nzeros() - do + do if (ktrw > nz) exit if (trw%ia(ktrw) > i) exit k = trw%ja(ktrw) - if ((jmin<=k).and.(k<=jmax)) then + if ((jmin<=k).and.(k<=jmax)) then row(k) = trw%val(ktrw) rowlevs(k) = 0 call heap%insert(k,info) @@ -674,10 +674,10 @@ contains use psb_base_mod - implicit none + implicit none ! Arguments - type(psb_i_heap), intent(inout) :: heap + type(psb_i_heap), intent(inout) :: heap integer(psb_ipk_), intent(in) :: i, fill_in integer(psb_ipk_), intent(inout) :: nidx,info integer(psb_ipk_), intent(inout) :: rowlevs(:) @@ -690,7 +690,7 @@ contains complex(psb_spk_) :: rwk info = psb_success_ - if (.not.allocated(idxs)) then + if (.not.allocated(idxs)) then allocate(idxs(200),stat=info) if (info /= psb_success_) return endif @@ -702,34 +702,34 @@ contains ! do ! Beware: (iret < 0) means that the heap is empty, not an error. - call heap%get_first(k,iret) + call heap%get_first(k,iret) if (iret < 0) return - ! + ! ! Just in case an index has been put on the heap more than once. ! if (k == lastk) cycle - lastk = k + lastk = k nidx = nidx + 1 - if (nidx>size(idxs)) then + if (nidx>size(idxs)) then call psb_realloc(nidx+psb_heap_resize,idxs,info) if (info /= psb_success_) return end if idxs(nidx) = k - - if ((row(k) /= czero).and.(rowlevs(k) <= fill_in).and.(ki) then + else if (j>i) then ! ! Copy the upper part of the row - ! - if (rowlevs(j) <= fill_in) then - l2 = l2 + 1 - if (size(uval) < l2) then - ! + ! + if (rowlevs(j) <= fill_in) then + l2 = l2 + 1 + if (size(uval) < l2) then + ! ! Figure out a good reallocation size! ! isz = max((l2/i)*m,int(1.2*l2),l2+100) - call psb_realloc(isz,uval,info) - if (info == psb_success_) call psb_realloc(isz,uja,info) + call psb_realloc(isz,uval,info) + if (info == psb_success_) call psb_realloc(isz,uja,info) if (info == psb_success_) call psb_realloc(isz,uplevs,info,pad=(m+1)) - if (info /= psb_success_) then + if (info /= psb_success_) then info=psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') goto 9999 @@ -941,11 +941,11 @@ contains end if uja(l2) = j uval(l2) = row(j) - uplevs(l2) = rowlevs(j) + uplevs(l2) = rowlevs(j) else if (ialg == psb_milu_n_) then ! ! MILU(k): add discarded entries to the diagonal one - ! + ! d(i) = d(i) + row(j) end if ! @@ -963,13 +963,13 @@ contains lirp(i+1) = l1 + 1 uirp(i+1) = l2 + 1 - ! + ! ! Check the pivot size ! if (abs(d(i)) < s_epstol) then ! ! Too small pivot: unstable factorization - ! + ! info = psb_err_pivot_too_small_ int_err(1) = i write(ch_err,'(g20.10)') d(i) diff --git a/prec/impl/psb_c_ilut_fact.f90 b/prec/impl/psb_c_ilut_fact.f90 index 06b8b477..633899de 100644 --- a/prec/impl/psb_c_ilut_fact.f90 +++ b/prec/impl/psb_c_ilut_fact.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,21 +27,21 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! +! ! Moved here from MLD2P4, original copyright below. -! -! -! +! +! +! ! MLD2P4 version 2.2 ! MultiLevel Domain Decomposition Parallel Preconditioners Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2008-2018 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Daniela di Serafino -! +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -53,7 +53,7 @@ ! 3. The name of the MLD2P4 group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -65,8 +65,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: psb_cilut_fact.f90 ! ! Subroutine: psb_cilut_fact @@ -87,7 +87,7 @@ ! CSR. The diagonal of the U factor is stored separately (actually, the ! inverse of the diagonal entries is stored; this is then managed in the ! solve stage associated to the ILU(k,t) factorization). -! +! ! ! Arguments: ! fill_in - integer, input. @@ -114,7 +114,7 @@ ! factorization. ! Note: its allocation is managed by the calling routine psb_ilu_bld, ! hence it cannot be only intent(out). -! info - integer, output. +! info - integer, output. ! Error code. ! blck - type(psb_cspmat_type), input, optional, target. ! The sparse matrix structure containing the remote rows of the @@ -122,9 +122,9 @@ ! to build an Additive Schwarz base preconditioner with overlap ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (see psb_fact_bld), then blck does not contain any row. -! +! subroutine psb_cilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) - + use psb_base_mod use psb_c_ilu_fact_mod, psb_protect_name => psb_cilut_fact @@ -141,7 +141,7 @@ subroutine psb_cilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) integer(psb_ipk_), intent(in), optional :: iscale ! Local Variables integer(psb_ipk_) :: l1, l2, m, err_act, iscale_ - + type(psb_cspmat_type), pointer :: blck_ type(psb_c_csr_sparse_mat) :: ll, uu real(psb_spk_) :: scale @@ -151,34 +151,34 @@ subroutine psb_cilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) info = psb_success_ call psb_erractionsave(err_act) - if (fill_in < 0) then + if (fill_in < 0) then info=psb_err_input_asize_invalid_i_ call psb_errpush(info,name, & & i_err=(/ione,fill_in,izero,izero,izero/)) goto 9999 end if - ! + ! ! Point to / allocate memory for the incomplete factorization ! - if (present(blck)) then + if (present(blck)) then blck_ => blck else - allocate(blck_,stat=info) - if (info == psb_success_) call blck_%csall(izero,izero,info,ione) + allocate(blck_,stat=info) + if (info == psb_success_) call blck_%allocate(izero,izero,info,ione,type='CSR') if (info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='csall' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if endif - if (present(iscale)) then + if (present(iscale)) then iscale_ = iscale else iscale_ = psb_ilu_scale_none_ end if - - select case(iscale_) + + select case(iscale_) case(psb_ilu_scale_none_) scale = sone case(psb_ilu_scale_maxval_) @@ -189,10 +189,10 @@ subroutine psb_cilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) call psb_errpush(info,name,i_err=(/ione*9,iscale_,izero,izero,izero/)) goto 9999 end select - + m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& - & (m > size(d)) ) then + & (m > size(d)) ) then write(0,*) 'Wrong allocation status for L,D,U? ',& & l%get_nrows(),size(d),u%get_nrows() info = -1 @@ -225,15 +225,15 @@ subroutine psb_cilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) call u%set_triangle() call u%set_unit() call u%set_upper() - + ! ! Nullify pointer / deallocate memory ! - if (present(blck)) then - blck_ => null() + if (present(blck)) then + blck_ => null() else call blck_%free() - deallocate(blck_,stat=info) + deallocate(blck_,stat=info) if(info.ne.0) then info=psb_err_from_subroutine_ ch_err='psb_sp_free' @@ -298,7 +298,7 @@ contains ! according to the CSR storage format. ! lirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the L factor in lval, according to the CSR storage format. + ! of the L factor in lval, according to the CSR storage format. ! uval - complex(psb_spk_), dimension(:), input/output. ! The U factor in the incomplete factorization. ! The entries of U are stored according to the CSR format. @@ -307,12 +307,12 @@ contains ! according to the CSR storage format. ! uirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the U factor in uval, according to the CSR storage format. + ! of the U factor in uval, according to the CSR storage format. ! l1 - integer, output ! The number of nonzero entries in lval. ! l2 - integer, output ! The number of nonzero entries in uval. - ! info - integer, output. + ! info - integer, output. ! Error code. ! subroutine psb_cilut_factint(fill_in,thres,a,b,& @@ -320,7 +320,7 @@ contains use psb_base_mod - implicit none + implicit none ! Arguments integer(psb_ipk_), intent(in) :: fill_in @@ -355,7 +355,7 @@ contains m = ma+mb ! - ! Allocate a temporary buffer for the ilut_copyin function + ! Allocate a temporary buffer for the ilut_copyin function ! call trw%allocate(izero,izero,ione) if (info == psb_success_) call psb_ensure_size(m+1,lirp,info) @@ -366,7 +366,7 @@ contains call psb_errpush(info,name,a_err='psb_sp_all') goto 9999 end if - + l1=0 l2=0 lirp(1) = 1 @@ -396,10 +396,10 @@ contains ! in the ilut_copyin function, and updated during the elimination, in ! the ilut_fact routine. The heap is ideal because at each step we need ! the lowest index, but we also need to insert new items, and the heap - ! allows to do both in log time. + ! allows to do both in log time. ! d(i) = czero - if (i<=ma) then + if (i<=ma) then call ilut_copyin(i,ma,a,i,ione,m,nlw,nup,jmaxup,nrmi,weight,& & row,heap,ktrw,trw,info) else @@ -414,7 +414,7 @@ contains & d,uja,uirp,uval,nidx,idxs,info) ! ! Copy the row into lval/d(i)/uval - ! + ! if (info == psb_success_) call ilut_copyout(fill_in,thres,i,m,& & nlw,nup,jmaxup,nrmi,row,nidx,idxs,& & l1,l2,lja,lirp,lval,d,uja,uirp,uval,info) @@ -429,7 +429,7 @@ contains ! ! Adjust diagonal accounting for scale factor ! - if (weight /= sone) then + if (weight /= sone) then d(1:m) = d(1:m)*weight end if @@ -469,7 +469,7 @@ contains ! - storing into a heap the column indices of the nonzero entries of the copied ! row; ! - computing the column index of the first entry with maximum absolute value - ! in the part of the row belonging to the upper triangle; + ! in the part of the row belonging to the upper triangle; ! - computing the 2-norm of the row. ! The output array row is such that it contains a full row of A, i.e. it contains ! also the zero entries of the row. This is useful for the elimination step @@ -482,11 +482,11 @@ contains ! ! This routine is used by psb_cilut_factint in the computation of the ILU(k,t) ! factorization of a local sparse matrix. - ! + ! ! ! Arguments: ! i - integer, input. - ! The local index of the row to be extracted from the + ! The local index of the row to be extracted from the ! sparse matrix structure a. ! m - integer, input. ! The number of rows of the local matrix stored into a. @@ -533,7 +533,7 @@ contains ! staging buffer trw. See below. ! trw - type(psb_cspmat_type), input/output. ! A staging buffer. If the matrix A is not in CSR format, we use - ! the psb_sp_getblk routine and store its output in trw; when we + ! the psb_sp_getblk routine and store its output in trw; when we ! need to call psb_sp_getblk we do it for a block of rows, and then ! we consume them from trw in successive calls to this routine, ! until we empty the buffer. Thus we will make a call to psb_sp_getblk @@ -542,7 +542,7 @@ contains subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,& & nrmi,weight,row,heap,ktrw,trw,info) use psb_base_mod - implicit none + implicit none type(psb_cspmat_type), intent(in) :: a type(psb_c_coo_sparse_mat), intent(inout) :: trw integer(psb_ipk_), intent(in) :: i, m,jmin,jmax,jd @@ -584,29 +584,29 @@ contains dmaxup = szero nrmi = szero - select type (aa=> a%a) - type is (psb_c_csr_sparse_mat) + select type (aa=> a%a) + type is (psb_c_csr_sparse_mat) ! ! Take a fast shortcut if the matrix is stored in CSR format ! do j = aa%irp(i), aa%irp(i+1) - 1 k = aa%ja(j) - if ((jmin<=k).and.(k<=jmax)) then - row(k) = aa%val(j)*weight + if ((jmin<=k).and.(k<=jmax)) then + row(k) = aa%val(j)*weight call heap%insert(k,info) if (info /= psb_success_) exit - if (kjd) then + if (kjd) then nup = nup + 1 - if (abs(row(k))>dmaxup) then + if (abs(row(k))>dmaxup) then jmaxup = k dmaxup = abs(row(k)) end if end if end if end do - if (info /= psb_success_) then + if (info /= psb_success_) then info=psb_err_from_subroutine_ call psb_errpush(info,name,a_err='psb_insert_heap') goto 9999 @@ -617,16 +617,16 @@ contains class default - + ! - ! Otherwise use psb_sp_getblk, slower but able (in principle) of + ! Otherwise use psb_sp_getblk, slower but able (in principle) of ! handling any format. In this case, a block of rows is extracted ! instead of a single row, for performance reasons, and these ! rows are copied one by one into the array row, through successive ! calls to ilut_copyin. ! - if ((mod(i,nrb) == 1).or.(nrb == 1)) then + if ((mod(i,nrb) == 1).or.(nrb == 1)) then irb = min(m-i+1,nrb) call aa%csget(i,i+irb-1,trw,info) if (info /= psb_success_) then @@ -639,18 +639,18 @@ contains kin = ktrw nz = trw%get_nzeros() - do + do if (ktrw > nz) exit if (trw%ia(ktrw) > i) exit k = trw%ja(ktrw) - if ((jmin<=k).and.(k<=jmax)) then + if ((jmin<=k).and.(k<=jmax)) then row(k) = trw%val(ktrw)*weight call heap%insert(k,info) if (info /= psb_success_) exit - if (kjd) then + if (kjd) then nup = nup + 1 - if (abs(row(k))>dmaxup) then + if (abs(row(k))>dmaxup) then jmaxup = k dmaxup = abs(row(k)) end if @@ -734,10 +734,10 @@ contains use psb_base_mod - implicit none + implicit none ! Arguments - type(psb_i_heap), intent(inout) :: heap + type(psb_i_heap), intent(inout) :: heap integer(psb_ipk_), intent(in) :: i integer(psb_ipk_), intent(inout) :: nidx,info real(psb_spk_), intent(in) :: thres,nrmi @@ -753,23 +753,23 @@ contains call psb_ensure_size(200*ione,idxs,info) if (info /= psb_success_) return nidx = 0 - lastk = -1 + lastk = -1 ! ! Do while there are indices to be processed ! do - call heap%get_first(k,iret) + call heap%get_first(k,iret) if (iret < 0) exit - ! + ! ! An index may have been put on the heap more than once. ! if (k == lastk) cycle - lastk = k + lastk = k lowert: if (k nidx) exit if (idxs(idxp) >= i) exit @@ -987,8 +987,8 @@ contains ! if (abs(witem) < thres*nrmi) cycle - nz = nz + 1 - xw(nz) = witem + nz = nz + 1 + xw(nz) = witem xwid(nz) = widx call heap%insert(witem,widx,info) if (info /= psb_success_) then @@ -1000,9 +1000,9 @@ contains ! ! Now we have to take out the first nlw+fill_in entries - ! + ! if (nz <= nlw+fill_in) then - ! + ! ! Just copy everything from xw, and it is already ordered ! else @@ -1014,7 +1014,7 @@ contains call psb_errpush(info,name,a_err='psb_heap_get_first') goto 9999 end if - + xw(k) = witem xwid(k) = widx end do @@ -1029,15 +1029,15 @@ contains ! Copy out the lower part of the row ! do k=1,nz - l1 = l1 + 1 + l1 = l1 + 1 if (size(lval) < l1) then - ! + ! ! Figure out a good reallocation size! - ! + ! isz = (max((l1/i)*m,int(1.2*l1),l1+100)) - call psb_realloc(isz,lval,info) - if (info == psb_success_) call psb_realloc(isz,lja,info) - if (info /= psb_success_) then + call psb_realloc(isz,lval,info) + if (info == psb_success_) call psb_realloc(isz,lja,info) + if (info /= psb_success_) then info=psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') goto 9999 @@ -1046,25 +1046,25 @@ contains lja(l1) = xwid(k) lval(l1) = xw(indx(k)) end do - + ! ! Make sure idxp points to the diagonal entry ! - if (idxp <= size(idxs)) then - if (idxs(idxp) < i) then - do + if (idxp <= size(idxs)) then + if (idxs(idxp) < i) then + do idxp = idxp + 1 if (idxp > nidx) exit if (idxs(idxp) >= i) exit end do end if end if - if (idxp > size(idxs)) then + if (idxp > size(idxs)) then !!$ write(0,*) 'Warning: missing diagonal element in the row ' else - if (idxs(idxp) > i) then + if (idxs(idxp) > i) then !!$ write(0,*) 'Warning: missing diagonal element in the row ' - else if (idxs(idxp) /= i) then + else if (idxs(idxp) /= i) then !!$ write(0,*) 'Warning: impossible error: diagonal has vanished' else ! @@ -1076,7 +1076,7 @@ contains if (abs(d(i)) < s_epstol) then ! ! Too small pivot: unstable factorization - ! + ! info = psb_err_pivot_too_small_ int_err(1) = i write(ch_err,'(g20.10)') d(i) @@ -1092,7 +1092,7 @@ contains end if ! - ! Now the upper part + ! Now the upper part ! call heap%init(info,dir=psb_asort_down_) @@ -1104,15 +1104,15 @@ contains nz = 0 do - + idxp = idxp + 1 if (idxp > nidx) exit widx = idxs(idxp) - if (widx <= i) then + if (widx <= i) then !!$ write(0,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) cycle end if - if (widx > m) then + if (widx > m) then !!$ write(0,*) 'Warning: impossible value',widx,i,idxp,idxs(idxp) cycle end if @@ -1120,12 +1120,12 @@ contains ! ! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway. ! - if ((widx /= jmaxup) .and. (abs(witem) < thres*nrmi)) then - cycle + if ((widx /= jmaxup) .and. (abs(witem) < thres*nrmi)) then + cycle end if nz = nz + 1 - xw(nz) = witem + xw(nz) = witem xwid(nz) = widx call heap%insert(witem,widx,info) if (info /= psb_success_) then @@ -1133,7 +1133,7 @@ contains call psb_errpush(info,name,a_err='psb_insert_heap') goto 9999 end if - + end do ! @@ -1141,7 +1141,7 @@ contains ! we include entry jmaxup. ! if (nz <= nup+fill_in) then - ! + ! ! Just copy everything from xw ! fndmaxup=.true. @@ -1155,13 +1155,13 @@ contains if (widx == jmaxup) fndmaxup=.true. end do end if - if ((i blck else - allocate(blck_,stat=info) - if (info == psb_success_) call blck_%csall(izero,izero,info,ione) + allocate(blck_,stat=info) + if (info == psb_success_) call blck_%allocate(izero,izero,info,ione,type='CSR') if (info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='csall' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if endif - if (present(upd)) then - upd_ = psb_toupper(upd) + if (present(upd)) then + upd_ = psb_toupper(upd) else upd_ = 'F' end if - + m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& - & (m > size(d)) ) then + & (m > size(d)) ) then write(0,*) 'Wrong allocation status for L,D,U? ',& & l%get_nrows(),size(d),u%get_nrows() info = -1 @@ -212,12 +212,12 @@ subroutine psb_dilu0_fact(ialg,a,l,u,d,info,blck, upd) call u%set_triangle() call u%set_unit() call u%set_upper() - + ! ! Nullify pointer / deallocate memory ! - if (present(blck)) then - blck_ => null() + if (present(blck)) then + blck_ => null() else call blck_%free() if(info.ne.0) then @@ -226,7 +226,7 @@ subroutine psb_dilu0_fact(ialg,a,l,u,d,info,blck, upd) call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - deallocate(blck_) + deallocate(blck_) endif call psb_erractionrestore(err_act) @@ -243,9 +243,9 @@ contains ! Version: real ! Note: internal subroutine of psb_dilu0_fact. ! - ! This routine computes either the ILU(0) or the MILU(0) factorization of the - ! diagonal blocks of a distributed matrix. - ! These factorizations are used to build the 'base preconditioner' + ! This routine computes either the ILU(0) or the MILU(0) factorization of the + ! diagonal blocks of a distributed matrix. + ! These factorizations are used to build the 'base preconditioner' ! (block-Jacobi preconditioner/solver, Additive Schwarz ! preconditioner) corresponding to a given level of a multilevel preconditioner. ! @@ -257,7 +257,7 @@ contains ! ! The routine copies and factors "on the fly" from the sparse matrix structures a ! and b into the arrays lval, uval, d (L, U without its diagonal, diagonal of U). - ! + ! ! ! Arguments: ! ialg - integer, input. @@ -296,7 +296,7 @@ contains ! according to the CSR storage format. ! lirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the L factor in lval, according to the CSR storage format. + ! of the L factor in lval, according to the CSR storage format. ! uval - real(psb_dpk_), dimension(:), input/output. ! The U factor in the incomplete factorization. ! The entries of U are stored according to the CSR format. @@ -305,18 +305,18 @@ contains ! according to the CSR storage format. ! uirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the U factor in uval, according to the CSR storage format. + ! of the U factor in uval, according to the CSR storage format. ! l1 - integer, output. ! The number of nonzero entries in lval. ! l2 - integer, output. ! The number of nonzero entries in uval. - ! info - integer, output. + ! info - integer, output. ! Error code. ! subroutine psb_dilu0_factint(ialg,a,b,& & d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,info) - implicit none + implicit none ! Arguments integer(psb_ipk_), intent(in) :: ialg @@ -332,7 +332,7 @@ contains real(psb_dpk_) :: dia,temp integer(psb_ipk_), parameter :: nrb=16 type(psb_d_coo_sparse_mat) :: trw - integer(psb_ipk_) :: int_err(5) + integer(psb_ipk_) :: int_err(5) character(len=20) :: name, ch_err name='psb_dilu0_factint' @@ -346,10 +346,10 @@ contains select case(ialg) case(psb_ilu_n_,psb_milu_n_) - ! Ok + ! Ok case default info=psb_err_input_asize_invalid_i_ - call psb_errpush(info,name,& + call psb_errpush(info,name,& & i_err=(/ione,ialg,izero,izero,izero/)) goto 9999 end select @@ -363,7 +363,7 @@ contains end if m = ma+mb - if (psb_toupper(upd) == 'F' ) then + if (psb_toupper(upd) == 'F' ) then lirp(1) = 1 uirp(1) = 1 l1 = 0 @@ -379,14 +379,14 @@ contains if (i <= ma) then ! ! Copy the i-th local row of the matrix, stored in a, - ! into lval/d(i)/uval + ! into lval/d(i)/uval ! call ilu_copyin(i,ma,a,i,ione,m,l1,lja,lval,& & d(i),l2,uja,uval,ktrw,trw,upd) else ! ! Copy the i-th local row of the matrix, stored in b - ! (as (i-ma)-th row), into lval/d(i)/uval + ! (as (i-ma)-th row), into lval/d(i)/uval ! call ilu_copyin(i-ma,mb,b,i,ione,m,l1,lja,lval,& & d(i),l2,uja,uval,ktrw,trw,upd) @@ -437,7 +437,7 @@ contains ! dia = dia - temp*uval(jj) cycle updateloop - ! + ! else if (j > i) then ! ! search u(i,*) (i-th row of U) for a matching index j @@ -454,23 +454,23 @@ contains end if enddo end if - ! - ! If we get here we missed the cycle updateloop, which means + ! + ! If we get here we missed the cycle updateloop, which means ! that this entry does not match; thus we accumulate on the ! diagonal for MILU(0). ! - if (ialg == psb_milu_n_) then + if (ialg == psb_milu_n_) then dia = dia - temp*uval(jj) end if enddo updateloop enddo - ! + ! ! Check the pivot size - ! + ! if (abs(dia) < d_epstol) then ! ! Too small pivot: unstable factorization - ! + ! info = psb_err_pivot_too_small_ int_err(1) = i write(ch_err,'(g20.10)') abs(dia) @@ -493,7 +493,7 @@ contains else write(0,*) 'Update not implemented ' info = 31 - call psb_errpush(info,name,& + call psb_errpush(info,name,& & i_err=(/ione*13,izero,izero,izero,izero/),a_err=upd) goto 9999 @@ -515,7 +515,7 @@ contains ! Version: real ! Note: internal subroutine of psb_dilu0_fact ! - ! This routine copies a row of a sparse matrix A, stored in the psb_dspmat_type + ! This routine copies a row of a sparse matrix A, stored in the psb_dspmat_type ! data structure a, into the arrays lval and uval and into the scalar variable ! dia, corresponding to the lower and upper triangles of A and to the diagonal ! entry of the row, respectively. The entries in lval and uval are stored @@ -529,14 +529,14 @@ contains ! ! The routine is used by psb_dilu0_factint in the computation of the ILU(0)/MILU(0) ! factorization of a local sparse matrix. - ! + ! ! TODO: modify the routine to allow copying into output L and U that are ! already filled with indices; this would allow computing an ILU(k) pattern, - ! then use the ILU(0) internal for subsequent calls with the same pattern. + ! then use the ILU(0) internal for subsequent calls with the same pattern. ! ! Arguments: ! i - integer, input. - ! The local index of the row to be extracted from the + ! The local index of the row to be extracted from the ! sparse matrix structure a. ! m - integer, input. ! The number of rows of the local matrix stored into a. @@ -570,13 +570,13 @@ contains ! to the CSR storage format. ! uval - real(psb_dpk_), dimension(:), input/output. ! The array where the entries of the row corresponding to the - ! upper triangle are copied. + ! upper triangle are copied. ! ktrw - integer, input/output. ! The index identifying the last entry taken from the ! staging buffer trw. See below. ! trw - type(psb_dspmat_type), input/output. ! A staging buffer. If the matrix A is not in CSR format, we use - ! the psb_sp_getblk routine and store its output in trw; when we + ! the psb_sp_getblk routine and store its output in trw; when we ! need to call psb_sp_getblk we do it for a block of rows, and then ! we consume them from trw in successive calls to this routine, ! until we empty the buffer. Thus we will make a call to psb_sp_getblk @@ -608,10 +608,10 @@ contains if (psb_errstatus_fatal()) then info = psb_err_internal_error_; goto 9999 end if - if (psb_toupper(upd) == 'F') then + if (psb_toupper(upd) == 'F') then - select type(aa => a%a) - type is (psb_d_csr_sparse_mat) + select type(aa => a%a) + type is (psb_d_csr_sparse_mat) ! ! Take a fast shortcut if the matrix is stored in CSR format @@ -633,19 +633,19 @@ contains end if enddo - class default + class default ! - ! Otherwise use psb_sp_getblk, slower but able (in principle) of + ! Otherwise use psb_sp_getblk, slower but able (in principle) of ! handling any format. In this case, a block of rows is extracted ! instead of a single row, for performance reasons, and these ! rows are copied one by one into lval, dia, uval, through ! successive calls to ilu_copyin. ! - if ((mod(i,nrb) == 1).or.(nrb == 1)) then + if ((mod(i,nrb) == 1).or.(nrb == 1)) then irb = min(m-i+1,nrb) - call aa%csget(i,i+irb-1,trw,info) + call aa%csget(i,i+irb-1,trw,info) if(info /= psb_success_) then info=psb_err_from_subroutine_ ch_err='csget' @@ -656,7 +656,7 @@ contains end if nz = trw%get_nzeros() - do + do if (ktrw > nz) exit if (trw%ia(ktrw) > i) exit k = trw%ja(ktrw) @@ -680,7 +680,7 @@ contains write(0,*) 'Update not implemented ' info = 31 - call psb_errpush(info,name,& + call psb_errpush(info,name,& & i_err=(/ione*13,izero,izero,izero,izero/),a_err=upd) goto 9999 diff --git a/prec/impl/psb_d_iluk_fact.f90 b/prec/impl/psb_d_iluk_fact.f90 index 6d644e42..544ec987 100644 --- a/prec/impl/psb_d_iluk_fact.f90 +++ b/prec/impl/psb_d_iluk_fact.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,21 +27,21 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! +! ! Moved here from MLD2P4, original copyright below. -! -! -! +! +! +! ! MLD2P4 version 2.2 ! MultiLevel Domain Decomposition Parallel Preconditioners Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2008-2018 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Daniela di Serafino -! +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -53,7 +53,7 @@ ! 3. The name of the MLD2P4 group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -65,8 +65,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: psb_diluk_fact.f90 ! ! Subroutine: psb_diluk_fact @@ -74,7 +74,7 @@ ! Contains: psb_diluk_factint, iluk_copyin, iluk_fact, iluk_copyout. ! ! This routine computes either the ILU(k) or the MILU(k) factorization of the -! diagonal blocks of a distributed matrix. These factorizations are used to +! diagonal blocks of a distributed matrix. These factorizations are used to ! build the 'base preconditioner' (block-Jacobi preconditioner/solver, ! Additive Schwarz preconditioner) corresponding to a certain level of a ! multilevel preconditioner. @@ -88,7 +88,7 @@ ! U factors is CSR. The diagonal of the U factor is stored separately (actually, ! the inverse of the diagonal entries is stored; this is then managed in the solve ! stage associated to the ILU(k)/MILU(k) factorization). -! +! ! ! Arguments: ! fill_in - integer, input. @@ -118,7 +118,7 @@ ! factorization. ! Note: its allocation is managed by the calling routine psb_ilu_bld, ! hence it cannot be only intent(out). -! info - integer, output. +! info - integer, output. ! Error code. ! blck - type(psb_dspmat_type), input, optional, target. ! The sparse matrix structure containing the remote rows of the @@ -126,7 +126,7 @@ ! to build an Additive Schwarz base preconditioner with overlap ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (see psb_fact_bld), then blck does not contain any row. -! +! subroutine psb_diluk_fact(fill_in,ialg,a,l,u,d,info,blck) use psb_base_mod @@ -143,7 +143,7 @@ subroutine psb_diluk_fact(fill_in,ialg,a,l,u,d,info,blck) real(psb_dpk_), intent(inout) :: d(:) ! Local Variables integer(psb_ipk_) :: l1, l2, m, err_act - + type(psb_dspmat_type), pointer :: blck_ type(psb_d_csr_sparse_mat) :: ll, uu character(len=20) :: name, ch_err @@ -152,17 +152,17 @@ subroutine psb_diluk_fact(fill_in,ialg,a,l,u,d,info,blck) info = psb_success_ call psb_erractionsave(err_act) - ! + ! ! Point to / allocate memory for the incomplete factorization ! - if (present(blck)) then + if (present(blck)) then blck_ => blck else - allocate(blck_,stat=info) - if (info == psb_success_) call blck_%csall(izero,izero,info,ione) + allocate(blck_,stat=info) + if (info == psb_success_) call blck_%allocate(izero,izero,info,ione,type='CSR') if (info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='csall' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if @@ -170,7 +170,7 @@ subroutine psb_diluk_fact(fill_in,ialg,a,l,u,d,info,blck) m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& - & (m > size(d)) ) then + & (m > size(d)) ) then write(0,*) 'Wrong allocation status for L,D,U? ',& & l%get_nrows(),size(d),u%get_nrows() info = -1 @@ -203,15 +203,15 @@ subroutine psb_diluk_fact(fill_in,ialg,a,l,u,d,info,blck) call u%set_triangle() call u%set_unit() call u%set_upper() - + ! ! Nullify pointer / deallocate memory ! - if (present(blck)) then - blck_ => null() + if (present(blck)) then + blck_ => null() else call blck_%free() - deallocate(blck_,stat=info) + deallocate(blck_,stat=info) if(info.ne.0) then info=psb_err_from_subroutine_ ch_err='psb_sp_free' @@ -280,7 +280,7 @@ contains ! according to the CSR storage format. ! lia2 - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the L factor in laspk, according to the CSR storage format. + ! of the L factor in laspk, according to the CSR storage format. ! uval - real(psb_dpk_), dimension(:), input/output. ! The U factor in the incomplete factorization. ! The entries of U are stored according to the CSR format. @@ -289,12 +289,12 @@ contains ! according to the CSR storage format. ! uirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the U factor in uval, according to the CSR storage format. + ! of the U factor in uval, according to the CSR storage format. ! l1 - integer, output ! The number of nonzero entries in laspk. ! l2 - integer, output ! The number of nonzero entries in uval. - ! info - integer, output. + ! info - integer, output. ! Error code. ! subroutine psb_diluk_factint(fill_in,ialg,a,b,& @@ -304,7 +304,7 @@ contains implicit none - ! Arguments + ! Arguments integer(psb_ipk_), intent(in) :: fill_in, ialg type(psb_dspmat_type),intent(in) :: a,b integer(psb_ipk_),intent(inout) :: l1,l2,info @@ -330,14 +330,14 @@ contains select case(ialg) case(psb_ilu_n_,psb_milu_n_) - ! Ok + ! Ok case default info=psb_err_input_asize_invalid_i_ call psb_errpush(info,name,& & i_err=(/itwo,ialg,izero,izero,izero/)) goto 9999 end select - if (fill_in < 0) then + if (fill_in < 0) then info=psb_err_input_asize_invalid_i_ call psb_errpush(info,name, & & i_err=(/ione,fill_in,izero,izero,izero/)) @@ -349,7 +349,7 @@ contains m = ma+mb ! - ! Allocate a temporary buffer for the iluk_copyin function + ! Allocate a temporary buffer for the iluk_copyin function ! call trw%allocate(izero,izero,ione) @@ -361,7 +361,7 @@ contains call psb_errpush(info,name,a_err='psb_sp_all') goto 9999 end if - + l1=0 l2=0 lirp(1) = 1 @@ -381,34 +381,34 @@ contains uplevs(:) = m+1 row(:) = dzero rowlevs(:) = -(m+1) - + ! ! Cycle over the matrix rows ! do i = 1, m - + ! ! At each iteration of the loop we keep in a heap the column indices ! affected by the factorization. The heap is initialized and filled ! in the iluk_copyin routine, and updated during the elimination, in ! the iluk_fact routine. The heap is ideal because at each step we need ! the lowest index, but we also need to insert new items, and the heap - ! allows to do both in log time. + ! allows to do both in log time. ! d(i) = dzero if (i<=ma) then ! - ! Copy into trw the i-th local row of the matrix, stored in a - ! + ! Copy into trw the i-th local row of the matrix, stored in a + ! call iluk_copyin(i,ma,a,ione,m,row,rowlevs,heap,ktrw,trw,info) else ! ! Copy into trw the i-th local row of the matrix, stored in b - ! (as (i-ma)-th row) - ! + ! (as (i-ma)-th row) + ! call iluk_copyin(i-ma,mb,b,ione,m,row,rowlevs,heap,ktrw,trw,info) endif - + ! Do an elimination step on the current row. It turns out we only ! need to keep track of fill levels for the upper triangle, hence we ! do not have a lowlevs variable. @@ -417,7 +417,7 @@ contains & d,uja,uirp,uval,uplevs,nidx,idxs,info) ! ! Copy the row into lval/d(i)/uval - ! + ! if (info == psb_success_) call iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,& & l1,l2,lja,lirp,lval,d,uja,uirp,uval,uplevs,info) if (info /= psb_success_) then @@ -474,11 +474,11 @@ contains ! ! This routine is used by psb_diluk_factint in the computation of the ! ILU(k)/MILU(k) factorization of a local sparse matrix. - ! + ! ! ! Arguments: ! i - integer, input. - ! The local index of the row to be extracted from the + ! The local index of the row to be extracted from the ! sparse matrix structure a. ! m - integer, input. ! The number of rows of the local matrix stored into a. @@ -510,7 +510,7 @@ contains ! staging buffer trw. See below. ! trw - type(psb_dspmat_type), input/output. ! A staging buffer. If the matrix A is not in CSR format, we use - ! the psb_sp_getblk routine and store its output in trw; when we + ! the psb_sp_getblk routine and store its output in trw; when we ! need to call psb_sp_getblk we do it for a block of rows, and then ! we consume them from trw in successive calls to this routine, ! until we empty the buffer. Thus we will make a call to psb_sp_getblk @@ -521,8 +521,8 @@ contains use psb_base_mod implicit none - - ! Arguments + + ! Arguments type(psb_dspmat_type), intent(in) :: a type(psb_d_coo_sparse_mat), intent(inout) :: trw integer(psb_ipk_), intent(in) :: i,m,jmin,jmax @@ -542,17 +542,17 @@ contains if (psb_errstatus_fatal()) then info = psb_err_internal_error_; goto 9999 end if - call heap%init(info) + call heap%init(info) - select type (aa=> a%a) - type is (psb_d_csr_sparse_mat) + select type (aa=> a%a) + type is (psb_d_csr_sparse_mat) ! ! Take a fast shortcut if the matrix is stored in CSR format ! - + do j = aa%irp(i), aa%irp(i+1) - 1 k = aa%ja(j) - if ((jmin<=k).and.(k<=jmax)) then + if ((jmin<=k).and.(k<=jmax)) then row(k) = aa%val(j) rowlevs(k) = 0 call heap%insert(k,info) @@ -562,14 +562,14 @@ contains class default ! - ! Otherwise use psb_sp_getblk, slower but able (in principle) of + ! Otherwise use psb_sp_getblk, slower but able (in principle) of ! handling any format. In this case, a block of rows is extracted ! instead of a single row, for performance reasons, and these ! rows are copied one by one into the array row, through successive ! calls to iluk_copyin. ! - if ((mod(i,nrb) == 1).or.(nrb == 1)) then + if ((mod(i,nrb) == 1).or.(nrb == 1)) then irb = min(m-i+1,nrb) call aa%csget(i,i+irb-1,trw,info) if (info /= psb_success_) then @@ -581,11 +581,11 @@ contains ktrw=1 end if nz = trw%get_nzeros() - do + do if (ktrw > nz) exit if (trw%ia(ktrw) > i) exit k = trw%ja(ktrw) - if ((jmin<=k).and.(k<=jmax)) then + if ((jmin<=k).and.(k<=jmax)) then row(k) = trw%val(ktrw) rowlevs(k) = 0 call heap%insert(k,info) @@ -674,10 +674,10 @@ contains use psb_base_mod - implicit none + implicit none ! Arguments - type(psb_i_heap), intent(inout) :: heap + type(psb_i_heap), intent(inout) :: heap integer(psb_ipk_), intent(in) :: i, fill_in integer(psb_ipk_), intent(inout) :: nidx,info integer(psb_ipk_), intent(inout) :: rowlevs(:) @@ -690,7 +690,7 @@ contains real(psb_dpk_) :: rwk info = psb_success_ - if (.not.allocated(idxs)) then + if (.not.allocated(idxs)) then allocate(idxs(200),stat=info) if (info /= psb_success_) return endif @@ -702,34 +702,34 @@ contains ! do ! Beware: (iret < 0) means that the heap is empty, not an error. - call heap%get_first(k,iret) + call heap%get_first(k,iret) if (iret < 0) return - ! + ! ! Just in case an index has been put on the heap more than once. ! if (k == lastk) cycle - lastk = k + lastk = k nidx = nidx + 1 - if (nidx>size(idxs)) then + if (nidx>size(idxs)) then call psb_realloc(nidx+psb_heap_resize,idxs,info) if (info /= psb_success_) return end if idxs(nidx) = k - - if ((row(k) /= dzero).and.(rowlevs(k) <= fill_in).and.(ki) then + else if (j>i) then ! ! Copy the upper part of the row - ! - if (rowlevs(j) <= fill_in) then - l2 = l2 + 1 - if (size(uval) < l2) then - ! + ! + if (rowlevs(j) <= fill_in) then + l2 = l2 + 1 + if (size(uval) < l2) then + ! ! Figure out a good reallocation size! ! isz = max((l2/i)*m,int(1.2*l2),l2+100) - call psb_realloc(isz,uval,info) - if (info == psb_success_) call psb_realloc(isz,uja,info) + call psb_realloc(isz,uval,info) + if (info == psb_success_) call psb_realloc(isz,uja,info) if (info == psb_success_) call psb_realloc(isz,uplevs,info,pad=(m+1)) - if (info /= psb_success_) then + if (info /= psb_success_) then info=psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') goto 9999 @@ -941,11 +941,11 @@ contains end if uja(l2) = j uval(l2) = row(j) - uplevs(l2) = rowlevs(j) + uplevs(l2) = rowlevs(j) else if (ialg == psb_milu_n_) then ! ! MILU(k): add discarded entries to the diagonal one - ! + ! d(i) = d(i) + row(j) end if ! @@ -963,13 +963,13 @@ contains lirp(i+1) = l1 + 1 uirp(i+1) = l2 + 1 - ! + ! ! Check the pivot size ! if (abs(d(i)) < d_epstol) then ! ! Too small pivot: unstable factorization - ! + ! info = psb_err_pivot_too_small_ int_err(1) = i write(ch_err,'(g20.10)') d(i) diff --git a/prec/impl/psb_d_ilut_fact.f90 b/prec/impl/psb_d_ilut_fact.f90 index bcd26396..6c2dc698 100644 --- a/prec/impl/psb_d_ilut_fact.f90 +++ b/prec/impl/psb_d_ilut_fact.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,21 +27,21 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! +! ! Moved here from MLD2P4, original copyright below. -! -! -! +! +! +! ! MLD2P4 version 2.2 ! MultiLevel Domain Decomposition Parallel Preconditioners Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2008-2018 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Daniela di Serafino -! +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -53,7 +53,7 @@ ! 3. The name of the MLD2P4 group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -65,8 +65,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: psb_dilut_fact.f90 ! ! Subroutine: psb_dilut_fact @@ -87,7 +87,7 @@ ! CSR. The diagonal of the U factor is stored separately (actually, the ! inverse of the diagonal entries is stored; this is then managed in the ! solve stage associated to the ILU(k,t) factorization). -! +! ! ! Arguments: ! fill_in - integer, input. @@ -114,7 +114,7 @@ ! factorization. ! Note: its allocation is managed by the calling routine psb_ilu_bld, ! hence it cannot be only intent(out). -! info - integer, output. +! info - integer, output. ! Error code. ! blck - type(psb_dspmat_type), input, optional, target. ! The sparse matrix structure containing the remote rows of the @@ -122,9 +122,9 @@ ! to build an Additive Schwarz base preconditioner with overlap ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (see psb_fact_bld), then blck does not contain any row. -! +! subroutine psb_dilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) - + use psb_base_mod use psb_d_ilu_fact_mod, psb_protect_name => psb_dilut_fact @@ -141,7 +141,7 @@ subroutine psb_dilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) integer(psb_ipk_), intent(in), optional :: iscale ! Local Variables integer(psb_ipk_) :: l1, l2, m, err_act, iscale_ - + type(psb_dspmat_type), pointer :: blck_ type(psb_d_csr_sparse_mat) :: ll, uu real(psb_dpk_) :: scale @@ -151,34 +151,34 @@ subroutine psb_dilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) info = psb_success_ call psb_erractionsave(err_act) - if (fill_in < 0) then + if (fill_in < 0) then info=psb_err_input_asize_invalid_i_ call psb_errpush(info,name, & & i_err=(/ione,fill_in,izero,izero,izero/)) goto 9999 end if - ! + ! ! Point to / allocate memory for the incomplete factorization ! - if (present(blck)) then + if (present(blck)) then blck_ => blck else - allocate(blck_,stat=info) - if (info == psb_success_) call blck_%csall(izero,izero,info,ione) + allocate(blck_,stat=info) + if (info == psb_success_) call blck_%allocate(izero,izero,info,ione,type='CSR') if (info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='csall' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if endif - if (present(iscale)) then + if (present(iscale)) then iscale_ = iscale else iscale_ = psb_ilu_scale_none_ end if - - select case(iscale_) + + select case(iscale_) case(psb_ilu_scale_none_) scale = sone case(psb_ilu_scale_maxval_) @@ -189,10 +189,10 @@ subroutine psb_dilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) call psb_errpush(info,name,i_err=(/ione*9,iscale_,izero,izero,izero/)) goto 9999 end select - + m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& - & (m > size(d)) ) then + & (m > size(d)) ) then write(0,*) 'Wrong allocation status for L,D,U? ',& & l%get_nrows(),size(d),u%get_nrows() info = -1 @@ -225,15 +225,15 @@ subroutine psb_dilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) call u%set_triangle() call u%set_unit() call u%set_upper() - + ! ! Nullify pointer / deallocate memory ! - if (present(blck)) then - blck_ => null() + if (present(blck)) then + blck_ => null() else call blck_%free() - deallocate(blck_,stat=info) + deallocate(blck_,stat=info) if(info.ne.0) then info=psb_err_from_subroutine_ ch_err='psb_sp_free' @@ -298,7 +298,7 @@ contains ! according to the CSR storage format. ! lirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the L factor in lval, according to the CSR storage format. + ! of the L factor in lval, according to the CSR storage format. ! uval - real(psb_dpk_), dimension(:), input/output. ! The U factor in the incomplete factorization. ! The entries of U are stored according to the CSR format. @@ -307,12 +307,12 @@ contains ! according to the CSR storage format. ! uirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the U factor in uval, according to the CSR storage format. + ! of the U factor in uval, according to the CSR storage format. ! l1 - integer, output ! The number of nonzero entries in lval. ! l2 - integer, output ! The number of nonzero entries in uval. - ! info - integer, output. + ! info - integer, output. ! Error code. ! subroutine psb_dilut_factint(fill_in,thres,a,b,& @@ -320,7 +320,7 @@ contains use psb_base_mod - implicit none + implicit none ! Arguments integer(psb_ipk_), intent(in) :: fill_in @@ -355,7 +355,7 @@ contains m = ma+mb ! - ! Allocate a temporary buffer for the ilut_copyin function + ! Allocate a temporary buffer for the ilut_copyin function ! call trw%allocate(izero,izero,ione) if (info == psb_success_) call psb_ensure_size(m+1,lirp,info) @@ -366,7 +366,7 @@ contains call psb_errpush(info,name,a_err='psb_sp_all') goto 9999 end if - + l1=0 l2=0 lirp(1) = 1 @@ -396,10 +396,10 @@ contains ! in the ilut_copyin function, and updated during the elimination, in ! the ilut_fact routine. The heap is ideal because at each step we need ! the lowest index, but we also need to insert new items, and the heap - ! allows to do both in log time. + ! allows to do both in log time. ! d(i) = czero - if (i<=ma) then + if (i<=ma) then call ilut_copyin(i,ma,a,i,ione,m,nlw,nup,jmaxup,nrmi,weight,& & row,heap,ktrw,trw,info) else @@ -414,7 +414,7 @@ contains & d,uja,uirp,uval,nidx,idxs,info) ! ! Copy the row into lval/d(i)/uval - ! + ! if (info == psb_success_) call ilut_copyout(fill_in,thres,i,m,& & nlw,nup,jmaxup,nrmi,row,nidx,idxs,& & l1,l2,lja,lirp,lval,d,uja,uirp,uval,info) @@ -429,7 +429,7 @@ contains ! ! Adjust diagonal accounting for scale factor ! - if (weight /= sone) then + if (weight /= sone) then d(1:m) = d(1:m)*weight end if @@ -469,7 +469,7 @@ contains ! - storing into a heap the column indices of the nonzero entries of the copied ! row; ! - computing the column index of the first entry with maximum absolute value - ! in the part of the row belonging to the upper triangle; + ! in the part of the row belonging to the upper triangle; ! - computing the 2-norm of the row. ! The output array row is such that it contains a full row of A, i.e. it contains ! also the zero entries of the row. This is useful for the elimination step @@ -482,11 +482,11 @@ contains ! ! This routine is used by psb_dilut_factint in the computation of the ILU(k,t) ! factorization of a local sparse matrix. - ! + ! ! ! Arguments: ! i - integer, input. - ! The local index of the row to be extracted from the + ! The local index of the row to be extracted from the ! sparse matrix structure a. ! m - integer, input. ! The number of rows of the local matrix stored into a. @@ -533,7 +533,7 @@ contains ! staging buffer trw. See below. ! trw - type(psb_dspmat_type), input/output. ! A staging buffer. If the matrix A is not in CSR format, we use - ! the psb_sp_getblk routine and store its output in trw; when we + ! the psb_sp_getblk routine and store its output in trw; when we ! need to call psb_sp_getblk we do it for a block of rows, and then ! we consume them from trw in successive calls to this routine, ! until we empty the buffer. Thus we will make a call to psb_sp_getblk @@ -542,7 +542,7 @@ contains subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,& & nrmi,weight,row,heap,ktrw,trw,info) use psb_base_mod - implicit none + implicit none type(psb_dspmat_type), intent(in) :: a type(psb_d_coo_sparse_mat), intent(inout) :: trw integer(psb_ipk_), intent(in) :: i, m,jmin,jmax,jd @@ -584,29 +584,29 @@ contains dmaxup = szero nrmi = szero - select type (aa=> a%a) - type is (psb_d_csr_sparse_mat) + select type (aa=> a%a) + type is (psb_d_csr_sparse_mat) ! ! Take a fast shortcut if the matrix is stored in CSR format ! do j = aa%irp(i), aa%irp(i+1) - 1 k = aa%ja(j) - if ((jmin<=k).and.(k<=jmax)) then - row(k) = aa%val(j)*weight + if ((jmin<=k).and.(k<=jmax)) then + row(k) = aa%val(j)*weight call heap%insert(k,info) if (info /= psb_success_) exit - if (kjd) then + if (kjd) then nup = nup + 1 - if (abs(row(k))>dmaxup) then + if (abs(row(k))>dmaxup) then jmaxup = k dmaxup = abs(row(k)) end if end if end if end do - if (info /= psb_success_) then + if (info /= psb_success_) then info=psb_err_from_subroutine_ call psb_errpush(info,name,a_err='psb_insert_heap') goto 9999 @@ -617,16 +617,16 @@ contains class default - + ! - ! Otherwise use psb_sp_getblk, slower but able (in principle) of + ! Otherwise use psb_sp_getblk, slower but able (in principle) of ! handling any format. In this case, a block of rows is extracted ! instead of a single row, for performance reasons, and these ! rows are copied one by one into the array row, through successive ! calls to ilut_copyin. ! - if ((mod(i,nrb) == 1).or.(nrb == 1)) then + if ((mod(i,nrb) == 1).or.(nrb == 1)) then irb = min(m-i+1,nrb) call aa%csget(i,i+irb-1,trw,info) if (info /= psb_success_) then @@ -639,18 +639,18 @@ contains kin = ktrw nz = trw%get_nzeros() - do + do if (ktrw > nz) exit if (trw%ia(ktrw) > i) exit k = trw%ja(ktrw) - if ((jmin<=k).and.(k<=jmax)) then + if ((jmin<=k).and.(k<=jmax)) then row(k) = trw%val(ktrw)*weight call heap%insert(k,info) if (info /= psb_success_) exit - if (kjd) then + if (kjd) then nup = nup + 1 - if (abs(row(k))>dmaxup) then + if (abs(row(k))>dmaxup) then jmaxup = k dmaxup = abs(row(k)) end if @@ -734,10 +734,10 @@ contains use psb_base_mod - implicit none + implicit none ! Arguments - type(psb_i_heap), intent(inout) :: heap + type(psb_i_heap), intent(inout) :: heap integer(psb_ipk_), intent(in) :: i integer(psb_ipk_), intent(inout) :: nidx,info real(psb_dpk_), intent(in) :: thres,nrmi @@ -753,23 +753,23 @@ contains call psb_ensure_size(200*ione,idxs,info) if (info /= psb_success_) return nidx = 0 - lastk = -1 + lastk = -1 ! ! Do while there are indices to be processed ! do - call heap%get_first(k,iret) + call heap%get_first(k,iret) if (iret < 0) exit - ! + ! ! An index may have been put on the heap more than once. ! if (k == lastk) cycle - lastk = k + lastk = k lowert: if (k nidx) exit if (idxs(idxp) >= i) exit @@ -987,8 +987,8 @@ contains ! if (abs(witem) < thres*nrmi) cycle - nz = nz + 1 - xw(nz) = witem + nz = nz + 1 + xw(nz) = witem xwid(nz) = widx call heap%insert(witem,widx,info) if (info /= psb_success_) then @@ -1000,9 +1000,9 @@ contains ! ! Now we have to take out the first nlw+fill_in entries - ! + ! if (nz <= nlw+fill_in) then - ! + ! ! Just copy everything from xw, and it is already ordered ! else @@ -1014,7 +1014,7 @@ contains call psb_errpush(info,name,a_err='psb_heap_get_first') goto 9999 end if - + xw(k) = witem xwid(k) = widx end do @@ -1029,15 +1029,15 @@ contains ! Copy out the lower part of the row ! do k=1,nz - l1 = l1 + 1 + l1 = l1 + 1 if (size(lval) < l1) then - ! + ! ! Figure out a good reallocation size! - ! + ! isz = (max((l1/i)*m,int(1.2*l1),l1+100)) - call psb_realloc(isz,lval,info) - if (info == psb_success_) call psb_realloc(isz,lja,info) - if (info /= psb_success_) then + call psb_realloc(isz,lval,info) + if (info == psb_success_) call psb_realloc(isz,lja,info) + if (info /= psb_success_) then info=psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') goto 9999 @@ -1046,25 +1046,25 @@ contains lja(l1) = xwid(k) lval(l1) = xw(indx(k)) end do - + ! ! Make sure idxp points to the diagonal entry ! - if (idxp <= size(idxs)) then - if (idxs(idxp) < i) then - do + if (idxp <= size(idxs)) then + if (idxs(idxp) < i) then + do idxp = idxp + 1 if (idxp > nidx) exit if (idxs(idxp) >= i) exit end do end if end if - if (idxp > size(idxs)) then + if (idxp > size(idxs)) then !!$ write(0,*) 'Warning: missing diagonal element in the row ' else - if (idxs(idxp) > i) then + if (idxs(idxp) > i) then !!$ write(0,*) 'Warning: missing diagonal element in the row ' - else if (idxs(idxp) /= i) then + else if (idxs(idxp) /= i) then !!$ write(0,*) 'Warning: impossible error: diagonal has vanished' else ! @@ -1076,7 +1076,7 @@ contains if (abs(d(i)) < d_epstol) then ! ! Too small pivot: unstable factorization - ! + ! info = psb_err_pivot_too_small_ int_err(1) = i write(ch_err,'(g20.10)') d(i) @@ -1092,7 +1092,7 @@ contains end if ! - ! Now the upper part + ! Now the upper part ! call heap%init(info,dir=psb_asort_down_) @@ -1104,15 +1104,15 @@ contains nz = 0 do - + idxp = idxp + 1 if (idxp > nidx) exit widx = idxs(idxp) - if (widx <= i) then + if (widx <= i) then !!$ write(0,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) cycle end if - if (widx > m) then + if (widx > m) then !!$ write(0,*) 'Warning: impossible value',widx,i,idxp,idxs(idxp) cycle end if @@ -1120,12 +1120,12 @@ contains ! ! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway. ! - if ((widx /= jmaxup) .and. (abs(witem) < thres*nrmi)) then - cycle + if ((widx /= jmaxup) .and. (abs(witem) < thres*nrmi)) then + cycle end if nz = nz + 1 - xw(nz) = witem + xw(nz) = witem xwid(nz) = widx call heap%insert(witem,widx,info) if (info /= psb_success_) then @@ -1133,7 +1133,7 @@ contains call psb_errpush(info,name,a_err='psb_insert_heap') goto 9999 end if - + end do ! @@ -1141,7 +1141,7 @@ contains ! we include entry jmaxup. ! if (nz <= nup+fill_in) then - ! + ! ! Just copy everything from xw ! fndmaxup=.true. @@ -1155,13 +1155,13 @@ contains if (widx == jmaxup) fndmaxup=.true. end do end if - if ((i blck else - allocate(blck_,stat=info) - if (info == psb_success_) call blck_%csall(izero,izero,info,ione) + allocate(blck_,stat=info) + if (info == psb_success_) call blck_%allocate(izero,izero,info,ione,type='CSR') if (info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='csall' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if endif - if (present(upd)) then - upd_ = psb_toupper(upd) + if (present(upd)) then + upd_ = psb_toupper(upd) else upd_ = 'F' end if - + m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& - & (m > size(d)) ) then + & (m > size(d)) ) then write(0,*) 'Wrong allocation status for L,D,U? ',& & l%get_nrows(),size(d),u%get_nrows() info = -1 @@ -212,12 +212,12 @@ subroutine psb_silu0_fact(ialg,a,l,u,d,info,blck, upd) call u%set_triangle() call u%set_unit() call u%set_upper() - + ! ! Nullify pointer / deallocate memory ! - if (present(blck)) then - blck_ => null() + if (present(blck)) then + blck_ => null() else call blck_%free() if(info.ne.0) then @@ -226,7 +226,7 @@ subroutine psb_silu0_fact(ialg,a,l,u,d,info,blck, upd) call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - deallocate(blck_) + deallocate(blck_) endif call psb_erractionrestore(err_act) @@ -243,9 +243,9 @@ contains ! Version: real ! Note: internal subroutine of psb_silu0_fact. ! - ! This routine computes either the ILU(0) or the MILU(0) factorization of the - ! diagonal blocks of a distributed matrix. - ! These factorizations are used to build the 'base preconditioner' + ! This routine computes either the ILU(0) or the MILU(0) factorization of the + ! diagonal blocks of a distributed matrix. + ! These factorizations are used to build the 'base preconditioner' ! (block-Jacobi preconditioner/solver, Additive Schwarz ! preconditioner) corresponding to a given level of a multilevel preconditioner. ! @@ -257,7 +257,7 @@ contains ! ! The routine copies and factors "on the fly" from the sparse matrix structures a ! and b into the arrays lval, uval, d (L, U without its diagonal, diagonal of U). - ! + ! ! ! Arguments: ! ialg - integer, input. @@ -296,7 +296,7 @@ contains ! according to the CSR storage format. ! lirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the L factor in lval, according to the CSR storage format. + ! of the L factor in lval, according to the CSR storage format. ! uval - real(psb_spk_), dimension(:), input/output. ! The U factor in the incomplete factorization. ! The entries of U are stored according to the CSR format. @@ -305,18 +305,18 @@ contains ! according to the CSR storage format. ! uirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the U factor in uval, according to the CSR storage format. + ! of the U factor in uval, according to the CSR storage format. ! l1 - integer, output. ! The number of nonzero entries in lval. ! l2 - integer, output. ! The number of nonzero entries in uval. - ! info - integer, output. + ! info - integer, output. ! Error code. ! subroutine psb_silu0_factint(ialg,a,b,& & d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,info) - implicit none + implicit none ! Arguments integer(psb_ipk_), intent(in) :: ialg @@ -332,7 +332,7 @@ contains real(psb_spk_) :: dia,temp integer(psb_ipk_), parameter :: nrb=16 type(psb_s_coo_sparse_mat) :: trw - integer(psb_ipk_) :: int_err(5) + integer(psb_ipk_) :: int_err(5) character(len=20) :: name, ch_err name='psb_silu0_factint' @@ -346,10 +346,10 @@ contains select case(ialg) case(psb_ilu_n_,psb_milu_n_) - ! Ok + ! Ok case default info=psb_err_input_asize_invalid_i_ - call psb_errpush(info,name,& + call psb_errpush(info,name,& & i_err=(/ione,ialg,izero,izero,izero/)) goto 9999 end select @@ -363,7 +363,7 @@ contains end if m = ma+mb - if (psb_toupper(upd) == 'F' ) then + if (psb_toupper(upd) == 'F' ) then lirp(1) = 1 uirp(1) = 1 l1 = 0 @@ -379,14 +379,14 @@ contains if (i <= ma) then ! ! Copy the i-th local row of the matrix, stored in a, - ! into lval/d(i)/uval + ! into lval/d(i)/uval ! call ilu_copyin(i,ma,a,i,ione,m,l1,lja,lval,& & d(i),l2,uja,uval,ktrw,trw,upd) else ! ! Copy the i-th local row of the matrix, stored in b - ! (as (i-ma)-th row), into lval/d(i)/uval + ! (as (i-ma)-th row), into lval/d(i)/uval ! call ilu_copyin(i-ma,mb,b,i,ione,m,l1,lja,lval,& & d(i),l2,uja,uval,ktrw,trw,upd) @@ -437,7 +437,7 @@ contains ! dia = dia - temp*uval(jj) cycle updateloop - ! + ! else if (j > i) then ! ! search u(i,*) (i-th row of U) for a matching index j @@ -454,23 +454,23 @@ contains end if enddo end if - ! - ! If we get here we missed the cycle updateloop, which means + ! + ! If we get here we missed the cycle updateloop, which means ! that this entry does not match; thus we accumulate on the ! diagonal for MILU(0). ! - if (ialg == psb_milu_n_) then + if (ialg == psb_milu_n_) then dia = dia - temp*uval(jj) end if enddo updateloop enddo - ! + ! ! Check the pivot size - ! + ! if (abs(dia) < s_epstol) then ! ! Too small pivot: unstable factorization - ! + ! info = psb_err_pivot_too_small_ int_err(1) = i write(ch_err,'(g20.10)') abs(dia) @@ -493,7 +493,7 @@ contains else write(0,*) 'Update not implemented ' info = 31 - call psb_errpush(info,name,& + call psb_errpush(info,name,& & i_err=(/ione*13,izero,izero,izero,izero/),a_err=upd) goto 9999 @@ -515,7 +515,7 @@ contains ! Version: real ! Note: internal subroutine of psb_silu0_fact ! - ! This routine copies a row of a sparse matrix A, stored in the psb_sspmat_type + ! This routine copies a row of a sparse matrix A, stored in the psb_sspmat_type ! data structure a, into the arrays lval and uval and into the scalar variable ! dia, corresponding to the lower and upper triangles of A and to the diagonal ! entry of the row, respectively. The entries in lval and uval are stored @@ -529,14 +529,14 @@ contains ! ! The routine is used by psb_silu0_factint in the computation of the ILU(0)/MILU(0) ! factorization of a local sparse matrix. - ! + ! ! TODO: modify the routine to allow copying into output L and U that are ! already filled with indices; this would allow computing an ILU(k) pattern, - ! then use the ILU(0) internal for subsequent calls with the same pattern. + ! then use the ILU(0) internal for subsequent calls with the same pattern. ! ! Arguments: ! i - integer, input. - ! The local index of the row to be extracted from the + ! The local index of the row to be extracted from the ! sparse matrix structure a. ! m - integer, input. ! The number of rows of the local matrix stored into a. @@ -570,13 +570,13 @@ contains ! to the CSR storage format. ! uval - real(psb_spk_), dimension(:), input/output. ! The array where the entries of the row corresponding to the - ! upper triangle are copied. + ! upper triangle are copied. ! ktrw - integer, input/output. ! The index identifying the last entry taken from the ! staging buffer trw. See below. ! trw - type(psb_sspmat_type), input/output. ! A staging buffer. If the matrix A is not in CSR format, we use - ! the psb_sp_getblk routine and store its output in trw; when we + ! the psb_sp_getblk routine and store its output in trw; when we ! need to call psb_sp_getblk we do it for a block of rows, and then ! we consume them from trw in successive calls to this routine, ! until we empty the buffer. Thus we will make a call to psb_sp_getblk @@ -608,10 +608,10 @@ contains if (psb_errstatus_fatal()) then info = psb_err_internal_error_; goto 9999 end if - if (psb_toupper(upd) == 'F') then + if (psb_toupper(upd) == 'F') then - select type(aa => a%a) - type is (psb_s_csr_sparse_mat) + select type(aa => a%a) + type is (psb_s_csr_sparse_mat) ! ! Take a fast shortcut if the matrix is stored in CSR format @@ -633,19 +633,19 @@ contains end if enddo - class default + class default ! - ! Otherwise use psb_sp_getblk, slower but able (in principle) of + ! Otherwise use psb_sp_getblk, slower but able (in principle) of ! handling any format. In this case, a block of rows is extracted ! instead of a single row, for performance reasons, and these ! rows are copied one by one into lval, dia, uval, through ! successive calls to ilu_copyin. ! - if ((mod(i,nrb) == 1).or.(nrb == 1)) then + if ((mod(i,nrb) == 1).or.(nrb == 1)) then irb = min(m-i+1,nrb) - call aa%csget(i,i+irb-1,trw,info) + call aa%csget(i,i+irb-1,trw,info) if(info /= psb_success_) then info=psb_err_from_subroutine_ ch_err='csget' @@ -656,7 +656,7 @@ contains end if nz = trw%get_nzeros() - do + do if (ktrw > nz) exit if (trw%ia(ktrw) > i) exit k = trw%ja(ktrw) @@ -680,7 +680,7 @@ contains write(0,*) 'Update not implemented ' info = 31 - call psb_errpush(info,name,& + call psb_errpush(info,name,& & i_err=(/ione*13,izero,izero,izero,izero/),a_err=upd) goto 9999 diff --git a/prec/impl/psb_s_iluk_fact.f90 b/prec/impl/psb_s_iluk_fact.f90 index 4b9f1f3f..6129663b 100644 --- a/prec/impl/psb_s_iluk_fact.f90 +++ b/prec/impl/psb_s_iluk_fact.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,21 +27,21 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! +! ! Moved here from MLD2P4, original copyright below. -! -! -! +! +! +! ! MLD2P4 version 2.2 ! MultiLevel Domain Decomposition Parallel Preconditioners Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2008-2018 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Daniela di Serafino -! +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -53,7 +53,7 @@ ! 3. The name of the MLD2P4 group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -65,8 +65,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: psb_siluk_fact.f90 ! ! Subroutine: psb_siluk_fact @@ -74,7 +74,7 @@ ! Contains: psb_siluk_factint, iluk_copyin, iluk_fact, iluk_copyout. ! ! This routine computes either the ILU(k) or the MILU(k) factorization of the -! diagonal blocks of a distributed matrix. These factorizations are used to +! diagonal blocks of a distributed matrix. These factorizations are used to ! build the 'base preconditioner' (block-Jacobi preconditioner/solver, ! Additive Schwarz preconditioner) corresponding to a certain level of a ! multilevel preconditioner. @@ -88,7 +88,7 @@ ! U factors is CSR. The diagonal of the U factor is stored separately (actually, ! the inverse of the diagonal entries is stored; this is then managed in the solve ! stage associated to the ILU(k)/MILU(k) factorization). -! +! ! ! Arguments: ! fill_in - integer, input. @@ -118,7 +118,7 @@ ! factorization. ! Note: its allocation is managed by the calling routine psb_ilu_bld, ! hence it cannot be only intent(out). -! info - integer, output. +! info - integer, output. ! Error code. ! blck - type(psb_sspmat_type), input, optional, target. ! The sparse matrix structure containing the remote rows of the @@ -126,7 +126,7 @@ ! to build an Additive Schwarz base preconditioner with overlap ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (see psb_fact_bld), then blck does not contain any row. -! +! subroutine psb_siluk_fact(fill_in,ialg,a,l,u,d,info,blck) use psb_base_mod @@ -143,7 +143,7 @@ subroutine psb_siluk_fact(fill_in,ialg,a,l,u,d,info,blck) real(psb_spk_), intent(inout) :: d(:) ! Local Variables integer(psb_ipk_) :: l1, l2, m, err_act - + type(psb_sspmat_type), pointer :: blck_ type(psb_s_csr_sparse_mat) :: ll, uu character(len=20) :: name, ch_err @@ -152,17 +152,17 @@ subroutine psb_siluk_fact(fill_in,ialg,a,l,u,d,info,blck) info = psb_success_ call psb_erractionsave(err_act) - ! + ! ! Point to / allocate memory for the incomplete factorization ! - if (present(blck)) then + if (present(blck)) then blck_ => blck else - allocate(blck_,stat=info) - if (info == psb_success_) call blck_%csall(izero,izero,info,ione) + allocate(blck_,stat=info) + if (info == psb_success_) call blck_%allocate(izero,izero,info,ione,type='CSR') if (info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='csall' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if @@ -170,7 +170,7 @@ subroutine psb_siluk_fact(fill_in,ialg,a,l,u,d,info,blck) m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& - & (m > size(d)) ) then + & (m > size(d)) ) then write(0,*) 'Wrong allocation status for L,D,U? ',& & l%get_nrows(),size(d),u%get_nrows() info = -1 @@ -203,15 +203,15 @@ subroutine psb_siluk_fact(fill_in,ialg,a,l,u,d,info,blck) call u%set_triangle() call u%set_unit() call u%set_upper() - + ! ! Nullify pointer / deallocate memory ! - if (present(blck)) then - blck_ => null() + if (present(blck)) then + blck_ => null() else call blck_%free() - deallocate(blck_,stat=info) + deallocate(blck_,stat=info) if(info.ne.0) then info=psb_err_from_subroutine_ ch_err='psb_sp_free' @@ -280,7 +280,7 @@ contains ! according to the CSR storage format. ! lia2 - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the L factor in laspk, according to the CSR storage format. + ! of the L factor in laspk, according to the CSR storage format. ! uval - real(psb_spk_), dimension(:), input/output. ! The U factor in the incomplete factorization. ! The entries of U are stored according to the CSR format. @@ -289,12 +289,12 @@ contains ! according to the CSR storage format. ! uirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the U factor in uval, according to the CSR storage format. + ! of the U factor in uval, according to the CSR storage format. ! l1 - integer, output ! The number of nonzero entries in laspk. ! l2 - integer, output ! The number of nonzero entries in uval. - ! info - integer, output. + ! info - integer, output. ! Error code. ! subroutine psb_siluk_factint(fill_in,ialg,a,b,& @@ -304,7 +304,7 @@ contains implicit none - ! Arguments + ! Arguments integer(psb_ipk_), intent(in) :: fill_in, ialg type(psb_sspmat_type),intent(in) :: a,b integer(psb_ipk_),intent(inout) :: l1,l2,info @@ -330,14 +330,14 @@ contains select case(ialg) case(psb_ilu_n_,psb_milu_n_) - ! Ok + ! Ok case default info=psb_err_input_asize_invalid_i_ call psb_errpush(info,name,& & i_err=(/itwo,ialg,izero,izero,izero/)) goto 9999 end select - if (fill_in < 0) then + if (fill_in < 0) then info=psb_err_input_asize_invalid_i_ call psb_errpush(info,name, & & i_err=(/ione,fill_in,izero,izero,izero/)) @@ -349,7 +349,7 @@ contains m = ma+mb ! - ! Allocate a temporary buffer for the iluk_copyin function + ! Allocate a temporary buffer for the iluk_copyin function ! call trw%allocate(izero,izero,ione) @@ -361,7 +361,7 @@ contains call psb_errpush(info,name,a_err='psb_sp_all') goto 9999 end if - + l1=0 l2=0 lirp(1) = 1 @@ -381,34 +381,34 @@ contains uplevs(:) = m+1 row(:) = szero rowlevs(:) = -(m+1) - + ! ! Cycle over the matrix rows ! do i = 1, m - + ! ! At each iteration of the loop we keep in a heap the column indices ! affected by the factorization. The heap is initialized and filled ! in the iluk_copyin routine, and updated during the elimination, in ! the iluk_fact routine. The heap is ideal because at each step we need ! the lowest index, but we also need to insert new items, and the heap - ! allows to do both in log time. + ! allows to do both in log time. ! d(i) = szero if (i<=ma) then ! - ! Copy into trw the i-th local row of the matrix, stored in a - ! + ! Copy into trw the i-th local row of the matrix, stored in a + ! call iluk_copyin(i,ma,a,ione,m,row,rowlevs,heap,ktrw,trw,info) else ! ! Copy into trw the i-th local row of the matrix, stored in b - ! (as (i-ma)-th row) - ! + ! (as (i-ma)-th row) + ! call iluk_copyin(i-ma,mb,b,ione,m,row,rowlevs,heap,ktrw,trw,info) endif - + ! Do an elimination step on the current row. It turns out we only ! need to keep track of fill levels for the upper triangle, hence we ! do not have a lowlevs variable. @@ -417,7 +417,7 @@ contains & d,uja,uirp,uval,uplevs,nidx,idxs,info) ! ! Copy the row into lval/d(i)/uval - ! + ! if (info == psb_success_) call iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,& & l1,l2,lja,lirp,lval,d,uja,uirp,uval,uplevs,info) if (info /= psb_success_) then @@ -474,11 +474,11 @@ contains ! ! This routine is used by psb_siluk_factint in the computation of the ! ILU(k)/MILU(k) factorization of a local sparse matrix. - ! + ! ! ! Arguments: ! i - integer, input. - ! The local index of the row to be extracted from the + ! The local index of the row to be extracted from the ! sparse matrix structure a. ! m - integer, input. ! The number of rows of the local matrix stored into a. @@ -510,7 +510,7 @@ contains ! staging buffer trw. See below. ! trw - type(psb_sspmat_type), input/output. ! A staging buffer. If the matrix A is not in CSR format, we use - ! the psb_sp_getblk routine and store its output in trw; when we + ! the psb_sp_getblk routine and store its output in trw; when we ! need to call psb_sp_getblk we do it for a block of rows, and then ! we consume them from trw in successive calls to this routine, ! until we empty the buffer. Thus we will make a call to psb_sp_getblk @@ -521,8 +521,8 @@ contains use psb_base_mod implicit none - - ! Arguments + + ! Arguments type(psb_sspmat_type), intent(in) :: a type(psb_s_coo_sparse_mat), intent(inout) :: trw integer(psb_ipk_), intent(in) :: i,m,jmin,jmax @@ -542,17 +542,17 @@ contains if (psb_errstatus_fatal()) then info = psb_err_internal_error_; goto 9999 end if - call heap%init(info) + call heap%init(info) - select type (aa=> a%a) - type is (psb_s_csr_sparse_mat) + select type (aa=> a%a) + type is (psb_s_csr_sparse_mat) ! ! Take a fast shortcut if the matrix is stored in CSR format ! - + do j = aa%irp(i), aa%irp(i+1) - 1 k = aa%ja(j) - if ((jmin<=k).and.(k<=jmax)) then + if ((jmin<=k).and.(k<=jmax)) then row(k) = aa%val(j) rowlevs(k) = 0 call heap%insert(k,info) @@ -562,14 +562,14 @@ contains class default ! - ! Otherwise use psb_sp_getblk, slower but able (in principle) of + ! Otherwise use psb_sp_getblk, slower but able (in principle) of ! handling any format. In this case, a block of rows is extracted ! instead of a single row, for performance reasons, and these ! rows are copied one by one into the array row, through successive ! calls to iluk_copyin. ! - if ((mod(i,nrb) == 1).or.(nrb == 1)) then + if ((mod(i,nrb) == 1).or.(nrb == 1)) then irb = min(m-i+1,nrb) call aa%csget(i,i+irb-1,trw,info) if (info /= psb_success_) then @@ -581,11 +581,11 @@ contains ktrw=1 end if nz = trw%get_nzeros() - do + do if (ktrw > nz) exit if (trw%ia(ktrw) > i) exit k = trw%ja(ktrw) - if ((jmin<=k).and.(k<=jmax)) then + if ((jmin<=k).and.(k<=jmax)) then row(k) = trw%val(ktrw) rowlevs(k) = 0 call heap%insert(k,info) @@ -674,10 +674,10 @@ contains use psb_base_mod - implicit none + implicit none ! Arguments - type(psb_i_heap), intent(inout) :: heap + type(psb_i_heap), intent(inout) :: heap integer(psb_ipk_), intent(in) :: i, fill_in integer(psb_ipk_), intent(inout) :: nidx,info integer(psb_ipk_), intent(inout) :: rowlevs(:) @@ -690,7 +690,7 @@ contains real(psb_spk_) :: rwk info = psb_success_ - if (.not.allocated(idxs)) then + if (.not.allocated(idxs)) then allocate(idxs(200),stat=info) if (info /= psb_success_) return endif @@ -702,34 +702,34 @@ contains ! do ! Beware: (iret < 0) means that the heap is empty, not an error. - call heap%get_first(k,iret) + call heap%get_first(k,iret) if (iret < 0) return - ! + ! ! Just in case an index has been put on the heap more than once. ! if (k == lastk) cycle - lastk = k + lastk = k nidx = nidx + 1 - if (nidx>size(idxs)) then + if (nidx>size(idxs)) then call psb_realloc(nidx+psb_heap_resize,idxs,info) if (info /= psb_success_) return end if idxs(nidx) = k - - if ((row(k) /= szero).and.(rowlevs(k) <= fill_in).and.(ki) then + else if (j>i) then ! ! Copy the upper part of the row - ! - if (rowlevs(j) <= fill_in) then - l2 = l2 + 1 - if (size(uval) < l2) then - ! + ! + if (rowlevs(j) <= fill_in) then + l2 = l2 + 1 + if (size(uval) < l2) then + ! ! Figure out a good reallocation size! ! isz = max((l2/i)*m,int(1.2*l2),l2+100) - call psb_realloc(isz,uval,info) - if (info == psb_success_) call psb_realloc(isz,uja,info) + call psb_realloc(isz,uval,info) + if (info == psb_success_) call psb_realloc(isz,uja,info) if (info == psb_success_) call psb_realloc(isz,uplevs,info,pad=(m+1)) - if (info /= psb_success_) then + if (info /= psb_success_) then info=psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') goto 9999 @@ -941,11 +941,11 @@ contains end if uja(l2) = j uval(l2) = row(j) - uplevs(l2) = rowlevs(j) + uplevs(l2) = rowlevs(j) else if (ialg == psb_milu_n_) then ! ! MILU(k): add discarded entries to the diagonal one - ! + ! d(i) = d(i) + row(j) end if ! @@ -963,13 +963,13 @@ contains lirp(i+1) = l1 + 1 uirp(i+1) = l2 + 1 - ! + ! ! Check the pivot size ! if (abs(d(i)) < s_epstol) then ! ! Too small pivot: unstable factorization - ! + ! info = psb_err_pivot_too_small_ int_err(1) = i write(ch_err,'(g20.10)') d(i) diff --git a/prec/impl/psb_s_ilut_fact.f90 b/prec/impl/psb_s_ilut_fact.f90 index 33b4374c..43cacf41 100644 --- a/prec/impl/psb_s_ilut_fact.f90 +++ b/prec/impl/psb_s_ilut_fact.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,21 +27,21 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! +! ! Moved here from MLD2P4, original copyright below. -! -! -! +! +! +! ! MLD2P4 version 2.2 ! MultiLevel Domain Decomposition Parallel Preconditioners Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2008-2018 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Daniela di Serafino -! +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -53,7 +53,7 @@ ! 3. The name of the MLD2P4 group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -65,8 +65,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: psb_silut_fact.f90 ! ! Subroutine: psb_silut_fact @@ -87,7 +87,7 @@ ! CSR. The diagonal of the U factor is stored separately (actually, the ! inverse of the diagonal entries is stored; this is then managed in the ! solve stage associated to the ILU(k,t) factorization). -! +! ! ! Arguments: ! fill_in - integer, input. @@ -114,7 +114,7 @@ ! factorization. ! Note: its allocation is managed by the calling routine psb_ilu_bld, ! hence it cannot be only intent(out). -! info - integer, output. +! info - integer, output. ! Error code. ! blck - type(psb_sspmat_type), input, optional, target. ! The sparse matrix structure containing the remote rows of the @@ -122,9 +122,9 @@ ! to build an Additive Schwarz base preconditioner with overlap ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (see psb_fact_bld), then blck does not contain any row. -! +! subroutine psb_silut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) - + use psb_base_mod use psb_s_ilu_fact_mod, psb_protect_name => psb_silut_fact @@ -141,7 +141,7 @@ subroutine psb_silut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) integer(psb_ipk_), intent(in), optional :: iscale ! Local Variables integer(psb_ipk_) :: l1, l2, m, err_act, iscale_ - + type(psb_sspmat_type), pointer :: blck_ type(psb_s_csr_sparse_mat) :: ll, uu real(psb_spk_) :: scale @@ -151,34 +151,34 @@ subroutine psb_silut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) info = psb_success_ call psb_erractionsave(err_act) - if (fill_in < 0) then + if (fill_in < 0) then info=psb_err_input_asize_invalid_i_ call psb_errpush(info,name, & & i_err=(/ione,fill_in,izero,izero,izero/)) goto 9999 end if - ! + ! ! Point to / allocate memory for the incomplete factorization ! - if (present(blck)) then + if (present(blck)) then blck_ => blck else - allocate(blck_,stat=info) - if (info == psb_success_) call blck_%csall(izero,izero,info,ione) + allocate(blck_,stat=info) + if (info == psb_success_) call blck_%allocate(izero,izero,info,ione,type='CSR') if (info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='csall' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if endif - if (present(iscale)) then + if (present(iscale)) then iscale_ = iscale else iscale_ = psb_ilu_scale_none_ end if - - select case(iscale_) + + select case(iscale_) case(psb_ilu_scale_none_) scale = sone case(psb_ilu_scale_maxval_) @@ -189,10 +189,10 @@ subroutine psb_silut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) call psb_errpush(info,name,i_err=(/ione*9,iscale_,izero,izero,izero/)) goto 9999 end select - + m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& - & (m > size(d)) ) then + & (m > size(d)) ) then write(0,*) 'Wrong allocation status for L,D,U? ',& & l%get_nrows(),size(d),u%get_nrows() info = -1 @@ -225,15 +225,15 @@ subroutine psb_silut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) call u%set_triangle() call u%set_unit() call u%set_upper() - + ! ! Nullify pointer / deallocate memory ! - if (present(blck)) then - blck_ => null() + if (present(blck)) then + blck_ => null() else call blck_%free() - deallocate(blck_,stat=info) + deallocate(blck_,stat=info) if(info.ne.0) then info=psb_err_from_subroutine_ ch_err='psb_sp_free' @@ -298,7 +298,7 @@ contains ! according to the CSR storage format. ! lirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the L factor in lval, according to the CSR storage format. + ! of the L factor in lval, according to the CSR storage format. ! uval - real(psb_spk_), dimension(:), input/output. ! The U factor in the incomplete factorization. ! The entries of U are stored according to the CSR format. @@ -307,12 +307,12 @@ contains ! according to the CSR storage format. ! uirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the U factor in uval, according to the CSR storage format. + ! of the U factor in uval, according to the CSR storage format. ! l1 - integer, output ! The number of nonzero entries in lval. ! l2 - integer, output ! The number of nonzero entries in uval. - ! info - integer, output. + ! info - integer, output. ! Error code. ! subroutine psb_silut_factint(fill_in,thres,a,b,& @@ -320,7 +320,7 @@ contains use psb_base_mod - implicit none + implicit none ! Arguments integer(psb_ipk_), intent(in) :: fill_in @@ -355,7 +355,7 @@ contains m = ma+mb ! - ! Allocate a temporary buffer for the ilut_copyin function + ! Allocate a temporary buffer for the ilut_copyin function ! call trw%allocate(izero,izero,ione) if (info == psb_success_) call psb_ensure_size(m+1,lirp,info) @@ -366,7 +366,7 @@ contains call psb_errpush(info,name,a_err='psb_sp_all') goto 9999 end if - + l1=0 l2=0 lirp(1) = 1 @@ -396,10 +396,10 @@ contains ! in the ilut_copyin function, and updated during the elimination, in ! the ilut_fact routine. The heap is ideal because at each step we need ! the lowest index, but we also need to insert new items, and the heap - ! allows to do both in log time. + ! allows to do both in log time. ! d(i) = czero - if (i<=ma) then + if (i<=ma) then call ilut_copyin(i,ma,a,i,ione,m,nlw,nup,jmaxup,nrmi,weight,& & row,heap,ktrw,trw,info) else @@ -414,7 +414,7 @@ contains & d,uja,uirp,uval,nidx,idxs,info) ! ! Copy the row into lval/d(i)/uval - ! + ! if (info == psb_success_) call ilut_copyout(fill_in,thres,i,m,& & nlw,nup,jmaxup,nrmi,row,nidx,idxs,& & l1,l2,lja,lirp,lval,d,uja,uirp,uval,info) @@ -429,7 +429,7 @@ contains ! ! Adjust diagonal accounting for scale factor ! - if (weight /= sone) then + if (weight /= sone) then d(1:m) = d(1:m)*weight end if @@ -469,7 +469,7 @@ contains ! - storing into a heap the column indices of the nonzero entries of the copied ! row; ! - computing the column index of the first entry with maximum absolute value - ! in the part of the row belonging to the upper triangle; + ! in the part of the row belonging to the upper triangle; ! - computing the 2-norm of the row. ! The output array row is such that it contains a full row of A, i.e. it contains ! also the zero entries of the row. This is useful for the elimination step @@ -482,11 +482,11 @@ contains ! ! This routine is used by psb_silut_factint in the computation of the ILU(k,t) ! factorization of a local sparse matrix. - ! + ! ! ! Arguments: ! i - integer, input. - ! The local index of the row to be extracted from the + ! The local index of the row to be extracted from the ! sparse matrix structure a. ! m - integer, input. ! The number of rows of the local matrix stored into a. @@ -533,7 +533,7 @@ contains ! staging buffer trw. See below. ! trw - type(psb_sspmat_type), input/output. ! A staging buffer. If the matrix A is not in CSR format, we use - ! the psb_sp_getblk routine and store its output in trw; when we + ! the psb_sp_getblk routine and store its output in trw; when we ! need to call psb_sp_getblk we do it for a block of rows, and then ! we consume them from trw in successive calls to this routine, ! until we empty the buffer. Thus we will make a call to psb_sp_getblk @@ -542,7 +542,7 @@ contains subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,& & nrmi,weight,row,heap,ktrw,trw,info) use psb_base_mod - implicit none + implicit none type(psb_sspmat_type), intent(in) :: a type(psb_s_coo_sparse_mat), intent(inout) :: trw integer(psb_ipk_), intent(in) :: i, m,jmin,jmax,jd @@ -584,29 +584,29 @@ contains dmaxup = szero nrmi = szero - select type (aa=> a%a) - type is (psb_s_csr_sparse_mat) + select type (aa=> a%a) + type is (psb_s_csr_sparse_mat) ! ! Take a fast shortcut if the matrix is stored in CSR format ! do j = aa%irp(i), aa%irp(i+1) - 1 k = aa%ja(j) - if ((jmin<=k).and.(k<=jmax)) then - row(k) = aa%val(j)*weight + if ((jmin<=k).and.(k<=jmax)) then + row(k) = aa%val(j)*weight call heap%insert(k,info) if (info /= psb_success_) exit - if (kjd) then + if (kjd) then nup = nup + 1 - if (abs(row(k))>dmaxup) then + if (abs(row(k))>dmaxup) then jmaxup = k dmaxup = abs(row(k)) end if end if end if end do - if (info /= psb_success_) then + if (info /= psb_success_) then info=psb_err_from_subroutine_ call psb_errpush(info,name,a_err='psb_insert_heap') goto 9999 @@ -617,16 +617,16 @@ contains class default - + ! - ! Otherwise use psb_sp_getblk, slower but able (in principle) of + ! Otherwise use psb_sp_getblk, slower but able (in principle) of ! handling any format. In this case, a block of rows is extracted ! instead of a single row, for performance reasons, and these ! rows are copied one by one into the array row, through successive ! calls to ilut_copyin. ! - if ((mod(i,nrb) == 1).or.(nrb == 1)) then + if ((mod(i,nrb) == 1).or.(nrb == 1)) then irb = min(m-i+1,nrb) call aa%csget(i,i+irb-1,trw,info) if (info /= psb_success_) then @@ -639,18 +639,18 @@ contains kin = ktrw nz = trw%get_nzeros() - do + do if (ktrw > nz) exit if (trw%ia(ktrw) > i) exit k = trw%ja(ktrw) - if ((jmin<=k).and.(k<=jmax)) then + if ((jmin<=k).and.(k<=jmax)) then row(k) = trw%val(ktrw)*weight call heap%insert(k,info) if (info /= psb_success_) exit - if (kjd) then + if (kjd) then nup = nup + 1 - if (abs(row(k))>dmaxup) then + if (abs(row(k))>dmaxup) then jmaxup = k dmaxup = abs(row(k)) end if @@ -734,10 +734,10 @@ contains use psb_base_mod - implicit none + implicit none ! Arguments - type(psb_i_heap), intent(inout) :: heap + type(psb_i_heap), intent(inout) :: heap integer(psb_ipk_), intent(in) :: i integer(psb_ipk_), intent(inout) :: nidx,info real(psb_spk_), intent(in) :: thres,nrmi @@ -753,23 +753,23 @@ contains call psb_ensure_size(200*ione,idxs,info) if (info /= psb_success_) return nidx = 0 - lastk = -1 + lastk = -1 ! ! Do while there are indices to be processed ! do - call heap%get_first(k,iret) + call heap%get_first(k,iret) if (iret < 0) exit - ! + ! ! An index may have been put on the heap more than once. ! if (k == lastk) cycle - lastk = k + lastk = k lowert: if (k nidx) exit if (idxs(idxp) >= i) exit @@ -987,8 +987,8 @@ contains ! if (abs(witem) < thres*nrmi) cycle - nz = nz + 1 - xw(nz) = witem + nz = nz + 1 + xw(nz) = witem xwid(nz) = widx call heap%insert(witem,widx,info) if (info /= psb_success_) then @@ -1000,9 +1000,9 @@ contains ! ! Now we have to take out the first nlw+fill_in entries - ! + ! if (nz <= nlw+fill_in) then - ! + ! ! Just copy everything from xw, and it is already ordered ! else @@ -1014,7 +1014,7 @@ contains call psb_errpush(info,name,a_err='psb_heap_get_first') goto 9999 end if - + xw(k) = witem xwid(k) = widx end do @@ -1029,15 +1029,15 @@ contains ! Copy out the lower part of the row ! do k=1,nz - l1 = l1 + 1 + l1 = l1 + 1 if (size(lval) < l1) then - ! + ! ! Figure out a good reallocation size! - ! + ! isz = (max((l1/i)*m,int(1.2*l1),l1+100)) - call psb_realloc(isz,lval,info) - if (info == psb_success_) call psb_realloc(isz,lja,info) - if (info /= psb_success_) then + call psb_realloc(isz,lval,info) + if (info == psb_success_) call psb_realloc(isz,lja,info) + if (info /= psb_success_) then info=psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') goto 9999 @@ -1046,25 +1046,25 @@ contains lja(l1) = xwid(k) lval(l1) = xw(indx(k)) end do - + ! ! Make sure idxp points to the diagonal entry ! - if (idxp <= size(idxs)) then - if (idxs(idxp) < i) then - do + if (idxp <= size(idxs)) then + if (idxs(idxp) < i) then + do idxp = idxp + 1 if (idxp > nidx) exit if (idxs(idxp) >= i) exit end do end if end if - if (idxp > size(idxs)) then + if (idxp > size(idxs)) then !!$ write(0,*) 'Warning: missing diagonal element in the row ' else - if (idxs(idxp) > i) then + if (idxs(idxp) > i) then !!$ write(0,*) 'Warning: missing diagonal element in the row ' - else if (idxs(idxp) /= i) then + else if (idxs(idxp) /= i) then !!$ write(0,*) 'Warning: impossible error: diagonal has vanished' else ! @@ -1076,7 +1076,7 @@ contains if (abs(d(i)) < s_epstol) then ! ! Too small pivot: unstable factorization - ! + ! info = psb_err_pivot_too_small_ int_err(1) = i write(ch_err,'(g20.10)') d(i) @@ -1092,7 +1092,7 @@ contains end if ! - ! Now the upper part + ! Now the upper part ! call heap%init(info,dir=psb_asort_down_) @@ -1104,15 +1104,15 @@ contains nz = 0 do - + idxp = idxp + 1 if (idxp > nidx) exit widx = idxs(idxp) - if (widx <= i) then + if (widx <= i) then !!$ write(0,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) cycle end if - if (widx > m) then + if (widx > m) then !!$ write(0,*) 'Warning: impossible value',widx,i,idxp,idxs(idxp) cycle end if @@ -1120,12 +1120,12 @@ contains ! ! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway. ! - if ((widx /= jmaxup) .and. (abs(witem) < thres*nrmi)) then - cycle + if ((widx /= jmaxup) .and. (abs(witem) < thres*nrmi)) then + cycle end if nz = nz + 1 - xw(nz) = witem + xw(nz) = witem xwid(nz) = widx call heap%insert(witem,widx,info) if (info /= psb_success_) then @@ -1133,7 +1133,7 @@ contains call psb_errpush(info,name,a_err='psb_insert_heap') goto 9999 end if - + end do ! @@ -1141,7 +1141,7 @@ contains ! we include entry jmaxup. ! if (nz <= nup+fill_in) then - ! + ! ! Just copy everything from xw ! fndmaxup=.true. @@ -1155,13 +1155,13 @@ contains if (widx == jmaxup) fndmaxup=.true. end do end if - if ((i blck else - allocate(blck_,stat=info) - if (info == psb_success_) call blck_%csall(izero,izero,info,ione) + allocate(blck_,stat=info) + if (info == psb_success_) call blck_%allocate(izero,izero,info,ione,type='CSR') if (info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='csall' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if endif - if (present(upd)) then - upd_ = psb_toupper(upd) + if (present(upd)) then + upd_ = psb_toupper(upd) else upd_ = 'F' end if - + m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& - & (m > size(d)) ) then + & (m > size(d)) ) then write(0,*) 'Wrong allocation status for L,D,U? ',& & l%get_nrows(),size(d),u%get_nrows() info = -1 @@ -212,12 +212,12 @@ subroutine psb_zilu0_fact(ialg,a,l,u,d,info,blck, upd) call u%set_triangle() call u%set_unit() call u%set_upper() - + ! ! Nullify pointer / deallocate memory ! - if (present(blck)) then - blck_ => null() + if (present(blck)) then + blck_ => null() else call blck_%free() if(info.ne.0) then @@ -226,7 +226,7 @@ subroutine psb_zilu0_fact(ialg,a,l,u,d,info,blck, upd) call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - deallocate(blck_) + deallocate(blck_) endif call psb_erractionrestore(err_act) @@ -243,9 +243,9 @@ contains ! Version: complex ! Note: internal subroutine of psb_zilu0_fact. ! - ! This routine computes either the ILU(0) or the MILU(0) factorization of the - ! diagonal blocks of a distributed matrix. - ! These factorizations are used to build the 'base preconditioner' + ! This routine computes either the ILU(0) or the MILU(0) factorization of the + ! diagonal blocks of a distributed matrix. + ! These factorizations are used to build the 'base preconditioner' ! (block-Jacobi preconditioner/solver, Additive Schwarz ! preconditioner) corresponding to a given level of a multilevel preconditioner. ! @@ -257,7 +257,7 @@ contains ! ! The routine copies and factors "on the fly" from the sparse matrix structures a ! and b into the arrays lval, uval, d (L, U without its diagonal, diagonal of U). - ! + ! ! ! Arguments: ! ialg - integer, input. @@ -296,7 +296,7 @@ contains ! according to the CSR storage format. ! lirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the L factor in lval, according to the CSR storage format. + ! of the L factor in lval, according to the CSR storage format. ! uval - complex(psb_dpk_), dimension(:), input/output. ! The U factor in the incomplete factorization. ! The entries of U are stored according to the CSR format. @@ -305,18 +305,18 @@ contains ! according to the CSR storage format. ! uirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the U factor in uval, according to the CSR storage format. + ! of the U factor in uval, according to the CSR storage format. ! l1 - integer, output. ! The number of nonzero entries in lval. ! l2 - integer, output. ! The number of nonzero entries in uval. - ! info - integer, output. + ! info - integer, output. ! Error code. ! subroutine psb_zilu0_factint(ialg,a,b,& & d,lval,lja,lirp,uval,uja,uirp,l1,l2,upd,info) - implicit none + implicit none ! Arguments integer(psb_ipk_), intent(in) :: ialg @@ -332,7 +332,7 @@ contains complex(psb_dpk_) :: dia,temp integer(psb_ipk_), parameter :: nrb=16 type(psb_z_coo_sparse_mat) :: trw - integer(psb_ipk_) :: int_err(5) + integer(psb_ipk_) :: int_err(5) character(len=20) :: name, ch_err name='psb_zilu0_factint' @@ -346,10 +346,10 @@ contains select case(ialg) case(psb_ilu_n_,psb_milu_n_) - ! Ok + ! Ok case default info=psb_err_input_asize_invalid_i_ - call psb_errpush(info,name,& + call psb_errpush(info,name,& & i_err=(/ione,ialg,izero,izero,izero/)) goto 9999 end select @@ -363,7 +363,7 @@ contains end if m = ma+mb - if (psb_toupper(upd) == 'F' ) then + if (psb_toupper(upd) == 'F' ) then lirp(1) = 1 uirp(1) = 1 l1 = 0 @@ -379,14 +379,14 @@ contains if (i <= ma) then ! ! Copy the i-th local row of the matrix, stored in a, - ! into lval/d(i)/uval + ! into lval/d(i)/uval ! call ilu_copyin(i,ma,a,i,ione,m,l1,lja,lval,& & d(i),l2,uja,uval,ktrw,trw,upd) else ! ! Copy the i-th local row of the matrix, stored in b - ! (as (i-ma)-th row), into lval/d(i)/uval + ! (as (i-ma)-th row), into lval/d(i)/uval ! call ilu_copyin(i-ma,mb,b,i,ione,m,l1,lja,lval,& & d(i),l2,uja,uval,ktrw,trw,upd) @@ -437,7 +437,7 @@ contains ! dia = dia - temp*uval(jj) cycle updateloop - ! + ! else if (j > i) then ! ! search u(i,*) (i-th row of U) for a matching index j @@ -454,23 +454,23 @@ contains end if enddo end if - ! - ! If we get here we missed the cycle updateloop, which means + ! + ! If we get here we missed the cycle updateloop, which means ! that this entry does not match; thus we accumulate on the ! diagonal for MILU(0). ! - if (ialg == psb_milu_n_) then + if (ialg == psb_milu_n_) then dia = dia - temp*uval(jj) end if enddo updateloop enddo - ! + ! ! Check the pivot size - ! + ! if (abs(dia) < d_epstol) then ! ! Too small pivot: unstable factorization - ! + ! info = psb_err_pivot_too_small_ int_err(1) = i write(ch_err,'(g20.10)') abs(dia) @@ -493,7 +493,7 @@ contains else write(0,*) 'Update not implemented ' info = 31 - call psb_errpush(info,name,& + call psb_errpush(info,name,& & i_err=(/ione*13,izero,izero,izero,izero/),a_err=upd) goto 9999 @@ -515,7 +515,7 @@ contains ! Version: complex ! Note: internal subroutine of psb_zilu0_fact ! - ! This routine copies a row of a sparse matrix A, stored in the psb_zspmat_type + ! This routine copies a row of a sparse matrix A, stored in the psb_zspmat_type ! data structure a, into the arrays lval and uval and into the scalar variable ! dia, corresponding to the lower and upper triangles of A and to the diagonal ! entry of the row, respectively. The entries in lval and uval are stored @@ -529,14 +529,14 @@ contains ! ! The routine is used by psb_zilu0_factint in the computation of the ILU(0)/MILU(0) ! factorization of a local sparse matrix. - ! + ! ! TODO: modify the routine to allow copying into output L and U that are ! already filled with indices; this would allow computing an ILU(k) pattern, - ! then use the ILU(0) internal for subsequent calls with the same pattern. + ! then use the ILU(0) internal for subsequent calls with the same pattern. ! ! Arguments: ! i - integer, input. - ! The local index of the row to be extracted from the + ! The local index of the row to be extracted from the ! sparse matrix structure a. ! m - integer, input. ! The number of rows of the local matrix stored into a. @@ -570,13 +570,13 @@ contains ! to the CSR storage format. ! uval - complex(psb_dpk_), dimension(:), input/output. ! The array where the entries of the row corresponding to the - ! upper triangle are copied. + ! upper triangle are copied. ! ktrw - integer, input/output. ! The index identifying the last entry taken from the ! staging buffer trw. See below. ! trw - type(psb_zspmat_type), input/output. ! A staging buffer. If the matrix A is not in CSR format, we use - ! the psb_sp_getblk routine and store its output in trw; when we + ! the psb_sp_getblk routine and store its output in trw; when we ! need to call psb_sp_getblk we do it for a block of rows, and then ! we consume them from trw in successive calls to this routine, ! until we empty the buffer. Thus we will make a call to psb_sp_getblk @@ -608,10 +608,10 @@ contains if (psb_errstatus_fatal()) then info = psb_err_internal_error_; goto 9999 end if - if (psb_toupper(upd) == 'F') then + if (psb_toupper(upd) == 'F') then - select type(aa => a%a) - type is (psb_z_csr_sparse_mat) + select type(aa => a%a) + type is (psb_z_csr_sparse_mat) ! ! Take a fast shortcut if the matrix is stored in CSR format @@ -633,19 +633,19 @@ contains end if enddo - class default + class default ! - ! Otherwise use psb_sp_getblk, slower but able (in principle) of + ! Otherwise use psb_sp_getblk, slower but able (in principle) of ! handling any format. In this case, a block of rows is extracted ! instead of a single row, for performance reasons, and these ! rows are copied one by one into lval, dia, uval, through ! successive calls to ilu_copyin. ! - if ((mod(i,nrb) == 1).or.(nrb == 1)) then + if ((mod(i,nrb) == 1).or.(nrb == 1)) then irb = min(m-i+1,nrb) - call aa%csget(i,i+irb-1,trw,info) + call aa%csget(i,i+irb-1,trw,info) if(info /= psb_success_) then info=psb_err_from_subroutine_ ch_err='csget' @@ -656,7 +656,7 @@ contains end if nz = trw%get_nzeros() - do + do if (ktrw > nz) exit if (trw%ia(ktrw) > i) exit k = trw%ja(ktrw) @@ -680,7 +680,7 @@ contains write(0,*) 'Update not implemented ' info = 31 - call psb_errpush(info,name,& + call psb_errpush(info,name,& & i_err=(/ione*13,izero,izero,izero,izero/),a_err=upd) goto 9999 diff --git a/prec/impl/psb_z_iluk_fact.f90 b/prec/impl/psb_z_iluk_fact.f90 index fe9e92d9..1a398cda 100644 --- a/prec/impl/psb_z_iluk_fact.f90 +++ b/prec/impl/psb_z_iluk_fact.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,21 +27,21 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! +! ! Moved here from MLD2P4, original copyright below. -! -! -! +! +! +! ! MLD2P4 version 2.2 ! MultiLevel Domain Decomposition Parallel Preconditioners Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2008-2018 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Daniela di Serafino -! +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -53,7 +53,7 @@ ! 3. The name of the MLD2P4 group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -65,8 +65,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: psb_ziluk_fact.f90 ! ! Subroutine: psb_ziluk_fact @@ -74,7 +74,7 @@ ! Contains: psb_ziluk_factint, iluk_copyin, iluk_fact, iluk_copyout. ! ! This routine computes either the ILU(k) or the MILU(k) factorization of the -! diagonal blocks of a distributed matrix. These factorizations are used to +! diagonal blocks of a distributed matrix. These factorizations are used to ! build the 'base preconditioner' (block-Jacobi preconditioner/solver, ! Additive Schwarz preconditioner) corresponding to a certain level of a ! multilevel preconditioner. @@ -88,7 +88,7 @@ ! U factors is CSR. The diagonal of the U factor is stored separately (actually, ! the inverse of the diagonal entries is stored; this is then managed in the solve ! stage associated to the ILU(k)/MILU(k) factorization). -! +! ! ! Arguments: ! fill_in - integer, input. @@ -118,7 +118,7 @@ ! factorization. ! Note: its allocation is managed by the calling routine psb_ilu_bld, ! hence it cannot be only intent(out). -! info - integer, output. +! info - integer, output. ! Error code. ! blck - type(psb_zspmat_type), input, optional, target. ! The sparse matrix structure containing the remote rows of the @@ -126,7 +126,7 @@ ! to build an Additive Schwarz base preconditioner with overlap ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (see psb_fact_bld), then blck does not contain any row. -! +! subroutine psb_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck) use psb_base_mod @@ -143,7 +143,7 @@ subroutine psb_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck) complex(psb_dpk_), intent(inout) :: d(:) ! Local Variables integer(psb_ipk_) :: l1, l2, m, err_act - + type(psb_zspmat_type), pointer :: blck_ type(psb_z_csr_sparse_mat) :: ll, uu character(len=20) :: name, ch_err @@ -152,17 +152,17 @@ subroutine psb_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck) info = psb_success_ call psb_erractionsave(err_act) - ! + ! ! Point to / allocate memory for the incomplete factorization ! - if (present(blck)) then + if (present(blck)) then blck_ => blck else - allocate(blck_,stat=info) - if (info == psb_success_) call blck_%csall(izero,izero,info,ione) + allocate(blck_,stat=info) + if (info == psb_success_) call blck_%allocate(izero,izero,info,ione,type='CSR') if (info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='csall' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if @@ -170,7 +170,7 @@ subroutine psb_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck) m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& - & (m > size(d)) ) then + & (m > size(d)) ) then write(0,*) 'Wrong allocation status for L,D,U? ',& & l%get_nrows(),size(d),u%get_nrows() info = -1 @@ -203,15 +203,15 @@ subroutine psb_ziluk_fact(fill_in,ialg,a,l,u,d,info,blck) call u%set_triangle() call u%set_unit() call u%set_upper() - + ! ! Nullify pointer / deallocate memory ! - if (present(blck)) then - blck_ => null() + if (present(blck)) then + blck_ => null() else call blck_%free() - deallocate(blck_,stat=info) + deallocate(blck_,stat=info) if(info.ne.0) then info=psb_err_from_subroutine_ ch_err='psb_sp_free' @@ -280,7 +280,7 @@ contains ! according to the CSR storage format. ! lia2 - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the L factor in laspk, according to the CSR storage format. + ! of the L factor in laspk, according to the CSR storage format. ! uval - complex(psb_dpk_), dimension(:), input/output. ! The U factor in the incomplete factorization. ! The entries of U are stored according to the CSR format. @@ -289,12 +289,12 @@ contains ! according to the CSR storage format. ! uirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the U factor in uval, according to the CSR storage format. + ! of the U factor in uval, according to the CSR storage format. ! l1 - integer, output ! The number of nonzero entries in laspk. ! l2 - integer, output ! The number of nonzero entries in uval. - ! info - integer, output. + ! info - integer, output. ! Error code. ! subroutine psb_ziluk_factint(fill_in,ialg,a,b,& @@ -304,7 +304,7 @@ contains implicit none - ! Arguments + ! Arguments integer(psb_ipk_), intent(in) :: fill_in, ialg type(psb_zspmat_type),intent(in) :: a,b integer(psb_ipk_),intent(inout) :: l1,l2,info @@ -330,14 +330,14 @@ contains select case(ialg) case(psb_ilu_n_,psb_milu_n_) - ! Ok + ! Ok case default info=psb_err_input_asize_invalid_i_ call psb_errpush(info,name,& & i_err=(/itwo,ialg,izero,izero,izero/)) goto 9999 end select - if (fill_in < 0) then + if (fill_in < 0) then info=psb_err_input_asize_invalid_i_ call psb_errpush(info,name, & & i_err=(/ione,fill_in,izero,izero,izero/)) @@ -349,7 +349,7 @@ contains m = ma+mb ! - ! Allocate a temporary buffer for the iluk_copyin function + ! Allocate a temporary buffer for the iluk_copyin function ! call trw%allocate(izero,izero,ione) @@ -361,7 +361,7 @@ contains call psb_errpush(info,name,a_err='psb_sp_all') goto 9999 end if - + l1=0 l2=0 lirp(1) = 1 @@ -381,34 +381,34 @@ contains uplevs(:) = m+1 row(:) = zzero rowlevs(:) = -(m+1) - + ! ! Cycle over the matrix rows ! do i = 1, m - + ! ! At each iteration of the loop we keep in a heap the column indices ! affected by the factorization. The heap is initialized and filled ! in the iluk_copyin routine, and updated during the elimination, in ! the iluk_fact routine. The heap is ideal because at each step we need ! the lowest index, but we also need to insert new items, and the heap - ! allows to do both in log time. + ! allows to do both in log time. ! d(i) = zzero if (i<=ma) then ! - ! Copy into trw the i-th local row of the matrix, stored in a - ! + ! Copy into trw the i-th local row of the matrix, stored in a + ! call iluk_copyin(i,ma,a,ione,m,row,rowlevs,heap,ktrw,trw,info) else ! ! Copy into trw the i-th local row of the matrix, stored in b - ! (as (i-ma)-th row) - ! + ! (as (i-ma)-th row) + ! call iluk_copyin(i-ma,mb,b,ione,m,row,rowlevs,heap,ktrw,trw,info) endif - + ! Do an elimination step on the current row. It turns out we only ! need to keep track of fill levels for the upper triangle, hence we ! do not have a lowlevs variable. @@ -417,7 +417,7 @@ contains & d,uja,uirp,uval,uplevs,nidx,idxs,info) ! ! Copy the row into lval/d(i)/uval - ! + ! if (info == psb_success_) call iluk_copyout(fill_in,ialg,i,m,row,rowlevs,nidx,idxs,& & l1,l2,lja,lirp,lval,d,uja,uirp,uval,uplevs,info) if (info /= psb_success_) then @@ -474,11 +474,11 @@ contains ! ! This routine is used by psb_ziluk_factint in the computation of the ! ILU(k)/MILU(k) factorization of a local sparse matrix. - ! + ! ! ! Arguments: ! i - integer, input. - ! The local index of the row to be extracted from the + ! The local index of the row to be extracted from the ! sparse matrix structure a. ! m - integer, input. ! The number of rows of the local matrix stored into a. @@ -510,7 +510,7 @@ contains ! staging buffer trw. See below. ! trw - type(psb_zspmat_type), input/output. ! A staging buffer. If the matrix A is not in CSR format, we use - ! the psb_sp_getblk routine and store its output in trw; when we + ! the psb_sp_getblk routine and store its output in trw; when we ! need to call psb_sp_getblk we do it for a block of rows, and then ! we consume them from trw in successive calls to this routine, ! until we empty the buffer. Thus we will make a call to psb_sp_getblk @@ -521,8 +521,8 @@ contains use psb_base_mod implicit none - - ! Arguments + + ! Arguments type(psb_zspmat_type), intent(in) :: a type(psb_z_coo_sparse_mat), intent(inout) :: trw integer(psb_ipk_), intent(in) :: i,m,jmin,jmax @@ -542,17 +542,17 @@ contains if (psb_errstatus_fatal()) then info = psb_err_internal_error_; goto 9999 end if - call heap%init(info) + call heap%init(info) - select type (aa=> a%a) - type is (psb_z_csr_sparse_mat) + select type (aa=> a%a) + type is (psb_z_csr_sparse_mat) ! ! Take a fast shortcut if the matrix is stored in CSR format ! - + do j = aa%irp(i), aa%irp(i+1) - 1 k = aa%ja(j) - if ((jmin<=k).and.(k<=jmax)) then + if ((jmin<=k).and.(k<=jmax)) then row(k) = aa%val(j) rowlevs(k) = 0 call heap%insert(k,info) @@ -562,14 +562,14 @@ contains class default ! - ! Otherwise use psb_sp_getblk, slower but able (in principle) of + ! Otherwise use psb_sp_getblk, slower but able (in principle) of ! handling any format. In this case, a block of rows is extracted ! instead of a single row, for performance reasons, and these ! rows are copied one by one into the array row, through successive ! calls to iluk_copyin. ! - if ((mod(i,nrb) == 1).or.(nrb == 1)) then + if ((mod(i,nrb) == 1).or.(nrb == 1)) then irb = min(m-i+1,nrb) call aa%csget(i,i+irb-1,trw,info) if (info /= psb_success_) then @@ -581,11 +581,11 @@ contains ktrw=1 end if nz = trw%get_nzeros() - do + do if (ktrw > nz) exit if (trw%ia(ktrw) > i) exit k = trw%ja(ktrw) - if ((jmin<=k).and.(k<=jmax)) then + if ((jmin<=k).and.(k<=jmax)) then row(k) = trw%val(ktrw) rowlevs(k) = 0 call heap%insert(k,info) @@ -674,10 +674,10 @@ contains use psb_base_mod - implicit none + implicit none ! Arguments - type(psb_i_heap), intent(inout) :: heap + type(psb_i_heap), intent(inout) :: heap integer(psb_ipk_), intent(in) :: i, fill_in integer(psb_ipk_), intent(inout) :: nidx,info integer(psb_ipk_), intent(inout) :: rowlevs(:) @@ -690,7 +690,7 @@ contains complex(psb_dpk_) :: rwk info = psb_success_ - if (.not.allocated(idxs)) then + if (.not.allocated(idxs)) then allocate(idxs(200),stat=info) if (info /= psb_success_) return endif @@ -702,34 +702,34 @@ contains ! do ! Beware: (iret < 0) means that the heap is empty, not an error. - call heap%get_first(k,iret) + call heap%get_first(k,iret) if (iret < 0) return - ! + ! ! Just in case an index has been put on the heap more than once. ! if (k == lastk) cycle - lastk = k + lastk = k nidx = nidx + 1 - if (nidx>size(idxs)) then + if (nidx>size(idxs)) then call psb_realloc(nidx+psb_heap_resize,idxs,info) if (info /= psb_success_) return end if idxs(nidx) = k - - if ((row(k) /= zzero).and.(rowlevs(k) <= fill_in).and.(ki) then + else if (j>i) then ! ! Copy the upper part of the row - ! - if (rowlevs(j) <= fill_in) then - l2 = l2 + 1 - if (size(uval) < l2) then - ! + ! + if (rowlevs(j) <= fill_in) then + l2 = l2 + 1 + if (size(uval) < l2) then + ! ! Figure out a good reallocation size! ! isz = max((l2/i)*m,int(1.2*l2),l2+100) - call psb_realloc(isz,uval,info) - if (info == psb_success_) call psb_realloc(isz,uja,info) + call psb_realloc(isz,uval,info) + if (info == psb_success_) call psb_realloc(isz,uja,info) if (info == psb_success_) call psb_realloc(isz,uplevs,info,pad=(m+1)) - if (info /= psb_success_) then + if (info /= psb_success_) then info=psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') goto 9999 @@ -941,11 +941,11 @@ contains end if uja(l2) = j uval(l2) = row(j) - uplevs(l2) = rowlevs(j) + uplevs(l2) = rowlevs(j) else if (ialg == psb_milu_n_) then ! ! MILU(k): add discarded entries to the diagonal one - ! + ! d(i) = d(i) + row(j) end if ! @@ -963,13 +963,13 @@ contains lirp(i+1) = l1 + 1 uirp(i+1) = l2 + 1 - ! + ! ! Check the pivot size ! if (abs(d(i)) < d_epstol) then ! ! Too small pivot: unstable factorization - ! + ! info = psb_err_pivot_too_small_ int_err(1) = i write(ch_err,'(g20.10)') d(i) diff --git a/prec/impl/psb_z_ilut_fact.f90 b/prec/impl/psb_z_ilut_fact.f90 index b7e8da05..291dc778 100644 --- a/prec/impl/psb_z_ilut_fact.f90 +++ b/prec/impl/psb_z_ilut_fact.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -27,21 +27,21 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! +! ! Moved here from MLD2P4, original copyright below. -! -! -! +! +! +! ! MLD2P4 version 2.2 ! MultiLevel Domain Decomposition Parallel Preconditioners Package ! based on PSBLAS (Parallel Sparse BLAS version 3.5) -! -! (C) Copyright 2008-2018 -! -! Salvatore Filippone -! Pasqua D'Ambra -! Daniela di Serafino -! +! +! (C) Copyright 2008-2018 +! +! Salvatore Filippone +! Pasqua D'Ambra +! Daniela di Serafino +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -53,7 +53,7 @@ ! 3. The name of the MLD2P4 group or the names of its contributors may ! not be used to endorse or promote products derived from this ! software without specific written permission. -! +! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED ! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR @@ -65,8 +65,8 @@ ! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE ! POSSIBILITY OF SUCH DAMAGE. -! -! +! +! ! File: psb_zilut_fact.f90 ! ! Subroutine: psb_zilut_fact @@ -87,7 +87,7 @@ ! CSR. The diagonal of the U factor is stored separately (actually, the ! inverse of the diagonal entries is stored; this is then managed in the ! solve stage associated to the ILU(k,t) factorization). -! +! ! ! Arguments: ! fill_in - integer, input. @@ -114,7 +114,7 @@ ! factorization. ! Note: its allocation is managed by the calling routine psb_ilu_bld, ! hence it cannot be only intent(out). -! info - integer, output. +! info - integer, output. ! Error code. ! blck - type(psb_zspmat_type), input, optional, target. ! The sparse matrix structure containing the remote rows of the @@ -122,9 +122,9 @@ ! to build an Additive Schwarz base preconditioner with overlap ! greater than 0. If the overlap is 0 or the matrix has been reordered ! (see psb_fact_bld), then blck does not contain any row. -! +! subroutine psb_zilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) - + use psb_base_mod use psb_z_ilu_fact_mod, psb_protect_name => psb_zilut_fact @@ -141,7 +141,7 @@ subroutine psb_zilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) integer(psb_ipk_), intent(in), optional :: iscale ! Local Variables integer(psb_ipk_) :: l1, l2, m, err_act, iscale_ - + type(psb_zspmat_type), pointer :: blck_ type(psb_z_csr_sparse_mat) :: ll, uu real(psb_dpk_) :: scale @@ -151,34 +151,34 @@ subroutine psb_zilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) info = psb_success_ call psb_erractionsave(err_act) - if (fill_in < 0) then + if (fill_in < 0) then info=psb_err_input_asize_invalid_i_ call psb_errpush(info,name, & & i_err=(/ione,fill_in,izero,izero,izero/)) goto 9999 end if - ! + ! ! Point to / allocate memory for the incomplete factorization ! - if (present(blck)) then + if (present(blck)) then blck_ => blck else - allocate(blck_,stat=info) - if (info == psb_success_) call blck_%csall(izero,izero,info,ione) + allocate(blck_,stat=info) + if (info == psb_success_) call blck_%allocate(izero,izero,info,ione,type='CSR') if (info /= psb_success_) then info=psb_err_from_subroutine_ - ch_err='csall' + ch_err='allocate' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if endif - if (present(iscale)) then + if (present(iscale)) then iscale_ = iscale else iscale_ = psb_ilu_scale_none_ end if - - select case(iscale_) + + select case(iscale_) case(psb_ilu_scale_none_) scale = sone case(psb_ilu_scale_maxval_) @@ -189,10 +189,10 @@ subroutine psb_zilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) call psb_errpush(info,name,i_err=(/ione*9,iscale_,izero,izero,izero/)) goto 9999 end select - + m = a%get_nrows() + blck_%get_nrows() if ((m /= l%get_nrows()).or.(m /= u%get_nrows()).or.& - & (m > size(d)) ) then + & (m > size(d)) ) then write(0,*) 'Wrong allocation status for L,D,U? ',& & l%get_nrows(),size(d),u%get_nrows() info = -1 @@ -225,15 +225,15 @@ subroutine psb_zilut_fact(fill_in,thres,a,l,u,d,info,blck,iscale) call u%set_triangle() call u%set_unit() call u%set_upper() - + ! ! Nullify pointer / deallocate memory ! - if (present(blck)) then - blck_ => null() + if (present(blck)) then + blck_ => null() else call blck_%free() - deallocate(blck_,stat=info) + deallocate(blck_,stat=info) if(info.ne.0) then info=psb_err_from_subroutine_ ch_err='psb_sp_free' @@ -298,7 +298,7 @@ contains ! according to the CSR storage format. ! lirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the L factor in lval, according to the CSR storage format. + ! of the L factor in lval, according to the CSR storage format. ! uval - complex(psb_dpk_), dimension(:), input/output. ! The U factor in the incomplete factorization. ! The entries of U are stored according to the CSR format. @@ -307,12 +307,12 @@ contains ! according to the CSR storage format. ! uirp - integer, dimension(:), input/output. ! The indices identifying the first nonzero entry of each row - ! of the U factor in uval, according to the CSR storage format. + ! of the U factor in uval, according to the CSR storage format. ! l1 - integer, output ! The number of nonzero entries in lval. ! l2 - integer, output ! The number of nonzero entries in uval. - ! info - integer, output. + ! info - integer, output. ! Error code. ! subroutine psb_zilut_factint(fill_in,thres,a,b,& @@ -320,7 +320,7 @@ contains use psb_base_mod - implicit none + implicit none ! Arguments integer(psb_ipk_), intent(in) :: fill_in @@ -355,7 +355,7 @@ contains m = ma+mb ! - ! Allocate a temporary buffer for the ilut_copyin function + ! Allocate a temporary buffer for the ilut_copyin function ! call trw%allocate(izero,izero,ione) if (info == psb_success_) call psb_ensure_size(m+1,lirp,info) @@ -366,7 +366,7 @@ contains call psb_errpush(info,name,a_err='psb_sp_all') goto 9999 end if - + l1=0 l2=0 lirp(1) = 1 @@ -396,10 +396,10 @@ contains ! in the ilut_copyin function, and updated during the elimination, in ! the ilut_fact routine. The heap is ideal because at each step we need ! the lowest index, but we also need to insert new items, and the heap - ! allows to do both in log time. + ! allows to do both in log time. ! d(i) = czero - if (i<=ma) then + if (i<=ma) then call ilut_copyin(i,ma,a,i,ione,m,nlw,nup,jmaxup,nrmi,weight,& & row,heap,ktrw,trw,info) else @@ -414,7 +414,7 @@ contains & d,uja,uirp,uval,nidx,idxs,info) ! ! Copy the row into lval/d(i)/uval - ! + ! if (info == psb_success_) call ilut_copyout(fill_in,thres,i,m,& & nlw,nup,jmaxup,nrmi,row,nidx,idxs,& & l1,l2,lja,lirp,lval,d,uja,uirp,uval,info) @@ -429,7 +429,7 @@ contains ! ! Adjust diagonal accounting for scale factor ! - if (weight /= sone) then + if (weight /= sone) then d(1:m) = d(1:m)*weight end if @@ -469,7 +469,7 @@ contains ! - storing into a heap the column indices of the nonzero entries of the copied ! row; ! - computing the column index of the first entry with maximum absolute value - ! in the part of the row belonging to the upper triangle; + ! in the part of the row belonging to the upper triangle; ! - computing the 2-norm of the row. ! The output array row is such that it contains a full row of A, i.e. it contains ! also the zero entries of the row. This is useful for the elimination step @@ -482,11 +482,11 @@ contains ! ! This routine is used by psb_zilut_factint in the computation of the ILU(k,t) ! factorization of a local sparse matrix. - ! + ! ! ! Arguments: ! i - integer, input. - ! The local index of the row to be extracted from the + ! The local index of the row to be extracted from the ! sparse matrix structure a. ! m - integer, input. ! The number of rows of the local matrix stored into a. @@ -533,7 +533,7 @@ contains ! staging buffer trw. See below. ! trw - type(psb_zspmat_type), input/output. ! A staging buffer. If the matrix A is not in CSR format, we use - ! the psb_sp_getblk routine and store its output in trw; when we + ! the psb_sp_getblk routine and store its output in trw; when we ! need to call psb_sp_getblk we do it for a block of rows, and then ! we consume them from trw in successive calls to this routine, ! until we empty the buffer. Thus we will make a call to psb_sp_getblk @@ -542,7 +542,7 @@ contains subroutine ilut_copyin(i,m,a,jd,jmin,jmax,nlw,nup,jmaxup,& & nrmi,weight,row,heap,ktrw,trw,info) use psb_base_mod - implicit none + implicit none type(psb_zspmat_type), intent(in) :: a type(psb_z_coo_sparse_mat), intent(inout) :: trw integer(psb_ipk_), intent(in) :: i, m,jmin,jmax,jd @@ -584,29 +584,29 @@ contains dmaxup = szero nrmi = szero - select type (aa=> a%a) - type is (psb_z_csr_sparse_mat) + select type (aa=> a%a) + type is (psb_z_csr_sparse_mat) ! ! Take a fast shortcut if the matrix is stored in CSR format ! do j = aa%irp(i), aa%irp(i+1) - 1 k = aa%ja(j) - if ((jmin<=k).and.(k<=jmax)) then - row(k) = aa%val(j)*weight + if ((jmin<=k).and.(k<=jmax)) then + row(k) = aa%val(j)*weight call heap%insert(k,info) if (info /= psb_success_) exit - if (kjd) then + if (kjd) then nup = nup + 1 - if (abs(row(k))>dmaxup) then + if (abs(row(k))>dmaxup) then jmaxup = k dmaxup = abs(row(k)) end if end if end if end do - if (info /= psb_success_) then + if (info /= psb_success_) then info=psb_err_from_subroutine_ call psb_errpush(info,name,a_err='psb_insert_heap') goto 9999 @@ -617,16 +617,16 @@ contains class default - + ! - ! Otherwise use psb_sp_getblk, slower but able (in principle) of + ! Otherwise use psb_sp_getblk, slower but able (in principle) of ! handling any format. In this case, a block of rows is extracted ! instead of a single row, for performance reasons, and these ! rows are copied one by one into the array row, through successive ! calls to ilut_copyin. ! - if ((mod(i,nrb) == 1).or.(nrb == 1)) then + if ((mod(i,nrb) == 1).or.(nrb == 1)) then irb = min(m-i+1,nrb) call aa%csget(i,i+irb-1,trw,info) if (info /= psb_success_) then @@ -639,18 +639,18 @@ contains kin = ktrw nz = trw%get_nzeros() - do + do if (ktrw > nz) exit if (trw%ia(ktrw) > i) exit k = trw%ja(ktrw) - if ((jmin<=k).and.(k<=jmax)) then + if ((jmin<=k).and.(k<=jmax)) then row(k) = trw%val(ktrw)*weight call heap%insert(k,info) if (info /= psb_success_) exit - if (kjd) then + if (kjd) then nup = nup + 1 - if (abs(row(k))>dmaxup) then + if (abs(row(k))>dmaxup) then jmaxup = k dmaxup = abs(row(k)) end if @@ -734,10 +734,10 @@ contains use psb_base_mod - implicit none + implicit none ! Arguments - type(psb_i_heap), intent(inout) :: heap + type(psb_i_heap), intent(inout) :: heap integer(psb_ipk_), intent(in) :: i integer(psb_ipk_), intent(inout) :: nidx,info real(psb_dpk_), intent(in) :: thres,nrmi @@ -753,23 +753,23 @@ contains call psb_ensure_size(200*ione,idxs,info) if (info /= psb_success_) return nidx = 0 - lastk = -1 + lastk = -1 ! ! Do while there are indices to be processed ! do - call heap%get_first(k,iret) + call heap%get_first(k,iret) if (iret < 0) exit - ! + ! ! An index may have been put on the heap more than once. ! if (k == lastk) cycle - lastk = k + lastk = k lowert: if (k nidx) exit if (idxs(idxp) >= i) exit @@ -987,8 +987,8 @@ contains ! if (abs(witem) < thres*nrmi) cycle - nz = nz + 1 - xw(nz) = witem + nz = nz + 1 + xw(nz) = witem xwid(nz) = widx call heap%insert(witem,widx,info) if (info /= psb_success_) then @@ -1000,9 +1000,9 @@ contains ! ! Now we have to take out the first nlw+fill_in entries - ! + ! if (nz <= nlw+fill_in) then - ! + ! ! Just copy everything from xw, and it is already ordered ! else @@ -1014,7 +1014,7 @@ contains call psb_errpush(info,name,a_err='psb_heap_get_first') goto 9999 end if - + xw(k) = witem xwid(k) = widx end do @@ -1029,15 +1029,15 @@ contains ! Copy out the lower part of the row ! do k=1,nz - l1 = l1 + 1 + l1 = l1 + 1 if (size(lval) < l1) then - ! + ! ! Figure out a good reallocation size! - ! + ! isz = (max((l1/i)*m,int(1.2*l1),l1+100)) - call psb_realloc(isz,lval,info) - if (info == psb_success_) call psb_realloc(isz,lja,info) - if (info /= psb_success_) then + call psb_realloc(isz,lval,info) + if (info == psb_success_) call psb_realloc(isz,lja,info) + if (info /= psb_success_) then info=psb_err_from_subroutine_ call psb_errpush(info,name,a_err='Allocate') goto 9999 @@ -1046,25 +1046,25 @@ contains lja(l1) = xwid(k) lval(l1) = xw(indx(k)) end do - + ! ! Make sure idxp points to the diagonal entry ! - if (idxp <= size(idxs)) then - if (idxs(idxp) < i) then - do + if (idxp <= size(idxs)) then + if (idxs(idxp) < i) then + do idxp = idxp + 1 if (idxp > nidx) exit if (idxs(idxp) >= i) exit end do end if end if - if (idxp > size(idxs)) then + if (idxp > size(idxs)) then !!$ write(0,*) 'Warning: missing diagonal element in the row ' else - if (idxs(idxp) > i) then + if (idxs(idxp) > i) then !!$ write(0,*) 'Warning: missing diagonal element in the row ' - else if (idxs(idxp) /= i) then + else if (idxs(idxp) /= i) then !!$ write(0,*) 'Warning: impossible error: diagonal has vanished' else ! @@ -1076,7 +1076,7 @@ contains if (abs(d(i)) < d_epstol) then ! ! Too small pivot: unstable factorization - ! + ! info = psb_err_pivot_too_small_ int_err(1) = i write(ch_err,'(g20.10)') d(i) @@ -1092,7 +1092,7 @@ contains end if ! - ! Now the upper part + ! Now the upper part ! call heap%init(info,dir=psb_asort_down_) @@ -1104,15 +1104,15 @@ contains nz = 0 do - + idxp = idxp + 1 if (idxp > nidx) exit widx = idxs(idxp) - if (widx <= i) then + if (widx <= i) then !!$ write(0,*) 'Warning: lower triangle in upper copy',widx,i,idxp,idxs(idxp) cycle end if - if (widx > m) then + if (widx > m) then !!$ write(0,*) 'Warning: impossible value',widx,i,idxp,idxs(idxp) cycle end if @@ -1120,12 +1120,12 @@ contains ! ! Dropping rule based on the 2-norm. But keep the jmaxup-th entry anyway. ! - if ((widx /= jmaxup) .and. (abs(witem) < thres*nrmi)) then - cycle + if ((widx /= jmaxup) .and. (abs(witem) < thres*nrmi)) then + cycle end if nz = nz + 1 - xw(nz) = witem + xw(nz) = witem xwid(nz) = widx call heap%insert(witem,widx,info) if (info /= psb_success_) then @@ -1133,7 +1133,7 @@ contains call psb_errpush(info,name,a_err='psb_insert_heap') goto 9999 end if - + end do ! @@ -1141,7 +1141,7 @@ contains ! we include entry jmaxup. ! if (nz <= nup+fill_in) then - ! + ! ! Just copy everything from xw ! fndmaxup=.true. @@ -1155,13 +1155,13 @@ contains if (widx == jmaxup) fndmaxup=.true. end do end if - if ((i