diff --git a/test/idx/Makefile b/test/idx/Makefile index f6179d09..7e8a9fcc 100644 --- a/test/idx/Makefile +++ b/test/idx/Makefile @@ -14,25 +14,16 @@ CCOPT= -g FINCLUDES=$(FMFLAG)$(MODDIR) $(FMFLAG). -EXEDIR=./runs +EXEDIR=./ -all: tryidxijk psb_d_pde3d test_gf731 +all: tryidxijk tryidxijk: tryidxijk.o $(FLINK) tryidxijk.o -o tryidxijk $(PSBLAS_LIB) $(LDLIBS) /bin/mv tryidxijk $(EXEDIR) -psb_d_pde3d: psb_d_pde3d.o - $(FLINK) psb_d_pde3d.o -o psb_d_pde3d $(PSBLAS_LIB) $(LDLIBS) - /bin/mv psb_d_pde3d $(EXEDIR) -test_gf731: test_gf731.o - $(FLINK) test_gf731.o -o test_gf731 $(PSBLAS_LIB) $(LDLIBS) - /bin/mv test_gf731 $(EXEDIR) - - - clean: - /bin/rm -f tryidxijk.o test_gf731.o psb_d_pde3d.o *$(.mod) \ + /bin/rm -f tryidxijk.o *$(.mod) \ $(EXEDIR)/tryidxijk $(EXEDIR)/psb_d_pde3d verycleanlib: (cd ../..; make veryclean) diff --git a/test/idx/psb_d_pde3d.f90 b/test/idx/psb_d_pde3d.f90 deleted file mode 100644 index 885bf00a..00000000 --- a/test/idx/psb_d_pde3d.f90 +++ /dev/null @@ -1,843 +0,0 @@ -! -! Parallel Sparse BLAS version 3.5 -! (C) Copyright 2006-2018 -! 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: -! 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: psb_d_pde3d.f90 -! -! Program: psb_d_pde3d -! 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 -! -! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) -! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f -! dxdx dydy dzdz dx dy dz -! -! with Dirichlet boundary conditions -! u = g -! -! on the unit cube 0<=x,y,z<=1. -! -! -! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. -! -! There are three choices available for data distribution: -! 1. A simple BLOCK distribution -! 2. A ditribution based on arbitrary assignment of indices to processes, -! typically from a graph partitioner -! 3. A 3D distribution in which the unit cube is partitioned -! into subcubes, each one assigned to a process. -! -! -module psb_d_pde3d_mod - - - use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_desc_type,& - & psb_dspmat_type, psb_d_vect_type, dzero,& - & psb_d_base_sparse_mat, psb_d_base_vect_type, psb_i_base_vect_type - - interface - function d_func_3d(x,y,z) result(val) - import :: psb_dpk_ - real(psb_dpk_), intent(in) :: x,y,z - real(psb_dpk_) :: val - end function d_func_3d - end interface - - interface psb_gen_pde3d - module procedure psb_d_gen_pde3d - end interface psb_gen_pde3d - -contains - - function d_null_func_3d(x,y,z) result(val) - - real(psb_dpk_), intent(in) :: x,y,z - real(psb_dpk_) :: val - - val = dzero - - end function d_null_func_3d - - - ! - ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. - ! - subroutine psb_d_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,& - & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,partition,nrl,iv) - use psb_base_mod - use psb_util_mod - ! - ! Discretizes the partial differential equation - ! - ! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) - ! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f - ! dxdx dydy dzdz dx dy dz - ! - ! with Dirichlet boundary conditions - ! u = g - ! - ! on the unit cube 0<=x,y,z<=1. - ! - ! - ! Note that if b1=b2=b3=c=0., the PDE is the Laplace equation. - ! - implicit none - procedure(d_func_3d) :: b1,b2,b3,c,a1,a2,a3,g - integer(psb_ipk_) :: idim - type(psb_dspmat_type) :: a - type(psb_d_vect_type) :: xv,bv - type(psb_desc_type) :: desc_a - integer(psb_ipk_) :: ictxt, info - character(len=*) :: afmt - procedure(d_func_3d), optional :: f - class(psb_d_base_sparse_mat), optional :: amold - class(psb_d_base_vect_type), optional :: vmold - class(psb_i_base_vect_type), optional :: imold - integer(psb_ipk_), optional :: partition, nrl,iv(:) - - ! Local variables. - - integer(psb_ipk_), parameter :: nb=20 - type(psb_d_csc_sparse_mat) :: acsc - type(psb_d_coo_sparse_mat) :: acoo - type(psb_d_csr_sparse_mat) :: acsr - real(psb_dpk_) :: zt(nb),x,y,z - integer(psb_ipk_) :: m,n,nnz,nr,nt,glob_row,nlr,i,j,ii,ib,k, partition_ - integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner - ! For 3D partition - integer(psb_ipk_) :: npx,npy,npz, npdims(3),iamx,iamy,iamz,mynx,myny,mynz - integer(psb_ipk_), allocatable :: bndx(:),bndy(:),bndz(:) - ! Process grid - integer(psb_ipk_) :: np, iam - integer(psb_ipk_) :: icoeff - integer(psb_lpk_), allocatable :: myidx(:) - integer(psb_ipk_), allocatable :: irow(:),icol(:) - real(psb_dpk_), allocatable :: val(:) - ! deltah dimension of each grid cell - ! deltat discretization time - real(psb_dpk_) :: deltah, sqdeltah, deltah2 - real(psb_dpk_), parameter :: rhs=dzero,one=done,zero=dzero - real(psb_dpk_) :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb - integer(psb_ipk_) :: err_act - procedure(d_func_3d), pointer :: f_ - character(len=20) :: name, ch_err,tmpfmt - - info = psb_success_ - name = 'create_matrix' - call psb_erractionsave(err_act) - - call psb_info(ictxt, iam, np) - - - if (present(f)) then - f_ => f - else - f_ => d_null_func_3d - end if - - deltah = done/(idim+2) - sqdeltah = deltah*deltah - deltah2 = (2*done)* deltah - - if (present(partition)) then - if ((1<= partition).and.(partition <= 3)) then - partition_ = partition - else - write(*,*) 'Invalid partition choice ',partition,' defaulting to 3' - partition_ = 3 - end if - else - partition_ = 3 - end if - - ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - t0 = psb_wtime() - m = idim*idim*idim - n = m - nnz = ((n*9)/(np)) - if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n - - select case(partition_) - case(1) - write(*,*) 'BLOCK partition ' - ! A BLOCK partition - if (present(nrl)) then - nr = nrl - else - ! - ! Using a simple BLOCK distribution. - ! - nt = (m+np-1)/np - nr = max(0,min(nt,m-(iam*nt))) - end if - - nt = nr - call psb_sum(ictxt,nt) - if (nt /= m) then - write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m - info = -1 - call psb_barrier(ictxt) - call psb_abort(ictxt) - return - end if - - ! - ! First example of use of CDALL: specify for each process a number of - ! contiguous rows - ! - call psb_cdall(ictxt,desc_a,info,nl=nr) - myidx = desc_a%get_global_indices() - nlr = size(myidx) - - case(2) - write(*,*) 'User Defined Partition' - ! A partition defined by the user through IV - - if (present(iv)) then - if (size(iv) /= m) then - write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m - info = -1 - call psb_barrier(ictxt) - call psb_abort(ictxt) - return - end if - else - write(psb_err_unit,*) iam, 'Initialization error: IV not present' - info = -1 - call psb_barrier(ictxt) - call psb_abort(ictxt) - return - end if - - ! - ! Second example of use of CDALL: specify for each row the - ! process that owns it - ! - call psb_cdall(ictxt,desc_a,info,vg=iv) - myidx = desc_a%get_global_indices() - nlr = size(myidx) - - case(3) - write(*,*) '3D coordinate planes Partition' - ! A 3-dimensional partition - - ! A nifty MPI function will split the process list - npdims = 0 - call mpi_dims_create(np,3,npdims,info) - npx = npdims(1) - npy = npdims(2) - npz = npdims(3) - - allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz)) - ! We can reuse idx2ijk for process indices as well. - call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0) - ! Now let's split the 3D cube in hexahedra - call dist1Didx(bndx,idim,npx) - mynx = bndx(iamx+1)-bndx(iamx) - call dist1Didx(bndy,idim,npy) - myny = bndy(iamy+1)-bndy(iamy) - call dist1Didx(bndz,idim,npz) - mynz = bndz(iamz+1)-bndz(iamz) - - ! How many indices do I own? - nlr = mynx*myny*mynz - allocate(myidx(nlr)) - ! Now, let's generate the list of indices I own - nr = 0 - do i=bndx(iamx),bndx(iamx+1)-1 - do j=bndy(iamy),bndy(iamy+1)-1 - do k=bndz(iamz),bndz(iamz+1)-1 - nr = nr + 1 - call ijk2idx(myidx(nr),i,j,k,idim,idim,idim) - end do - end do - end do - if (nr /= nlr) then - write(psb_err_unit,*) iam,iamx,iamy,iamz, 'Initialization error: NR vs NLR ',& - & nr,nlr,mynx,myny,mynz - info = -1 - call psb_barrier(ictxt) - call psb_abort(ictxt) - end if - - ! - ! Third example of use of CDALL: specify for each process - ! the set of global indices it owns. - ! - call psb_cdall(ictxt,desc_a,info,vll=myidx) - - case default - write(psb_err_unit,*) iam, 'Initialization error: should not get here' - info = -1 - call psb_barrier(ictxt) - call psb_abort(ictxt) - return - end select - - 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(xv,desc_a,info) - if (info == psb_success_) call psb_geall(bv,desc_a,info) - - 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),stat=info) - if (info /= psb_success_ ) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - endif - - - ! 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) - icoeff = 1 - do k=1,ib - i=ii+k-1 - ! local matrix pointer - glob_row=myidx(i) - ! compute gridpoint coordinates - call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim) - ! x, y, z coordinates - x = (ix-1)*deltah - y = (iy-1)*deltah - z = (iz-1)*deltah - zt(k) = f_(x,y,z) -!!$ write(*,*) 'idx2ijk ',ix,iy,iz,glob_row,x,y,z,a1(x,y,z) -!!$ return - ! internal point: build discretization - ! - ! term depending on (x-1,y,z) - ! - val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2 - if (ix == 1) then - zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k) - else - icol(icoeff) = (ix-2)*idim*idim+(iy-1)*idim+(iz) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - ! term depending on (x,y-1,z) - val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2 - if (iy == 1) then - zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k) - else - icol(icoeff) = (ix-1)*idim*idim+(iy-2)*idim+(iz) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - ! term depending on (x,y,z-1) - val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2 - if (iz == 1) then - zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k) - else - icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz-1) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - - ! term depending on (x,y,z) - val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & - & + c(x,y,z) - icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz) - irow(icoeff) = glob_row - icoeff = icoeff+1 - ! term depending on (x,y,z+1) - val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2 - if (iz == idim) then - zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k) - else - icol(icoeff) = (ix-1)*idim*idim+(iy-1)*idim+(iz+1) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - ! term depending on (x,y+1,z) - val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2 - if (iy == idim) then - zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k) - else - icol(icoeff) = (ix-1)*idim*idim+(iy)*idim+(iz) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - ! term depending on (x+1,y,z) - val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2 - if (ix==idim) then - zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k) - else - icol(icoeff) = (ix)*idim*idim+(iy-1)*idim+(iz) - irow(icoeff) = glob_row - icoeff = icoeff+1 - endif - - end do - call psb_spins(icoeff-1,irow,icol,val,a,desc_a,info) - if(info /= psb_success_) exit - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),bv,desc_a,info) - if(info /= psb_success_) exit - zt(:)=dzero - call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info) - if(info /= psb_success_) exit - end do - - tgen = psb_wtime()-t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='insert rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - deallocate(val,irow,icol) - - call psb_barrier(ictxt) - t1 = psb_wtime() - call psb_cdasb(desc_a,info,mold=imold) - tcdasb = psb_wtime()-t1 - call psb_barrier(ictxt) - t1 = psb_wtime() - if (info == psb_success_) then - if (present(amold)) then - call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) - else - call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) - end if - end if - call psb_barrier(ictxt) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='asb rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - if (info == psb_success_) call psb_geasb(xv,desc_a,info,mold=vmold) - if (info == psb_success_) call psb_geasb(bv,desc_a,info,mold=vmold) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='asb rout.' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - tasb = psb_wtime()-t1 - call psb_barrier(ictxt) - ttot = psb_wtime() - t0 - - call psb_amx(ictxt,talc) - call psb_amx(ictxt,tgen) - call psb_amx(ictxt,tasb) - call psb_amx(ictxt,ttot) - if(iam == psb_root_) then - tmpfmt = a%get_fmt() - write(psb_out_unit,'("The matrix has been generated and assembled in ",a3," format.")')& - & tmpfmt - write(psb_out_unit,'("-allocation time : ",es12.5)') talc - write(psb_out_unit,'("-coeff. gen. time : ",es12.5)') tgen - write(psb_out_unit,'("-desc asbly time : ",es12.5)') tcdasb - write(psb_out_unit,'("- mat asbly time : ",es12.5)') tasb - write(psb_out_unit,'("-total time : ",es12.5)') ttot - - end if - call psb_erractionrestore(err_act) - return - -9999 call psb_error_handler(ictxt,err_act) - - return - end subroutine psb_d_gen_pde3d - - -end module psb_d_pde3d_mod - -program psb_d_pde3d - use psb_base_mod - use psb_prec_mod - use psb_krylov_mod - use psb_util_mod - use psb_d_pde3d_mod - implicit none - - ! input parameters - character(len=20) :: kmethd, ptype - character(len=5) :: afmt - integer(psb_ipk_) :: idim - - ! miscellaneous - real(psb_dpk_), parameter :: one = done - real(psb_dpk_) :: t1, t2, tprec - - ! sparse matrix and preconditioner - type(psb_dspmat_type) :: a - type(psb_dprec_type) :: prec - ! descriptor - type(psb_desc_type) :: desc_a - ! dense vectors - type(psb_d_vect_type) :: xxv,bv - ! parallel environment - integer(psb_ipk_) :: ictxt, iam, np - - ! solver parameters - integer(psb_ipk_) :: iter, itmax,itrace, istopc, irst - integer(psb_epk_) :: amatsize, precsize, descsize, d2size - real(psb_dpk_) :: err, eps - - ! other variables - integer(psb_ipk_) :: info, i - character(len=20) :: name,ch_err - character(len=40) :: fname - - 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='pde3d90' - call psb_set_errverbosity(itwo) - call psb_cd_set_large_threshold(itwo) - ! - ! Hello world - ! - if (iam == psb_root_) then - write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ - write(*,*) 'This is the ',trim(name),' sample program' - end if - ! - ! 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 - ! - !write(*,*) 'Check a1:',a1(dzero,dzero,dzero) - call psb_barrier(ictxt) - t1 = psb_wtime() - call psb_gen_pde3d(ictxt,idim,a,bv,xxv,desc_a,afmt,& - & a1,a2,a3,b1,b2,b3,c,g,info,partition=1) - call psb_barrier(ictxt) - t2 = psb_wtime() - t1 - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_gen_pde3d' - 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 prec%init(ptype,info) - - call psb_barrier(ictxt) - t1 = psb_wtime() - call prec%build(a,desc_a,info) - 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,'(" ")') - call prec%descr() - ! - ! 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,bv,xxv,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 = a%sizeof() - descsize = desc_a%sizeof() - precsize = prec%sizeof() - 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,'("Number of processes : ",i0)')np - write(psb_out_unit,'("Time to solve system : ",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 PREC: ",i12)')precsize - write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize - write(psb_out_unit,'("Storage type for DESC_A: ",a)') desc_a%get_fmt() - end if - - - ! - ! cleanup storage and exit - ! - call psb_gefree(bv,desc_a,info) - call psb_gefree(xxv,desc_a,info) - call psb_spfree(a,desc_a,info) - call prec%free(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 - - call psb_exit(ictxt) - stop - -9999 call psb_error(ictxt) - - stop - -contains - ! - ! get iteration parameters from standard input - ! - subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst) - integer(psb_ipk_) :: ictxt - character(len=*) :: kmethd, ptype, afmt - integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst - integer(psb_ipk_) :: np, iam - integer(psb_ipk_) :: ip, inp_unit - character(len=1024) :: filename - - call psb_info(ictxt, iam, np) - - if (iam == 0) then - if (command_argument_count()>0) then - call get_command_argument(1,filename) - inp_unit = 30 - open(inp_unit,file=filename,action='read',iostat=info) - if (info /= 0) then - write(psb_err_unit,*) 'Could not open file ',filename,' for input' - call psb_abort(ictxt) - stop - else - write(psb_err_unit,*) 'Opened file ',trim(filename),' for input' - end if - else - inp_unit=psb_inp_unit - end if - read(inp_unit,*) ip - if (ip >= 3) then - read(inp_unit,*) kmethd - read(inp_unit,*) ptype - read(inp_unit,*) afmt - - read(inp_unit,*) idim - if (ip >= 4) then - read(inp_unit,*) istopc - else - istopc=1 - endif - if (ip >= 5) then - read(inp_unit,*) itmax - else - itmax=500 - endif - if (ip >= 6) then - read(inp_unit,*) itrace - else - itrace=-1 - endif - if (ip >= 7) then - read(inp_unit,*) irst - else - irst=1 - endif - ! broadcast parameters to all processors - - - 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(izero) - call psb_abort(ictxt) - stop 1 - endif - if (inp_unit /= psb_inp_unit) then - close(inp_unit) - end if - - end if - ! broadcast parameters to all processors - call psb_bcast(ictxt,kmethd) - call psb_bcast(ictxt,afmt) - call psb_bcast(ictxt,ptype) - call psb_bcast(ictxt,idim) - call psb_bcast(ictxt,istopc) - call psb_bcast(ictxt,itmax) - call psb_bcast(ictxt,itrace) - call psb_bcast(ictxt,irst) - - return - - end subroutine get_parms - ! - ! print an error message - ! - subroutine pr_usage(iout) - integer(psb_ipk_) :: iout - write(iout,*)'incorrect parameter(s) found' - write(iout,*)' usage: pde3d90 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 - ! - ! functions parametrizing the differential equation - ! - function b1(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: b1 - real(psb_dpk_), intent(in) :: x,y,z - b1=done/sqrt((3*done)) - end function b1 - function b2(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: b2 - real(psb_dpk_), intent(in) :: x,y,z - b2=done/sqrt((3*done)) - end function b2 - function b3(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: b3 - real(psb_dpk_), intent(in) :: x,y,z - b3=done/sqrt((3*done)) - end function b3 - function c(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: c - real(psb_dpk_), intent(in) :: x,y,z - c=dzero - end function c - function a1(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a1 - real(psb_dpk_), intent(in) :: x,y,z - a1=done/80 - end function a1 - function a2(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a2 - real(psb_dpk_), intent(in) :: x,y,z - a2=done/80 - end function a2 - function a3(x,y,z) - use psb_base_mod, only : psb_dpk_ - real(psb_dpk_) :: a3 - real(psb_dpk_), intent(in) :: x,y,z - a3=done/80 - end function a3 - function g(x,y,z) - use psb_base_mod, only : psb_dpk_, done, dzero - real(psb_dpk_) :: g - real(psb_dpk_), intent(in) :: x,y,z - g = dzero - if (x == done) then - g = done - else if (x == dzero) then - g = exp(y**2-z**2) - end if - end function g - - -end program psb_d_pde3d - -