From c8fe04999431028859c29d676483603da79e9b8c Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 22 Nov 2010 13:58:47 +0000 Subject: [PATCH] psblas3: prec/psb_c_base_prec_mod.f03 prec/psb_c_bjacprec.f03 prec/psb_c_diagprec.f03 prec/psb_c_nullprec.f03 prec/psb_cprc_aply.f90 prec/psb_cprecbld.f90 prec/psb_d_base_prec_mod.f03 prec/psb_d_bjacprec.f03 prec/psb_d_diagprec.f03 prec/psb_d_nullprec.f03 prec/psb_dprc_aply.f90 prec/psb_dprecbld.f90 prec/psb_prec_mod.f90 prec/psb_s_base_prec_mod.f03 prec/psb_s_bjacprec.f03 prec/psb_s_diagprec.f03 prec/psb_s_nullprec.f03 prec/psb_sprc_aply.f90 prec/psb_sprecbld.f90 prec/psb_z_base_prec_mod.f03 prec/psb_z_bjacprec.f03 prec/psb_z_diagprec.f03 prec/psb_z_nullprec.f03 prec/psb_zprc_aply.f90 prec/psb_zprecbld.f90 test/newfmt test/newfmt/Makefile test/newfmt/ppde.f90 test/newfmt/runs test/newfmt/runs/Makefile test/newfmt/runs/ppde.inp test/newfmt/spde.f90 1. New test/newfmt directory in which to test for new storage versions. 2. New MOLD/AFMT arguments to PRECBLD for storing the matrices. In preconditioners such as DIAG and NULL they are ignored. --- prec/psb_c_base_prec_mod.f03 | 4 +- prec/psb_c_bjacprec.f03 | 23 +- prec/psb_c_diagprec.f03 | 10 +- prec/psb_c_nullprec.f03 | 10 +- prec/psb_cprc_aply.f90 | 195 ---------- prec/psb_cprecbld.f90 | 16 +- prec/psb_d_base_prec_mod.f03 | 10 +- prec/psb_d_bjacprec.f03 | 23 +- prec/psb_d_diagprec.f03 | 10 +- prec/psb_d_nullprec.f03 | 10 +- prec/psb_dprc_aply.f90 | 194 ---------- prec/psb_dprecbld.f90 | 17 +- prec/psb_prec_mod.f90 | 28 +- prec/psb_s_base_prec_mod.f03 | 4 +- prec/psb_s_bjacprec.f03 | 24 +- prec/psb_s_diagprec.f03 | 10 +- prec/psb_s_nullprec.f03 | 10 +- prec/psb_sprc_aply.f90 | 194 ---------- prec/psb_sprecbld.f90 | 17 +- prec/psb_z_base_prec_mod.f03 | 10 +- prec/psb_z_bjacprec.f03 | 28 +- prec/psb_z_diagprec.f03 | 10 +- prec/psb_z_nullprec.f03 | 10 +- prec/psb_zprc_aply.f90 | 195 ---------- prec/psb_zprecbld.f90 | 17 +- test/newfmt/Makefile | 38 ++ test/newfmt/ppde.f90 | 671 +++++++++++++++++++++++++++++++++++ test/newfmt/runs/Makefile | 12 + test/newfmt/runs/ppde.inp | 11 + test/newfmt/spde.f90 | 664 ++++++++++++++++++++++++++++++++++ 30 files changed, 1557 insertions(+), 918 deletions(-) delete mode 100644 prec/psb_cprc_aply.f90 delete mode 100644 prec/psb_dprc_aply.f90 delete mode 100644 prec/psb_sprc_aply.f90 delete mode 100644 prec/psb_zprc_aply.f90 create mode 100644 test/newfmt/Makefile create mode 100644 test/newfmt/ppde.f90 create mode 100644 test/newfmt/runs/Makefile create mode 100644 test/newfmt/runs/ppde.inp create mode 100644 test/newfmt/spde.f90 diff --git a/prec/psb_c_base_prec_mod.f03 b/prec/psb_c_base_prec_mod.f03 index 68e0525e..d251ef02 100644 --- a/prec/psb_c_base_prec_mod.f03 +++ b/prec/psb_c_base_prec_mod.f03 @@ -132,7 +132,7 @@ contains return end subroutine psb_c_base_precinit - subroutine psb_c_base_precbld(a,desc_a,prec,info,upd) + subroutine psb_c_base_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod Implicit None @@ -142,6 +142,8 @@ contains class(psb_c_base_prec_type),intent(inout) :: prec integer, intent(out) :: info character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_c_base_sparse_mat), intent(in), optional :: mold Integer :: err_act, nrow character(len=20) :: name='c_base_precbld' diff --git a/prec/psb_c_bjacprec.f03 b/prec/psb_c_bjacprec.f03 index 02a0c73b..ffcc5133 100644 --- a/prec/psb_c_bjacprec.f03 +++ b/prec/psb_c_bjacprec.f03 @@ -215,17 +215,19 @@ contains end subroutine psb_c_bjac_precinit - subroutine psb_c_bjac_precbld(a,desc_a,prec,info,upd) + subroutine psb_c_bjac_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod use psb_prec_mod Implicit None type(psb_cspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_c_bjac_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_c_base_sparse_mat), intent(in), optional :: mold ! .. Local Scalars .. integer :: i, m @@ -325,12 +327,6 @@ contains goto 9999 end if -!!$ call prec%av(psb_l_pr_)%print(30+me) -!!$ call prec%av(psb_u_pr_)%print(40+me) -!!$ do i=1,n_row -!!$ write(50+me,*) i,prec%d(i) -!!$ end do - case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' @@ -344,6 +340,13 @@ contains goto 9999 end select + if (present(mold)) then + call prec%av(psb_l_pr_)%cscnv(info,mold=mold) + call prec%av(psb_u_pr_)%cscnv(info,mold=mold) + else if (present(afmt)) then + call prec%av(psb_l_pr_)%cscnv(info,type=afmt) + call prec%av(psb_u_pr_)%cscnv(info,type=afmt) + end if call psb_erractionrestore(err_act) return diff --git a/prec/psb_c_diagprec.f03 b/prec/psb_c_diagprec.f03 index d4b91992..5038b374 100644 --- a/prec/psb_c_diagprec.f03 +++ b/prec/psb_c_diagprec.f03 @@ -153,16 +153,18 @@ contains end subroutine psb_c_diag_precinit - subroutine psb_c_diag_precbld(a,desc_a,prec,info,upd) + subroutine psb_c_diag_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod Implicit None type(psb_cspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_c_diag_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_c_base_sparse_mat), intent(in), optional :: mold Integer :: err_act, nrow,i character(len=20) :: name='c_diag_precbld' diff --git a/prec/psb_c_nullprec.f03 b/prec/psb_c_nullprec.f03 index f961cae5..b9f42bdd 100644 --- a/prec/psb_c_nullprec.f03 +++ b/prec/psb_c_nullprec.f03 @@ -103,16 +103,18 @@ contains return end subroutine psb_c_null_precinit - subroutine psb_c_null_precbld(a,desc_a,prec,info,upd) + subroutine psb_c_null_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod Implicit None type(psb_cspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_c_null_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_c_base_sparse_mat), intent(in), optional :: mold Integer :: err_act, nrow character(len=20) :: name='c_null_precbld' diff --git a/prec/psb_cprc_aply.f90 b/prec/psb_cprc_aply.f90 deleted file mode 100644 index 0462aa92..00000000 --- a/prec/psb_cprc_aply.f90 +++ /dev/null @@ -1,195 +0,0 @@ -!!$ -!!$ Parallel Sparse BLAS version 3.0 -!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 -!!$ Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari CNRS-IRIT, Toulouse -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the PSBLAS group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -subroutine psb_cprc_aply(prec,x,y,desc_data,info,trans, work) - - use psb_sparse_mod - use psb_prec_mod, psb_protect_name => psb_cprc_aply - - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_cprec_type), intent(in) :: prec - complex(psb_spk_),intent(in) :: x(:) - complex(psb_spk_),intent(inout) :: y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - complex(psb_spk_), intent(inout), optional, target :: work(:) - - ! Local variables - character :: trans_ - complex(psb_spk_), pointer :: work_(:) - integer :: ictxt,np,me,err_act - character(len=20) :: name - - name='psb_prc_aply' - info = psb_success_ - call psb_erractionsave(err_act) - - ictxt=desc_data%matrix_data(psb_ctxt_) - call psb_info(ictxt, me, np) - - if (present(trans)) then - trans_=trans - else - trans_='N' - end if - - if (present(work)) then - work_ => work - else - allocate(work_(4*desc_data%matrix_data(psb_n_col_)),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - - end if - - if (.not.allocated(prec%prec)) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner") - goto 9999 - end if - call prec%prec%apply(cone,x,czero,y,desc_data,info,trans_,work=work_) - if (present(work)) then - else - deallocate(work_) - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - -end subroutine psb_cprc_aply - - -!!$ -!!$ -!!$ MD2P4 -!!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS -!!$ for -!!$ Parallel Sparse BLAS version 3.0 -!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 -!!$ Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari CNRS-IRIT, Toulouse -!!$ Daniela di Serafino Second University of Naples -!!$ Pasqua D'Ambra ICAR-CNR -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the MD2P4 group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MD2P4 GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -subroutine psb_cprc_aply1(prec,x,desc_data,info,trans) - - use psb_sparse_mod - use psb_prec_mod, psb_protect_name => psb_cprc_aply1 - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_cprec_type), intent(in) :: prec - complex(psb_spk_),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - - ! Local variables - character :: trans_ - integer :: ictxt,np,me, err_act - complex(psb_spk_), pointer :: WW(:), w1(:) - character(len=20) :: name - name='psb_prc_aply1' - info = psb_success_ - call psb_erractionsave(err_act) - - - ictxt=desc_data%matrix_data(psb_ctxt_) - call psb_info(ictxt, me, np) - if (present(trans)) then - trans_=psb_toupper(trans) - else - trans_='N' - end if - - if (.not.allocated(prec%prec)) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner") - goto 9999 - end if - allocate(ww(size(x)),w1(size(x)),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - call prec%prec%apply(cone,x,czero,ww,desc_data,info,trans_,work=w1) - if(info /= psb_success_) goto 9999 - x(:) = ww(:) - deallocate(ww,W1) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_errpush(info,name) - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return -end subroutine psb_cprc_aply1 diff --git a/prec/psb_cprecbld.f90 b/prec/psb_cprecbld.f90 index d72578cf..938f3387 100644 --- a/prec/psb_cprecbld.f90 +++ b/prec/psb_cprecbld.f90 @@ -29,7 +29,7 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_cprecbld(a,desc_a,p,info,upd) +subroutine psb_cprecbld(a,desc_a,p,info,upd,mold,afmt) use psb_sparse_mod use psb_prec_mod, psb_protect_name => psb_cprecbld @@ -40,6 +40,8 @@ subroutine psb_cprecbld(a,desc_a,p,info,upd) type(psb_cprec_type),intent(inout) :: p integer, intent(out) :: info character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_c_base_sparse_mat), intent(in), optional :: mold ! Local scalars @@ -63,16 +65,6 @@ subroutine psb_cprecbld(a,desc_a,p,info,upd) call psb_info(ictxt, me, np) - if (present(upd)) then - upd_ = psb_toupper(upd) - else - upd_='F' - endif - if ((upd_ == 'F').or.(upd_ == 'T')) then - ! ok - else - upd_='F' - endif n_row = psb_cd_get_local_rows(desc_a) n_col = psb_cd_get_local_cols(desc_a) mglob = psb_cd_get_global_rows(desc_a) @@ -88,7 +80,7 @@ subroutine psb_cprecbld(a,desc_a,p,info,upd) goto 9999 end if - call p%prec%precbld(a,desc_a,info,upd) + call p%prec%precbld(a,desc_a,info,upd=upd,afmt=afmt,mold=mold) if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) diff --git a/prec/psb_d_base_prec_mod.f03 b/prec/psb_d_base_prec_mod.f03 index eb48ea79..0ff49b56 100644 --- a/prec/psb_d_base_prec_mod.f03 +++ b/prec/psb_d_base_prec_mod.f03 @@ -138,16 +138,18 @@ contains return end subroutine psb_d_base_precinit - subroutine psb_d_base_precbld(a,desc_a,prec,info,upd) + subroutine psb_d_base_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod Implicit None type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_d_base_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_d_base_sparse_mat), intent(in), optional :: mold Integer :: err_act, nrow character(len=20) :: name='d_base_precbld' diff --git a/prec/psb_d_bjacprec.f03 b/prec/psb_d_bjacprec.f03 index 9b1e8d54..43a66d86 100644 --- a/prec/psb_d_bjacprec.f03 +++ b/prec/psb_d_bjacprec.f03 @@ -210,17 +210,19 @@ contains end subroutine psb_d_bjac_precinit - subroutine psb_d_bjac_precbld(a,desc_a,prec,info,upd) + subroutine psb_d_bjac_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod use psb_prec_mod Implicit None type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_d_bjac_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_d_base_sparse_mat), intent(in), optional :: mold ! .. Local Scalars .. integer :: i, m @@ -322,12 +324,6 @@ contains goto 9999 end if -!!$ call prec%av(psb_l_pr_)%print(30+me) -!!$ call prec%av(psb_u_pr_)%print(40+me) -!!$ do i=1,n_row -!!$ write(50+me,*) i,prec%d(i) -!!$ end do - case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' @@ -341,6 +337,13 @@ contains goto 9999 end select + if (present(mold)) then + call prec%av(psb_l_pr_)%cscnv(info,mold=mold) + call prec%av(psb_u_pr_)%cscnv(info,mold=mold) + else if (present(afmt)) then + call prec%av(psb_l_pr_)%cscnv(info,type=afmt) + call prec%av(psb_u_pr_)%cscnv(info,type=afmt) + end if call psb_erractionrestore(err_act) return diff --git a/prec/psb_d_diagprec.f03 b/prec/psb_d_diagprec.f03 index 03b7c9be..2b371cf2 100644 --- a/prec/psb_d_diagprec.f03 +++ b/prec/psb_d_diagprec.f03 @@ -130,16 +130,18 @@ contains end subroutine psb_d_diag_precinit - subroutine psb_d_diag_precbld(a,desc_a,prec,info,upd) + subroutine psb_d_diag_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod Implicit None type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_d_diag_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_d_base_sparse_mat), intent(in), optional :: mold Integer :: err_act, nrow,i character(len=20) :: name='d_diag_precbld' diff --git a/prec/psb_d_nullprec.f03 b/prec/psb_d_nullprec.f03 index 2a7308fc..8d05df38 100644 --- a/prec/psb_d_nullprec.f03 +++ b/prec/psb_d_nullprec.f03 @@ -103,16 +103,18 @@ contains return end subroutine psb_d_null_precinit - subroutine psb_d_null_precbld(a,desc_a,prec,info,upd) + subroutine psb_d_null_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod Implicit None type(psb_dspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_d_null_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_d_base_sparse_mat), intent(in), optional :: mold Integer :: err_act, nrow character(len=20) :: name='d_null_precbld' diff --git a/prec/psb_dprc_aply.f90 b/prec/psb_dprc_aply.f90 deleted file mode 100644 index edca9b9e..00000000 --- a/prec/psb_dprc_aply.f90 +++ /dev/null @@ -1,194 +0,0 @@ -!!$ -!!$ Parallel Sparse BLAS version 3.0 -!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 -!!$ Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari CNRS-IRIT, Toulouse -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the PSBLAS group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -subroutine psb_dprc_aply(prec,x,y,desc_data,info,trans, work) - - use psb_sparse_mod - use psb_prec_mod, psb_protect_name => psb_dprc_aply - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_dprec_type), intent(in) :: prec - real(psb_dpk_),intent(in) :: x(:) - real(psb_dpk_),intent(inout) :: y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - real(psb_dpk_), intent(inout), optional, target :: work(:) - - ! Local variables - character :: trans_ - real(psb_dpk_), pointer :: work_(:) - integer :: ictxt,np,me,err_act - character(len=20) :: name - - name='psb_prc_aply' - info = psb_success_ - call psb_erractionsave(err_act) - - ictxt=desc_data%matrix_data(psb_ctxt_) - call psb_info(ictxt, me, np) - - if (present(trans)) then - trans_=trans - else - trans_='N' - end if - - if (present(work)) then - work_ => work - else - allocate(work_(4*desc_data%matrix_data(psb_n_col_)),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - - end if - - if (.not.allocated(prec%prec)) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner") - goto 9999 - end if - call prec%prec%apply(done,x,dzero,y,desc_data,info,trans_,work=work_) - if (present(work)) then - else - deallocate(work_) - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - -end subroutine psb_dprc_aply - - -!!$ -!!$ -!!$ MD2P4 -!!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS -!!$ for -!!$ Parallel Sparse BLAS version 3.0 -!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 -!!$ Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari CNRS-IRIT, Toulouse -!!$ Daniela di Serafino Second University of Naples -!!$ Pasqua D'Ambra ICAR-CNR -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the MD2P4 group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MD2P4 GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -subroutine psb_dprc_aply1(prec,x,desc_data,info,trans) - - use psb_sparse_mod - use psb_prec_mod, psb_protect_name => psb_dprc_aply1 - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_dprec_type), intent(in) :: prec - real(psb_dpk_),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - - ! Local variables - character :: trans_ - integer :: ictxt,np,me, err_act - real(psb_dpk_), pointer :: WW(:), w1(:) - character(len=20) :: name - name='psb_prc_aply1' - info = psb_success_ - call psb_erractionsave(err_act) - - - ictxt=desc_data%matrix_data(psb_ctxt_) - call psb_info(ictxt, me, np) - if (present(trans)) then - trans_=psb_toupper(trans) - else - trans_='N' - end if - - if (.not.allocated(prec%prec)) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner") - goto 9999 - end if - allocate(ww(size(x)),w1(size(x)),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - call prec%prec%apply(done,x,dzero,ww,desc_data,info,trans_,work=w1) - if(info /= psb_success_) goto 9999 - x(:) = ww(:) - deallocate(ww,W1) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_errpush(info,name) - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return -end subroutine psb_dprc_aply1 diff --git a/prec/psb_dprecbld.f90 b/prec/psb_dprecbld.f90 index 0555615b..d1f263d5 100644 --- a/prec/psb_dprecbld.f90 +++ b/prec/psb_dprecbld.f90 @@ -29,7 +29,7 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_dprecbld(a,desc_a,p,info,upd) +subroutine psb_dprecbld(a,desc_a,p,info,upd,mold,afmt) use psb_sparse_mod use psb_prec_mod, psb_protect_name => psb_dprecbld @@ -40,6 +40,8 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) type(psb_dprec_type),intent(inout) :: p integer, intent(out) :: info character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_d_base_sparse_mat), intent(in), optional :: mold ! Local scalars Integer :: err, n_row, n_col,ictxt,& @@ -62,16 +64,6 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) call psb_info(ictxt, me, np) - if (present(upd)) then - upd_ = psb_toupper(upd) - else - upd_='F' - endif - if ((upd_ == 'F').or.(upd_ == 'T')) then - ! ok - else - upd_='F' - endif n_row = psb_cd_get_local_rows(desc_a) n_col = psb_cd_get_local_cols(desc_a) mglob = psb_cd_get_global_rows(desc_a) @@ -87,7 +79,8 @@ subroutine psb_dprecbld(a,desc_a,p,info,upd) goto 9999 end if - call p%prec%precbld(a,desc_a,info,upd) + call p%prec%precbld(a,desc_a,info,upd=upd,afmt=afmt,mold=mold) + if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) diff --git a/prec/psb_prec_mod.f90 b/prec/psb_prec_mod.f90 index a2f0dd71..1ce8e469 100644 --- a/prec/psb_prec_mod.f90 +++ b/prec/psb_prec_mod.f90 @@ -34,8 +34,9 @@ module psb_prec_mod use psb_prec_type interface psb_precbld - subroutine psb_sprecbld(a,desc_a,prec,info,upd) - use psb_sparse_mod, only : psb_desc_type, psb_sspmat_type, psb_spk_ + subroutine psb_sprecbld(a,desc_a,prec,info,upd,mold,afmt) + use psb_sparse_mod, only : psb_desc_type, psb_sspmat_type,& + & psb_s_base_sparse_mat, psb_spk_ use psb_prec_type, only : psb_sprec_type implicit none type(psb_sspmat_type), intent(in), target :: a @@ -43,9 +44,12 @@ module psb_prec_mod type(psb_sprec_type), intent(inout) :: prec integer, intent(out) :: info character, intent(in),optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_s_base_sparse_mat), intent(in), optional :: mold end subroutine psb_sprecbld - subroutine psb_dprecbld(a,desc_a,prec,info,upd) - use psb_sparse_mod, only : psb_desc_type, psb_dspmat_type, psb_dpk_ + subroutine psb_dprecbld(a,desc_a,prec,info,upd,mold,afmt) + use psb_sparse_mod, only : psb_desc_type, psb_dspmat_type,& + & psb_d_base_sparse_mat, psb_dpk_ use psb_prec_type, only : psb_dprec_type implicit none type(psb_dspmat_type), intent(in), target :: a @@ -53,9 +57,12 @@ module psb_prec_mod type(psb_dprec_type), intent(inout) :: prec integer, intent(out) :: info character, intent(in),optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_d_base_sparse_mat), intent(in), optional :: mold end subroutine psb_dprecbld - subroutine psb_cprecbld(a,desc_a,prec,info,upd) - use psb_sparse_mod, only : psb_desc_type, psb_cspmat_type, psb_spk_ + subroutine psb_cprecbld(a,desc_a,prec,info,upd,mold,afmt) + use psb_sparse_mod, only : psb_desc_type, psb_cspmat_type,& + & psb_c_base_sparse_mat, psb_spk_ use psb_prec_type, only : psb_cprec_type implicit none type(psb_cspmat_type), intent(in), target :: a @@ -63,9 +70,12 @@ module psb_prec_mod type(psb_cprec_type), intent(inout) :: prec integer, intent(out) :: info character, intent(in),optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_c_base_sparse_mat), intent(in), optional :: mold end subroutine psb_cprecbld - subroutine psb_zprecbld(a,desc_a,prec,info,upd) - use psb_sparse_mod, only : psb_desc_type, psb_zspmat_type, psb_dpk_ + subroutine psb_zprecbld(a,desc_a,prec,info,upd,mold,afmt) + use psb_sparse_mod, only : psb_desc_type, psb_zspmat_type,& + & psb_z_base_sparse_mat, psb_dpk_ use psb_prec_type, only : psb_zprec_type implicit none type(psb_zspmat_type), intent(in), target :: a @@ -73,6 +83,8 @@ module psb_prec_mod type(psb_zprec_type), intent(inout) :: prec integer, intent(out) :: info character, intent(in),optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_z_base_sparse_mat), intent(in), optional :: mold end subroutine psb_zprecbld end interface diff --git a/prec/psb_s_base_prec_mod.f03 b/prec/psb_s_base_prec_mod.f03 index bf377d5d..b1445bfe 100644 --- a/prec/psb_s_base_prec_mod.f03 +++ b/prec/psb_s_base_prec_mod.f03 @@ -132,7 +132,7 @@ contains return end subroutine psb_s_base_precinit - subroutine psb_s_base_precbld(a,desc_a,prec,info,upd) + subroutine psb_s_base_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod Implicit None @@ -142,6 +142,8 @@ contains class(psb_s_base_prec_type),intent(inout) :: prec integer, intent(out) :: info character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_s_base_sparse_mat), intent(in), optional :: mold Integer :: err_act, nrow character(len=20) :: name='s_base_precbld' diff --git a/prec/psb_s_bjacprec.f03 b/prec/psb_s_bjacprec.f03 index 52086c2d..72f1e69e 100644 --- a/prec/psb_s_bjacprec.f03 +++ b/prec/psb_s_bjacprec.f03 @@ -209,17 +209,19 @@ contains end subroutine psb_s_bjac_precinit - subroutine psb_s_bjac_precbld(a,desc_a,prec,info,upd) + subroutine psb_s_bjac_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod use psb_prec_mod Implicit None type(psb_sspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_s_bjac_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_s_base_sparse_mat), intent(in), optional :: mold ! .. Local Scalars .. integer :: i, m @@ -319,12 +321,6 @@ contains goto 9999 end if -!!$ call prec%av(psb_l_pr_)%print(30+me) -!!$ call prec%av(psb_u_pr_)%print(40+me) -!!$ do i=1,n_row -!!$ write(50+me,*) i,prec%d(i) -!!$ end do - case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' @@ -338,6 +334,14 @@ contains goto 9999 end select + if (present(mold)) then + call prec%av(psb_l_pr_)%cscnv(info,mold=mold) + call prec%av(psb_u_pr_)%cscnv(info,mold=mold) + else if (present(afmt)) then + call prec%av(psb_l_pr_)%cscnv(info,type=afmt) + call prec%av(psb_u_pr_)%cscnv(info,type=afmt) + end if + call psb_erractionrestore(err_act) return diff --git a/prec/psb_s_diagprec.f03 b/prec/psb_s_diagprec.f03 index aff2b24f..ff886bf9 100644 --- a/prec/psb_s_diagprec.f03 +++ b/prec/psb_s_diagprec.f03 @@ -130,16 +130,18 @@ contains end subroutine psb_s_diag_precinit - subroutine psb_s_diag_precbld(a,desc_a,prec,info,upd) + subroutine psb_s_diag_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod Implicit None type(psb_sspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_s_diag_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_s_base_sparse_mat), intent(in), optional :: mold Integer :: err_act, nrow,i character(len=20) :: name='s_diag_precbld' diff --git a/prec/psb_s_nullprec.f03 b/prec/psb_s_nullprec.f03 index 893e25a2..aec3c9ba 100644 --- a/prec/psb_s_nullprec.f03 +++ b/prec/psb_s_nullprec.f03 @@ -102,16 +102,18 @@ contains return end subroutine psb_s_null_precinit - subroutine psb_s_null_precbld(a,desc_a,prec,info,upd) + subroutine psb_s_null_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod Implicit None type(psb_sspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_s_null_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_s_base_sparse_mat), intent(in), optional :: mold Integer :: err_act, nrow character(len=20) :: name='s_null_precbld' diff --git a/prec/psb_sprc_aply.f90 b/prec/psb_sprc_aply.f90 deleted file mode 100644 index 241d83e7..00000000 --- a/prec/psb_sprc_aply.f90 +++ /dev/null @@ -1,194 +0,0 @@ -!!$ -!!$ Parallel Sparse BLAS version 3.0 -!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 -!!$ Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari CNRS-IRIT, Toulouse -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the PSBLAS group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -subroutine psb_sprc_aply(prec,x,y,desc_data,info,trans, work) - - use psb_sparse_mod - use psb_prec_mod, psb_protect_name => psb_sprc_aply - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_sprec_type), intent(in) :: prec - real(psb_spk_),intent(in) :: x(:) - real(psb_spk_),intent(inout) :: y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - real(psb_spk_), intent(inout), optional, target :: work(:) - - ! Local variables - character :: trans_ - real(psb_spk_), pointer :: work_(:) - integer :: ictxt,np,me,err_act - character(len=20) :: name - - name='psb_prc_aply' - info = psb_success_ - call psb_erractionsave(err_act) - - ictxt=desc_data%matrix_data(psb_ctxt_) - call psb_info(ictxt, me, np) - - if (present(trans)) then - trans_=trans - else - trans_='N' - end if - - if (present(work)) then - work_ => work - else - allocate(work_(4*desc_data%matrix_data(psb_n_col_)),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - - end if - - if (.not.allocated(prec%prec)) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner") - goto 9999 - end if - call prec%prec%apply(sone,x,szero,y,desc_data,info,trans_,work=work_) - if (present(work)) then - else - deallocate(work_) - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - -end subroutine psb_sprc_aply - - -!!$ -!!$ -!!$ MD2P4 -!!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS -!!$ for -!!$ Parallel Sparse BLAS version 3.0 -!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 -!!$ Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari CNRS-IRIT, Toulouse -!!$ Daniela di Serafino Second University of Naples -!!$ Pasqua D'Ambra ICAR-CNR -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the MD2P4 group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MD2P4 GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -subroutine psb_sprc_aply1(prec,x,desc_data,info,trans) - - use psb_sparse_mod - use psb_prec_mod, psb_protect_name => psb_sprc_aply1 - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_sprec_type), intent(in) :: prec - real(psb_spk_),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - - ! Local variables - character :: trans_ - integer :: ictxt,np,me, err_act - real(psb_spk_), pointer :: WW(:), w1(:) - character(len=20) :: name - name='psb_prc_aply1' - info = psb_success_ - call psb_erractionsave(err_act) - - - ictxt=desc_data%matrix_data(psb_ctxt_) - call psb_info(ictxt, me, np) - if (present(trans)) then - trans_=psb_toupper(trans) - else - trans_='N' - end if - - if (.not.allocated(prec%prec)) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner") - goto 9999 - end if - allocate(ww(size(x)),w1(size(x)),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - call prec%prec%apply(sone,x,szero,ww,desc_data,info,trans_,work=w1) - if(info /= psb_success_) goto 9999 - x(:) = ww(:) - deallocate(ww,W1) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_errpush(info,name) - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return -end subroutine psb_sprc_aply1 diff --git a/prec/psb_sprecbld.f90 b/prec/psb_sprecbld.f90 index 1ec1ccde..99f0f4cf 100644 --- a/prec/psb_sprecbld.f90 +++ b/prec/psb_sprecbld.f90 @@ -29,7 +29,7 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_sprecbld(a,desc_a,p,info,upd) +subroutine psb_sprecbld(a,desc_a,p,info,upd,mold,afmt) use psb_sparse_mod use psb_prec_mod, psb_protect_name => psb_sprecbld @@ -40,6 +40,8 @@ subroutine psb_sprecbld(a,desc_a,p,info,upd) type(psb_sprec_type),intent(inout) :: p integer, intent(out) :: info character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_s_base_sparse_mat), intent(in), optional :: mold ! Local scalars Integer :: err, n_row, n_col,ictxt,& @@ -62,16 +64,6 @@ subroutine psb_sprecbld(a,desc_a,p,info,upd) call psb_info(ictxt, me, np) - if (present(upd)) then - upd_ = psb_toupper(upd) - else - upd_='F' - endif - if ((upd_ == 'F').or.(upd_ == 'T')) then - ! ok - else - upd_='F' - endif n_row = psb_cd_get_local_rows(desc_a) n_col = psb_cd_get_local_cols(desc_a) mglob = psb_cd_get_global_rows(desc_a) @@ -87,7 +79,8 @@ subroutine psb_sprecbld(a,desc_a,p,info,upd) goto 9999 end if - call p%prec%precbld(a,desc_a,info,upd) + call p%prec%precbld(a,desc_a,info,upd=upd,afmt=afmt,mold=mold) + if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) diff --git a/prec/psb_z_base_prec_mod.f03 b/prec/psb_z_base_prec_mod.f03 index f259b12b..3c99415c 100644 --- a/prec/psb_z_base_prec_mod.f03 +++ b/prec/psb_z_base_prec_mod.f03 @@ -133,16 +133,18 @@ contains return end subroutine psb_z_base_precinit - subroutine psb_z_base_precbld(a,desc_a,prec,info,upd) + subroutine psb_z_base_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod Implicit None type(psb_zspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_z_base_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_z_base_sparse_mat), intent(in), optional :: mold Integer :: err_act, nrow character(len=20) :: name='z_base_precbld' diff --git a/prec/psb_z_bjacprec.f03 b/prec/psb_z_bjacprec.f03 index a8913609..bdcaa581 100644 --- a/prec/psb_z_bjacprec.f03 +++ b/prec/psb_z_bjacprec.f03 @@ -4,8 +4,8 @@ module psb_z_bjacprec type, extends(psb_z_base_prec_type) :: psb_z_bjac_prec_type integer, allocatable :: iprcparm(:) - type(psb_zspmat_type), allocatable :: av(:) - complex(psb_dpk_), allocatable :: d(:) + type(psb_zspmat_type), allocatable :: av(:) + complex(psb_dpk_), allocatable :: d(:) contains procedure, pass(prec) :: apply => psb_z_bjac_apply procedure, pass(prec) :: precbld => psb_z_bjac_precbld @@ -215,17 +215,19 @@ contains end subroutine psb_z_bjac_precinit - subroutine psb_z_bjac_precbld(a,desc_a,prec,info,upd) + subroutine psb_z_bjac_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod use psb_prec_mod Implicit None type(psb_zspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_z_bjac_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_z_base_sparse_mat), intent(in), optional :: mold ! .. Local Scalars .. integer :: i, m @@ -325,12 +327,6 @@ contains goto 9999 end if -!!$ call prec%av(psb_l_pr_)%print(30+me) -!!$ call prec%av(psb_u_pr_)%print(40+me) -!!$ do i=1,n_row -!!$ write(50+me,*) i,prec%d(i) -!!$ end do - case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' @@ -344,6 +340,14 @@ contains goto 9999 end select + if (present(mold)) then + call prec%av(psb_l_pr_)%cscnv(info,mold=mold) + call prec%av(psb_u_pr_)%cscnv(info,mold=mold) + else if (present(afmt)) then + call prec%av(psb_l_pr_)%cscnv(info,type=afmt) + call prec%av(psb_u_pr_)%cscnv(info,type=afmt) + end if + call psb_erractionrestore(err_act) return diff --git a/prec/psb_z_diagprec.f03 b/prec/psb_z_diagprec.f03 index 71b20948..72cd34ac 100644 --- a/prec/psb_z_diagprec.f03 +++ b/prec/psb_z_diagprec.f03 @@ -153,16 +153,18 @@ contains end subroutine psb_z_diag_precinit - subroutine psb_z_diag_precbld(a,desc_a,prec,info,upd) + subroutine psb_z_diag_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod Implicit None type(psb_zspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_z_diag_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_z_base_sparse_mat), intent(in), optional :: mold Integer :: err_act, nrow,i character(len=20) :: name='z_diag_precbld' diff --git a/prec/psb_z_nullprec.f03 b/prec/psb_z_nullprec.f03 index a8c0e728..f4d7322a 100644 --- a/prec/psb_z_nullprec.f03 +++ b/prec/psb_z_nullprec.f03 @@ -104,16 +104,18 @@ contains return end subroutine psb_z_null_precinit - subroutine psb_z_null_precbld(a,desc_a,prec,info,upd) + subroutine psb_z_null_precbld(a,desc_a,prec,info,upd,mold,afmt) use psb_sparse_mod Implicit None type(psb_zspmat_type), intent(in), target :: a - type(psb_desc_type), intent(in), target :: desc_a + type(psb_desc_type), intent(in), target :: desc_a class(psb_z_null_prec_type),intent(inout) :: prec - integer, intent(out) :: info - character, intent(in), optional :: upd + integer, intent(out) :: info + character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_z_base_sparse_mat), intent(in), optional :: mold Integer :: err_act, nrow character(len=20) :: name='z_null_precbld' diff --git a/prec/psb_zprc_aply.f90 b/prec/psb_zprc_aply.f90 deleted file mode 100644 index 5033bd92..00000000 --- a/prec/psb_zprc_aply.f90 +++ /dev/null @@ -1,195 +0,0 @@ -!!$ -!!$ Parallel Sparse BLAS version 3.0 -!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 -!!$ Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari CNRS-IRIT, Toulouse -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the PSBLAS group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -subroutine psb_zprc_aply(prec,x,y,desc_data,info,trans, work) - - use psb_sparse_mod - use psb_prec_mod, psb_protect_name => psb_zprc_aply - - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_zprec_type), intent(in) :: prec - complex(psb_dpk_),intent(in) :: x(:) - complex(psb_dpk_),intent(inout) :: y(:) - integer, intent(out) :: info - character(len=1), optional :: trans - complex(psb_dpk_), intent(inout), optional, target :: work(:) - - ! Local variables - character :: trans_ - complex(psb_dpk_), pointer :: work_(:) - integer :: ictxt,np,me,err_act - character(len=20) :: name - - name='psb_prc_aply' - info = psb_success_ - call psb_erractionsave(err_act) - - ictxt=desc_data%matrix_data(psb_ctxt_) - call psb_info(ictxt, me, np) - - if (present(trans)) then - trans_=trans - else - trans_='N' - end if - - if (present(work)) then - work_ => work - else - allocate(work_(4*desc_data%matrix_data(psb_n_col_)),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - - end if - - if (.not.allocated(prec%prec)) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner") - goto 9999 - end if - call prec%prec%apply(zone,x,zzero,y,desc_data,info,trans_,work=work_) - if (present(work)) then - else - deallocate(work_) - end if - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return - -end subroutine psb_zprc_aply - - -!!$ -!!$ -!!$ MD2P4 -!!$ Multilevel Domain Decomposition Parallel Preconditioner Package for PSBLAS -!!$ for -!!$ Parallel Sparse BLAS version 3.0 -!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 -!!$ Salvatore Filippone University of Rome Tor Vergata -!!$ Alfredo Buttari CNRS-IRIT, Toulouse -!!$ Daniela di Serafino Second University of Naples -!!$ Pasqua D'Ambra ICAR-CNR -!!$ -!!$ Redistribution and use in source and binary forms, with or without -!!$ modification, are permitted provided that the following conditions -!!$ are met: -!!$ 1. Redistributions of source code must retain the above copyright -!!$ notice, this list of conditions and the following disclaimer. -!!$ 2. Redistributions in binary form must reproduce the above copyright -!!$ notice, this list of conditions, and the following disclaimer in the -!!$ documentation and/or other materials provided with the distribution. -!!$ 3. The name of the MD2P4 group or the names of its contributors may -!!$ not be used to endorse or promote products derived from this -!!$ software without specific written permission. -!!$ -!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS -!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED -!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MD2P4 GROUP OR ITS CONTRIBUTORS -!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR -!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF -!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS -!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN -!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) -!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -!!$ POSSIBILITY OF SUCH DAMAGE. -!!$ -!!$ -subroutine psb_zprc_aply1(prec,x,desc_data,info,trans) - - use psb_sparse_mod - use psb_prec_mod, psb_protect_name => psb_zprc_aply1 - implicit none - - type(psb_desc_type),intent(in) :: desc_data - type(psb_zprec_type), intent(in) :: prec - complex(psb_dpk_),intent(inout) :: x(:) - integer, intent(out) :: info - character(len=1), optional :: trans - - ! Local variables - character :: trans_ - integer :: ictxt,np,me, err_act - complex(psb_dpk_), pointer :: WW(:), w1(:) - character(len=20) :: name - name='psb_prc_aply1' - info = psb_success_ - call psb_erractionsave(err_act) - - - ictxt=desc_data%matrix_data(psb_ctxt_) - call psb_info(ictxt, me, np) - if (present(trans)) then - trans_=psb_toupper(trans) - else - trans_='N' - end if - - if (.not.allocated(prec%prec)) then - info = 1124 - call psb_errpush(info,name,a_err="preconditioner") - goto 9999 - end if - allocate(ww(size(x)),w1(size(x)),stat=info) - if (info /= psb_success_) then - call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 - end if - call prec%prec%apply(zone,x,zzero,ww,desc_data,info,trans_,work=w1) - if(info /= psb_success_) goto 9999 - x(:) = ww(:) - deallocate(ww,W1) - - call psb_erractionrestore(err_act) - return - -9999 continue - call psb_errpush(info,name) - call psb_erractionrestore(err_act) - if (err_act == psb_act_abort_) then - call psb_error() - return - end if - return -end subroutine psb_zprc_aply1 diff --git a/prec/psb_zprecbld.f90 b/prec/psb_zprecbld.f90 index ef4d80f1..61151fdb 100644 --- a/prec/psb_zprecbld.f90 +++ b/prec/psb_zprecbld.f90 @@ -29,7 +29,7 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -subroutine psb_zprecbld(a,desc_a,p,info,upd) +subroutine psb_zprecbld(a,desc_a,p,info,upd,mold,afmt) use psb_sparse_mod use psb_prec_mod, psb_protect_name => psb_zprecbld @@ -40,6 +40,8 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd) type(psb_zprec_type),intent(inout) :: p integer, intent(out) :: info character, intent(in), optional :: upd + character(len=*), intent(in), optional :: afmt + class(psb_z_base_sparse_mat), intent(in), optional :: mold ! Local scalars @@ -63,16 +65,6 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd) call psb_info(ictxt, me, np) - if (present(upd)) then - upd_ = psb_toupper(upd) - else - upd_='F' - endif - if ((upd_ == 'F').or.(upd_ == 'T')) then - ! ok - else - upd_='F' - endif n_row = psb_cd_get_local_rows(desc_a) n_col = psb_cd_get_local_cols(desc_a) mglob = psb_cd_get_global_rows(desc_a) @@ -88,7 +80,8 @@ subroutine psb_zprecbld(a,desc_a,p,info,upd) goto 9999 end if - call p%prec%precbld(a,desc_a,info,upd) + call p%prec%precbld(a,desc_a,info,upd=upd,afmt=afmt,mold=mold) + if (info /= psb_success_) goto 9999 call psb_erractionrestore(err_act) diff --git a/test/newfmt/Makefile b/test/newfmt/Makefile new file mode 100644 index 00000000..c8bc953e --- /dev/null +++ b/test/newfmt/Makefile @@ -0,0 +1,38 @@ +include ../../Make.inc +# +# Libraries used +# +LIBDIR=../../lib/ +OPTDIR=../../opt +PSBLAS_LIB= -L$(LIBDIR) -lpsb_util -lpsb_krylov -lpsb_prec -lpsb_base -L$(OPTDIR) -lopt +LDLIBS=$(PSBLDLIBS) +# +# Compilers and such +# +CCOPT= -g +FINCLUDES=$(FMFLAG)$(LIBDIR) $(FMFLAG). $(FMFLAG)$(OPTDIR) + + +EXEDIR=./runs + +all: ppde spde + +ppde: ppde.o + $(F90LINK) ppde.o -o ppde $(PSBLAS_LIB) $(LDLIBS) + /bin/mv ppde $(EXEDIR) + + +spde: spde.o + $(F90LINK) spde.o -o spde $(PSBLAS_LIB) $(LDLIBS) + /bin/mv spde $(EXEDIR) + + +clean: + /bin/rm -f ppde.o spde.o $(EXEDIR)/ppde +verycleanlib: + (cd ../..; make veryclean) +lib: + (cd ../../; make library) + + + diff --git a/test/newfmt/ppde.f90 b/test/newfmt/ppde.f90 new file mode 100644 index 00000000..3e7621b3 --- /dev/null +++ b/test/newfmt/ppde.f90 @@ -0,0 +1,671 @@ +!!$ +!!$ Parallel Sparse BLAS version 2.3.1 +!!$ (C) Copyright 2006, 2007, 2008, 2009, 2010 +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the PSBLAS group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! File: ppde.f90 +! +! Program: ppde +! This sample program solves a linear system obtained by discretizing a +! PDE with Dirichlet BCs. +! +! +! The PDE is a general second order equation in 3d +! +! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) +! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u = 0 +! dxdx dydy dzdz dx dy dz +! +! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. +! +! Example taken from: +! C.T.Kelley +! Iterative Methods for Linear and Nonlinear Equations +! SIAM 1995 +! +! In this sample program the index space of the discretized +! computational domain is first numbered sequentially in a standard way, +! then the corresponding vector is distributed according to a BLOCK +! data distribution. +! +! Boundary conditions are set in a very simple way, by adding +! equations of the form +! +! u(x,y) = exp(-x^2-y^2-z^2) +! +! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. +! +program ppde + use psb_sparse_mod + use psb_prec_mod + use psb_krylov_mod + use psb_d_ell_mat_mod + implicit none + + ! input parameters + character(len=20) :: kmethd, ptype + character(len=5) :: afmt + integer :: idim + + ! miscellaneous + real(psb_dpk_), parameter :: one = 1.d0 + real(psb_dpk_) :: t1, t2, tprec + + ! sparse matrix and preconditioner + type(psb_dspmat_type) :: a + type(psb_dprec_type) :: prec + type(psb_d_ell_sparse_mat) :: aell + ! descriptor + type(psb_desc_type) :: desc_a + ! dense matrices + real(psb_dpk_), allocatable :: b(:), x(:) + ! blacs parameters + integer :: ictxt, iam, np + + ! solver parameters + integer :: iter, itmax,itrace, istopc, irst + integer(psb_long_int_k_) :: amatsize, precsize, descsize + real(psb_dpk_) :: err, eps + + ! other variables + integer :: info, i + character(len=20) :: name,ch_err + + info=psb_success_ + + + call psb_init(ictxt) + call psb_info(ictxt,iam,np) + + if (iam < 0) then + ! This should not happen, but just in case + call psb_exit(ictxt) + stop + endif + if(psb_get_errstatus() /= 0) goto 9999 + name='pde90' + call psb_set_errverbosity(2) + ! + ! get parameters + ! + call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst) + + ! + ! allocate and fill in the coefficient matrix, rhs and initial guess + ! + call psb_barrier(ictxt) + t1 = psb_wtime() + call create_matrix(idim,a,b,x,desc_a,ictxt,afmt,info) + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='create_matrix' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + if (iam == psb_root_) write(psb_out_unit,'("Overall matrix creation time : ",es12.5)')t2 + if (iam == psb_root_) write(psb_out_unit,'(" ")') + ! + ! prepare the preconditioner. + ! + if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ",a)')ptype + call psb_precinit(prec,ptype,info) + + call psb_barrier(ictxt) + t1 = psb_wtime() + call psb_precbld(a,desc_a,prec,info,mold=aell) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_precbld' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + tprec = psb_wtime()-t1 + + call psb_amx(ictxt,tprec) + + if (iam == psb_root_) write(psb_out_unit,'("Preconditioner time : ",es12.5)')tprec + if (iam == psb_root_) write(psb_out_unit,'(" ")') + ! + ! iterative method parameters + ! + if(iam == psb_root_) write(psb_out_unit,'("Calling iterative method ",a)')kmethd + call psb_barrier(ictxt) + t1 = psb_wtime() + eps = 1.d-9 + call psb_krylov(kmethd,a,prec,b,x,eps,desc_a,info,& + & itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='solver routine' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + call psb_barrier(ictxt) + t2 = psb_wtime() - t1 + call psb_amx(ictxt,t2) + + amatsize = psb_sizeof(a) + descsize = psb_sizeof(desc_a) + precsize = psb_sizeof(prec) + call psb_sum(ictxt,amatsize) + call psb_sum(ictxt,descsize) + call psb_sum(ictxt,precsize) + + if (iam == psb_root_) then + write(psb_out_unit,'(" ")') + write(psb_out_unit,'("Time to solve matrix : ",es12.5)')t2 + write(psb_out_unit,'("Time per iteration : ",es12.5)')t2/iter + write(psb_out_unit,'("Number of iterations : ",i0)')iter + write(psb_out_unit,'("Convergence indicator on exit : ",es12.5)')err + write(psb_out_unit,'("Info on exit : ",i0)')info + write(psb_out_unit,'("Total memory occupation for A: ",i12)')amatsize + write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize + write(psb_out_unit,'("Total memory occupation for PREC: ",i12)')precsize + end if + + ! + ! cleanup storage and exit + ! + call psb_gefree(b,desc_a,info) + call psb_gefree(x,desc_a,info) + call psb_spfree(a,desc_a,info) + call psb_precfree(prec,info) + call psb_cdfree(desc_a,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='free routine' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + +9999 continue + if(info /= psb_success_) then + call psb_error(ictxt) + end if + call psb_exit(ictxt) + stop + +contains + ! + ! get iteration parameters from standard input + ! + subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst) + integer :: ictxt + character(len=*) :: kmethd, ptype, afmt + integer :: idim, istopc,itmax,itrace,irst + integer :: np, iam + integer :: intbuf(10), ip + + call psb_info(ictxt, iam, np) + + if (iam == 0) then + read(psb_inp_unit,*) ip + if (ip >= 3) then + read(psb_inp_unit,*) kmethd + read(psb_inp_unit,*) ptype + read(psb_inp_unit,*) afmt + + ! broadcast parameters to all processors + call psb_bcast(ictxt,kmethd) + call psb_bcast(ictxt,afmt) + call psb_bcast(ictxt,ptype) + + + read(psb_inp_unit,*) idim + if (ip >= 4) then + read(psb_inp_unit,*) istopc + else + istopc=1 + endif + if (ip >= 5) then + read(psb_inp_unit,*) itmax + else + itmax=500 + endif + if (ip >= 6) then + read(psb_inp_unit,*) itrace + else + itrace=-1 + endif + if (ip >= 7) then + read(psb_inp_unit,*) irst + else + irst=1 + endif + ! broadcast parameters to all processors + + intbuf(1) = idim + intbuf(2) = istopc + intbuf(3) = itmax + intbuf(4) = itrace + intbuf(5) = irst + call psb_bcast(ictxt,intbuf(1:5)) + + write(psb_out_unit,'("Solving matrix : ell1")') + write(psb_out_unit,'("Grid dimensions : ",i4,"x",i4,"x",i4)')idim,idim,idim + write(psb_out_unit,'("Number of processors : ",i0)')np + write(psb_out_unit,'("Data distribution : BLOCK")') + write(psb_out_unit,'("Preconditioner : ",a)') ptype + write(psb_out_unit,'("Iterative method : ",a)') kmethd + write(psb_out_unit,'(" ")') + else + ! wrong number of parameter, print an error message and exit + call pr_usage(0) + call psb_abort(ictxt) + stop 1 + endif + else + call psb_bcast(ictxt,kmethd) + call psb_bcast(ictxt,afmt) + call psb_bcast(ictxt,ptype) + call psb_bcast(ictxt,intbuf(1:5)) + idim = intbuf(1) + istopc = intbuf(2) + itmax = intbuf(3) + itrace = intbuf(4) + irst = intbuf(5) + end if + return + + end subroutine get_parms + ! + ! print an error message + ! + subroutine pr_usage(iout) + integer :: iout + write(iout,*)'incorrect parameter(s) found' + write(iout,*)' usage: pde90 methd prec dim & + &[istop itmax itrace]' + write(iout,*)' where:' + write(iout,*)' methd: cgstab cgs rgmres bicgstabl' + write(iout,*)' prec : bjac diag none' + write(iout,*)' dim number of points along each axis' + write(iout,*)' the size of the resulting linear ' + write(iout,*)' system is dim**3' + write(iout,*)' istop stopping criterion 1, 2 ' + write(iout,*)' itmax maximum number of iterations [500] ' + write(iout,*)' itrace <=0 (no tracing, default) or ' + write(iout,*)' >= 1 do tracing every itrace' + write(iout,*)' iterations ' + end subroutine pr_usage + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine create_matrix(idim,a,b,xv,desc_a,ictxt,afmt,info) + ! + ! discretize the partial diferential equation + ! + ! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) + ! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u + ! dxdx dydy dzdz dx dy dz + ! + ! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. + ! + ! Boundary conditions are set in a very simple way, by adding + ! equations of the form + ! + ! u(x,y) = exp(-x^2-y^2-z^2) + ! + ! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. + ! + use psb_sparse_mod + use psb_d_ell_mat_mod + implicit none + integer :: idim + integer, parameter :: nb=20 + real(psb_dpk_), allocatable :: b(:),xv(:) + type(psb_desc_type) :: desc_a + integer :: ictxt, info + character :: afmt*5 + type(psb_dspmat_type) :: a + type(psb_d_csc_sparse_mat) :: acsc + type(psb_d_coo_sparse_mat) :: acoo + type(psb_d_csr_sparse_mat) :: acsr + type(psb_d_ell_sparse_mat) :: aell + real(psb_dpk_) :: zt(nb),x,y,z + integer :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer :: ix,iy,iz,ia,indx_owner + integer :: np, iam, nr, nt + integer :: element + integer, allocatable :: irow(:),icol(:),myidx(:) + real(psb_dpk_), allocatable :: val(:) + ! deltah dimension of each grid cell + + ! deltat discretization time + real(psb_dpk_) :: deltah + real(psb_dpk_),parameter :: rhs=0.d0,one=1.d0,zero=0.d0 + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen + real(psb_dpk_) :: a1, a2, a3, a4, b1, b2, b3 + external :: a1, a2, a3, a4, b1, b2, b3 + integer :: err_act + + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + deltah = 1.d0/(idim-1) + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim*idim + n = m + nnz = ((n*9)/(np)) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + call psb_barrier(ictxt) + t0 = psb_wtime() + call psb_cdall(ictxt,desc_a,info,nl=nr) + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) + ! define rhs from boundary conditions; also build initial guess + if (info == psb_success_) call psb_geall(b,desc_a,info) + if (info == psb_success_) call psb_geall(xv,desc_a,info) + nlr = psb_cd_get_local_rows(desc_a) + call psb_barrier(ictxt) + talc = psb_wtime()-t0 + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! we build an auxiliary matrix consisting of one row at a + ! time; just a small matrix. might be extended to generate + ! a bunch of rows per call. + ! + allocate(val(20*nb),irow(20*nb),& + &icol(20*nb),myidx(nlr),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + + do i=1,nlr + myidx(i) = i + end do + + + call psb_loc_to_glob(myidx,desc_a,info) + + ! loop over rows belonging to current process in a block + ! distribution. + + call psb_barrier(ictxt) + t1 = psb_wtime() + do ii=1, nlr,nb + ib = min(nb,nlr-ii+1) + element = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + if (mod(glob_row,(idim*idim)) == 0) then + ix = glob_row/(idim*idim) + else + ix = glob_row/(idim*idim)+1 + endif + if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then + iy = (glob_row-(ix-1)*idim*idim)/idim + else + iy = (glob_row-(ix-1)*idim*idim)/idim+1 + endif + iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim + ! x, y, x coordinates + x = ix*deltah + y = iy*deltah + z = iz*deltah + + ! check on boundary points + zt(k) = 0.d0 + ! internal point: build discretization + ! + ! term depending on (x-1,y,z) + ! + if (ix == 1) then + val(element)=-b1(x,y,z)-a1(x,y,z) + val(element) = val(element)/(deltah*& + & deltah) + zt(k) = exp(-y**2-z**2)*(-val(element)) + else + val(element)=-b1(x,y,z)-a1(x,y,z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element) = (ix-2)*idim*idim+(iy-1)*idim+(iz) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x,y-1,z) + if (iy == 1) then + val(element)=-b2(x,y,z)-a2(x,y,z) + val(element) = val(element)/(deltah*& + & deltah) + zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) + else + val(element)=-b2(x,y,z)-a2(x,y,z) + val(element) = val(element)/(deltah*deltah) + icol(element) = (ix-1)*idim*idim+(iy-2)*idim+(iz) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x,y,z-1) + if (iz == 1) then + val(element)=-b3(x,y,z)-a3(x,y,z) + val(element) = val(element)/(deltah*deltah) + zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) + else + val(element)=-b3(x,y,z)-a3(x,y,z) + val(element) = val(element)/(deltah*deltah) + icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x,y,z) + val(element)=2*b1(x,y,z) + 2*b2(x,y,z)& + & + 2*b3(x,y,z) + a1(x,y,z)& + & + a2(x,y,z) + a3(x,y,z) + val(element) = val(element)/(deltah*deltah) + icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz) + irow(element) = glob_row + element = element+1 + ! term depending on (x,y,z+1) + if (iz == idim) then + val(element)=-b1(x,y,z) + val(element) = val(element)/(deltah*deltah) + zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) + else + val(element)=-b1(x,y,z) + val(element) = val(element)/(deltah*deltah) + icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x,y+1,z) + if (iy == idim) then + val(element)=-b2(x,y,z) + val(element) = val(element)/(deltah*deltah) + zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) + else + val(element)=-b2(x,y,z) + val(element) = val(element)/(deltah*deltah) + icol(element) = (ix-1)*idim*idim+(iy)*idim+(iz) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x+1,y,z) + if (ix= 3) then + read(psb_inp_unit,*) kmethd + read(psb_inp_unit,*) ptype + read(psb_inp_unit,*) afmt + + ! broadcast parameters to all processors + call psb_bcast(ictxt,kmethd) + call psb_bcast(ictxt,afmt) + call psb_bcast(ictxt,ptype) + + + read(psb_inp_unit,*) idim + if (ip >= 4) then + read(psb_inp_unit,*) istopc + else + istopc=1 + endif + if (ip >= 5) then + read(psb_inp_unit,*) itmax + else + itmax=500 + endif + if (ip >= 6) then + read(psb_inp_unit,*) itrace + else + itrace=-1 + endif + if (ip >= 7) then + read(psb_inp_unit,*) irst + else + irst=1 + endif + ! broadcast parameters to all processors + + intbuf(1) = idim + intbuf(2) = istopc + intbuf(3) = itmax + intbuf(4) = itrace + intbuf(5) = irst + call psb_bcast(ictxt,intbuf(1:5)) + + write(psb_out_unit,'("Solving matrix : ell1")') + write(psb_out_unit,'("Grid dimensions : ",i4,"x",i4,"x",i4)')idim,idim,idim + write(psb_out_unit,'("Number of processors : ",i0)')np + write(psb_out_unit,'("Data distribution : BLOCK")') + write(psb_out_unit,'("Preconditioner : ",a)') ptype + write(psb_out_unit,'("Iterative method : ",a)') kmethd + write(psb_out_unit,'(" ")') + else + ! wrong number of parameter, print an error message and exit + call pr_usage(0) + call psb_abort(ictxt) + stop 1 + endif + else + call psb_bcast(ictxt,kmethd) + call psb_bcast(ictxt,afmt) + call psb_bcast(ictxt,ptype) + call psb_bcast(ictxt,intbuf(1:5)) + idim = intbuf(1) + istopc = intbuf(2) + itmax = intbuf(3) + itrace = intbuf(4) + irst = intbuf(5) + end if + return + + end subroutine get_parms + ! + ! print an error message + ! + subroutine pr_usage(iout) + integer :: iout + write(iout,*)'incorrect parameter(s) found' + write(iout,*)' usage: pde90 methd prec dim & + &[istop itmax itrace]' + write(iout,*)' where:' + write(iout,*)' methd: cgstab cgs rgmres bicgstabl' + write(iout,*)' prec : bjac diag none' + write(iout,*)' dim number of points along each axis' + write(iout,*)' the size of the resulting linear ' + write(iout,*)' system is dim**3' + write(iout,*)' istop stopping criterion 1, 2 ' + write(iout,*)' itmax maximum number of iterations [500] ' + write(iout,*)' itrace <=0 (no tracing, default) or ' + write(iout,*)' >= 1 do tracing every itrace' + write(iout,*)' iterations ' + end subroutine pr_usage + + ! + ! subroutine to allocate and fill in the coefficient matrix and + ! the rhs. + ! + subroutine create_matrix(idim,a,b,xv,desc_a,ictxt,afmt,info) + ! + ! discretize the partial diferential equation + ! + ! b1 dd(u) b2 dd(u) b3 dd(u) a1 d(u) a2 d(u) a3 d(u) + ! - ------ - ------ - ------ - ----- - ------ - ------ + a4 u + ! dxdx dydy dzdz dx dy dz + ! + ! with Dirichlet boundary conditions, on the unit cube 0<=x,y,z<=1. + ! + ! Boundary conditions are set in a very simple way, by adding + ! equations of the form + ! + ! u(x,y) = exp(-x^2-y^2-z^2) + ! + ! Note that if a1=a2=a3=a4=0., the PDE is the well-known Laplace equation. + ! + use psb_sparse_mod + use psb_mat_mod + implicit none + integer :: idim + integer, parameter :: nb=20 + real(psb_spk_), allocatable :: b(:),xv(:) + type(psb_desc_type) :: desc_a + integer :: ictxt, info + character :: afmt*5 + type(psb_sspmat_type) :: a + type(psb_s_coo_sparse_mat) :: acoo + type(psb_s_csr_sparse_mat) :: acsr + real(psb_spk_) :: zt(nb),x,y,z + integer :: m,n,nnz,glob_row,nlr,i,ii,ib,k + integer :: ix,iy,iz,ia,indx_owner + integer :: np, iam, nr, nt + integer :: element + integer, allocatable :: irow(:),icol(:),myidx(:) + real(psb_spk_), allocatable :: val(:) + ! deltah dimension of each grid cell + ! deltat discretization time + real(psb_spk_) :: deltah + real(psb_spk_),parameter :: rhs=0.0,one=1.0,zero=0.0 + real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen + real(psb_spk_) :: a1, a2, a3, a4, b1, b2, b3 + external :: a1, a2, a3, a4, b1, b2, b3 + integer :: err_act + + character(len=20) :: name, ch_err,tmpfmt + + info = psb_success_ + name = 'create_matrix' + call psb_erractionsave(err_act) + + call psb_info(ictxt, iam, np) + + deltah = 1.d0/(idim-1) + + ! initialize array descriptor and sparse matrix storage. provide an + ! estimate of the number of non zeroes + + m = idim*idim*idim + n = m + nnz = ((n*9)/(np)) + if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n + + ! + ! Using a simple BLOCK distribution. + ! + nt = (m+np-1)/np + nr = max(0,min(nt,m-(iam*nt))) + + nt = nr + call psb_sum(ictxt,nt) + if (nt /= m) write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m + call psb_barrier(ictxt) + t0 = psb_wtime() + call psb_cdall(ictxt,desc_a,info,nl=nr) + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) + ! define rhs from boundary conditions; also build initial guess + if (info == psb_success_) call psb_geall(b,desc_a,info) + if (info == psb_success_) call psb_geall(xv,desc_a,info) + nlr = psb_cd_get_local_rows(desc_a) + call psb_barrier(ictxt) + talc = psb_wtime()-t0 + + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='allocation rout.' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + ! we build an auxiliary matrix consisting of one row at a + ! time; just a small matrix. might be extended to generate + ! a bunch of rows per call. + ! + allocate(val(20*nb),irow(20*nb),& + &icol(20*nb),myidx(nlr),stat=info) + if (info /= psb_success_ ) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + endif + + do i=1,nlr + myidx(i) = i + end do + + + call psb_loc_to_glob(myidx,desc_a,info) + + ! loop over rows belonging to current process in a block + ! distribution. + + call psb_barrier(ictxt) + t1 = psb_wtime() + do ii=1, nlr,nb + ib = min(nb,nlr-ii+1) + element = 1 + do k=1,ib + i=ii+k-1 + ! local matrix pointer + glob_row=myidx(i) + ! compute gridpoint coordinates + if (mod(glob_row,(idim*idim)) == 0) then + ix = glob_row/(idim*idim) + else + ix = glob_row/(idim*idim)+1 + endif + if (mod((glob_row-(ix-1)*idim*idim),idim) == 0) then + iy = (glob_row-(ix-1)*idim*idim)/idim + else + iy = (glob_row-(ix-1)*idim*idim)/idim+1 + endif + iz = glob_row-(ix-1)*idim*idim-(iy-1)*idim + ! x, y, x coordinates + x = ix*deltah + y = iy*deltah + z = iz*deltah + + ! check on boundary points + zt(k) = 0.d0 + ! internal point: build discretization + ! + ! term depending on (x-1,y,z) + ! + if (ix == 1) then + val(element)=-b1(x,y,z)-a1(x,y,z) + val(element) = val(element)/(deltah*& + & deltah) + zt(k) = exp(-y**2-z**2)*(-val(element)) + else + val(element)=-b1(x,y,z)-a1(x,y,z) + val(element) = val(element)/(deltah*& + & deltah) + icol(element) = (ix-2)*idim*idim+(iy-1)*idim+(iz) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x,y-1,z) + if (iy == 1) then + val(element)=-b2(x,y,z)-a2(x,y,z) + val(element) = val(element)/(deltah*& + & deltah) + zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) + else + val(element)=-b2(x,y,z)-a2(x,y,z) + val(element) = val(element)/(deltah*deltah) + icol(element) = (ix-1)*idim*idim+(iy-2)*idim+(iz) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x,y,z-1) + if (iz == 1) then + val(element)=-b3(x,y,z)-a3(x,y,z) + val(element) = val(element)/(deltah*deltah) + zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) + else + val(element)=-b3(x,y,z)-a3(x,y,z) + val(element) = val(element)/(deltah*deltah) + icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x,y,z) + val(element)=2*b1(x,y,z) + 2*b2(x,y,z)& + & + 2*b3(x,y,z) + a1(x,y,z)& + & + a2(x,y,z) + a3(x,y,z) + val(element) = val(element)/(deltah*deltah) + icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz) + irow(element) = glob_row + element = element+1 + ! term depending on (x,y,z+1) + if (iz == idim) then + val(element)=-b1(x,y,z) + val(element) = val(element)/(deltah*deltah) + zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) + else + val(element)=-b1(x,y,z) + val(element) = val(element)/(deltah*deltah) + icol(element) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x,y+1,z) + if (iy == idim) then + val(element)=-b2(x,y,z) + val(element) = val(element)/(deltah*deltah) + zt(k) = exp(-y**2-z**2)*exp(-x)*(-val(element)) + else + val(element)=-b2(x,y,z) + val(element) = val(element)/(deltah*deltah) + icol(element) = (ix-1)*idim*idim+(iy)*idim+(iz) + irow(element) = glob_row + element = element+1 + endif + ! term depending on (x+1,y,z) + if (ix