New example programs

examples/fileread/mld_dexample_1lev.f90

@@ -193,12 +193,15 @@ program mld_dexample_1lev
 
   ! set RAS
 
-  call P%init('AS',info)
+!!$  call P%init('AS',info)
-
+!!$
-  ! set number of overlaps
+!!$  ! set number of overlaps
-
+!!$
-  call P%set('SUB_OVR',2,info)
+!!$  call P%set('SUB_OVR',2,info)
-    
+!!$    
+  call P%init('FBGS',info)
+  !call prec%set('smoother_sweeps', p_choice%outer_sweeps, info)
+      
   ! build the preconditioner
 
   t1 = psb_wtime()

examples/pdegen/mld_dexample_1lev.f90

@@ -93,7 +93,7 @@ program mld_dexample_1lev
   real(psb_dpk_) :: resmx, resmxp
   real(psb_dpk_) :: t1, t2, tprec
   character(len=5)   :: afmt='CSR'
-  character(len=20)  :: name
+  character(len=20)  :: name, kmethod
 
   ! initialize the parallel environment
   call psb_init(ictxt)
@@ -125,7 +125,7 @@ program mld_dexample_1lev
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
-  call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
+  call mld_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
        & a1,a2,a3,b1,b2,b3,c,g,info)  
   call psb_barrier(ictxt)
   t2 = psb_wtime() - t1
@@ -166,18 +166,19 @@ program mld_dexample_1lev
   call x%zero()
   call psb_geasb(x,desc_A,info)
 
-  ! solve Ax=b with preconditioned BiCGSTAB
+  ! solve Ax=b with preconditioned Krylov method: BiCGSTAB
+  kmethod = 'BiCGSTAB'
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
 
-  call psb_krylov('BICGSTAB',A,P,b,x,tol,desc_A,info,itmax,iter,err,itrace=1,istop=2)
+  call psb_krylov(kmethod,A,P,b,x,tol,desc_A,info,itmax,iter,err,itrace=1,istop=2)
 
   t2 = psb_wtime() - t1
   call psb_amx(ictxt,t2)
 
   call psb_geall(r,desc_A,info)
-  call r%set(dzero)
+  call r%zero()
   call psb_geasb(r,desc_A,info)
   call psb_geaxpby(done,b,dzero,r,desc_A,info)
   call psb_spmm(-done,A,x,done,r,desc_A,info)
@@ -197,6 +198,7 @@ program mld_dexample_1lev
     write(*,'(" ")')
     write(*,'("Matrix from PDE example")')
     write(*,'("Computed solution on ",i8," processors")')np
+    write(*,'("Krylov method             : ",a)') kmethod
     write(*,'("Iterations to convergence : ",i6)')iter
     write(*,'("Error estimate on exit    : ",es12.5)')err
     write(*,'("Time to build prec.       : ",es12.5)')tprec

examples/pdegen/mld_dexample_ml.f90

@@ -146,7 +146,7 @@ program mld_dexample_ml
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
-  call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
+  call mld_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
        & a1,a2,a3,b1,b2,b3,c,g,info)  
   call psb_barrier(ictxt)
   t2 = psb_wtime() - t1
@@ -219,7 +219,7 @@ program mld_dexample_ml
   call x%zero()
   call psb_geasb(x,desc_A,info)
 
-  ! solve Ax=b with preconditioned CG
+  ! solve Ax=b with preconditioned Krylov method
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
@@ -250,6 +250,7 @@ program mld_dexample_ml
     write(*,'(" ")')
     write(*,'("Matrix from PDE example")')
     write(*,'("Computed solution on ",i8," processors")')np
+    write(*,'("Krylov method             : ",a)') kmethod
     write(*,'("Iterations to convergence : ",i6)')iter
     write(*,'("Error estimate on exit    : ",es12.5)')err
     write(*,'("Time to build prec.       : ",es12.5)')tprec

examples/pdegen/mld_dpde_mod.f90

@@ -35,7 +35,355 @@
 !    POSSIBILITY OF SUCH DAMAGE.
 !   
 module mld_d_pde_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 mld_gen_pde3d
+    module procedure  mld_d_gen_pde3d
+  end interface mld_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 mld_d_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,&
+       & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,nrl,iv)
+    use psb_base_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 :: 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,glob_row,nlr,i,ii,ib,k
+    integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
+    integer(psb_ipk_) :: np, iam, nr, nt
+    integer(psb_ipk_) :: icoeff
+    integer(psb_ipk_), allocatable     :: irow(:),icol(:),myidx(:)
+    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   = 1.d0/(idim+2)
+    sqdeltah = deltah*deltah
+    deltah2  = 2.d0* deltah
+
+    ! 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
+
+    if (.not.present(iv)) then 
+      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
+    else
+      if (size(iv) /= m) then
+        write(psb_err_unit,*) iam, 'Initialization error IV',size(iv),m
+        info = -1
+        call psb_barrier(ictxt)
+        call psb_abort(ictxt)
+        return    
+      end if
+
+    end if
+    call psb_barrier(ictxt)
+    t0 = psb_wtime()
+    if (present(iv)) then 
+      call psb_cdall(ictxt,desc_a,info,vg=iv)
+    else
+      call psb_cdall(ictxt,desc_a,info,nl=nr)
+    end if
+
+    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
+
+    myidx = desc_a%get_global_indices()
+    nlr = size(myidx)
+
+    ! 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
+        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-1)*deltah
+        y = (iy-1)*deltah
+        z = (iz-1)*deltah
+        zt(k) = f_(x,y,z)
+        ! 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.d0*(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(:)=0.d0
+      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 mld_d_gen_pde3d
+
+
   !
   ! functions parametrizing the differential equation 
   !  

examples/pdegen/mld_sexample_1lev.f90

@@ -93,7 +93,7 @@ program mld_sexample_1lev
   real(psb_spk_) :: resmx, resmxp
   real(psb_dpk_) :: t1, t2, tprec
   character(len=5)   :: afmt='CSR'
-  character(len=20)  :: name
+  character(len=20)  :: name, kmethod
 
   ! initialize the parallel environment
   call psb_init(ictxt)
@@ -125,7 +125,7 @@ program mld_sexample_1lev
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
-  call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
+  call mld_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
        & a1,a2,a3,b1,b2,b3,c,g,info)  
   call psb_barrier(ictxt)
   t2 = psb_wtime() - t1
@@ -166,18 +166,19 @@ program mld_sexample_1lev
   call x%zero()
   call psb_geasb(x,desc_A,info)
 
-  ! solve Ax=b with preconditioned BiCGSTAB
+  ! solve Ax=b with preconditioned Krylov method: BiCGSTAB
+  kmethod = 'BiCGSTAB'
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
 
-  call psb_krylov('BICGSTAB',A,P,b,x,tol,desc_A,info,itmax,iter,err,itrace=1,istop=2)
+  call psb_krylov(kmethod,A,P,b,x,tol,desc_A,info,itmax,iter,err,itrace=1,istop=2)
 
   t2 = psb_wtime() - t1
   call psb_amx(ictxt,t2)
 
   call psb_geall(r,desc_A,info)
-  call r%set(szero)
+  call r%zero()
   call psb_geasb(r,desc_A,info)
   call psb_geaxpby(sone,b,szero,r,desc_A,info)
   call psb_spmm(-sone,A,x,sone,r,desc_A,info)
@@ -197,6 +198,7 @@ program mld_sexample_1lev
     write(*,'(" ")')
     write(*,'("Matrix from PDE example")')
     write(*,'("Computed solution on ",i8," processors")')np
+    write(*,'("Krylov method             : ",a)') kmethod
     write(*,'("Iterations to convergence : ",i6)')iter
     write(*,'("Error estimate on exit    : ",es12.5)')err
     write(*,'("Time to build prec.       : ",es12.5)')tprec

examples/pdegen/mld_sexample_ml.f90

@@ -146,7 +146,7 @@ program mld_sexample_ml
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
-  call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
+  call mld_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
        & a1,a2,a3,b1,b2,b3,c,g,info)  
   call psb_barrier(ictxt)
   t2 = psb_wtime() - t1
@@ -219,7 +219,7 @@ program mld_sexample_ml
   call x%zero()
   call psb_geasb(x,desc_A,info)
 
-  ! solve Ax=b with preconditioned CG
+  ! solve Ax=b with preconditioned Krylov method
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
@@ -250,6 +250,7 @@ program mld_sexample_ml
     write(*,'(" ")')
     write(*,'("Matrix from PDE example")')
     write(*,'("Computed solution on ",i8," processors")')np
+    write(*,'("Krylov method             : ",a)') kmethod
     write(*,'("Iterations to convergence : ",i6)')iter
     write(*,'("Error estimate on exit    : ",es12.5)')err
     write(*,'("Time to build prec.       : ",es12.5)')tprec

examples/pdegen/mld_spde_mod.f90

@@ -35,7 +35,355 @@
 !    POSSIBILITY OF SUCH DAMAGE.
 !   
 module mld_s_pde_mod
+  use psb_base_mod, only : psb_spk_, psb_ipk_, psb_desc_type,&
+       &  psb_sspmat_type, psb_s_vect_type, szero,&
+       &  psb_s_base_sparse_mat, psb_s_base_vect_type, psb_i_base_vect_type
+
+  interface 
+    function s_func_3d(x,y,z) result(val)
+      import :: psb_spk_
+      real(psb_spk_), intent(in) :: x,y,z
+      real(psb_spk_) :: val
+    end function s_func_3d
+  end interface 
+
+  interface mld_gen_pde3d
+    module procedure  mld_s_gen_pde3d
+  end interface mld_gen_pde3d
+
+  
 contains
+
+  function s_null_func_3d(x,y,z) result(val)
+
+    real(psb_spk_), intent(in) :: x,y,z
+    real(psb_spk_) :: val
+    
+    val = szero
+
+  end function s_null_func_3d
+
+  !
+  !  subroutine to allocate and fill in the coefficient matrix and
+  !  the rhs. 
+  !
+  subroutine mld_s_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,&
+       & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,nrl,iv)
+    use psb_base_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(s_func_3d)  :: b1,b2,b3,c,a1,a2,a3,g
+    integer(psb_ipk_)     :: idim
+    type(psb_sspmat_type) :: a
+    type(psb_s_vect_type) :: xv,bv
+    type(psb_desc_type)   :: desc_a
+    integer(psb_ipk_)     :: ictxt, info
+    character(len=*)      :: afmt
+    procedure(s_func_3d), optional :: f
+    class(psb_s_base_sparse_mat), optional :: amold
+    class(psb_s_base_vect_type), optional :: vmold 
+    class(psb_i_base_vect_type), optional :: imold
+    integer(psb_ipk_), optional :: nrl,iv(:)
+
+    ! Local variables.
+
+    integer(psb_ipk_), parameter :: nb=20
+    type(psb_s_csc_sparse_mat)  :: acsc
+    type(psb_s_coo_sparse_mat)  :: acoo
+    type(psb_s_csr_sparse_mat)  :: acsr
+    real(psb_spk_)           :: zt(nb),x,y,z
+    integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k
+    integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
+    integer(psb_ipk_) :: np, iam, nr, nt
+    integer(psb_ipk_) :: icoeff
+    integer(psb_ipk_), allocatable     :: irow(:),icol(:),myidx(:)
+    real(psb_spk_), allocatable :: val(:)
+    ! deltah dimension of each grid cell
+    ! deltat discretization time
+    real(psb_spk_)            :: deltah, sqdeltah, deltah2
+    real(psb_spk_), parameter :: rhs=szero,one=sone,zero=szero
+    real(psb_dpk_)    :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
+    integer(psb_ipk_) :: err_act
+    procedure(s_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_ => s_null_func_3d
+    end if
+
+    deltah   = 1.d0/(idim+2)
+    sqdeltah = deltah*deltah
+    deltah2  = 2.d0* deltah
+
+    ! 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
+
+    if (.not.present(iv)) then 
+      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
+    else
+      if (size(iv) /= m) then
+        write(psb_err_unit,*) iam, 'Initialization error IV',size(iv),m
+        info = -1
+        call psb_barrier(ictxt)
+        call psb_abort(ictxt)
+        return    
+      end if
+
+    end if
+    call psb_barrier(ictxt)
+    t0 = psb_wtime()
+    if (present(iv)) then 
+      call psb_cdall(ictxt,desc_a,info,vg=iv)
+    else
+      call psb_cdall(ictxt,desc_a,info,nl=nr)
+    end if
+
+    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
+
+    myidx = desc_a%get_global_indices()
+    nlr = size(myidx)
+
+    ! 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
+        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-1)*deltah
+        y = (iy-1)*deltah
+        z = (iz-1)*deltah
+        zt(k) = f_(x,y,z)
+        ! 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(szero,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,szero,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,szero)*(-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.d0*(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,sone)*(-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,sone,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(sone,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(:)=0.d0
+      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 mld_s_gen_pde3d
+
+
   !
   ! functions parametrizing the differential equation 
   !  

tests/fileread/runs/dfs.inp

@@ -1,11 +1,11 @@
-mld_mat.mtx                 ! This matrix (and others) from: http://math.nist.gov/MatrixMarket/ or 
+poisson_200x200.mtx         ! This matrix (and others) from: http://math.nist.gov/MatrixMarket/ or 
-mld_rhs.mtx                 ! rhs          | http://www.cise.ufl.edu/research/sparse/matrices/index.html
+NONE                        ! rhs          | http://www.cise.ufl.edu/research/sparse/matrices/index.html
 NONE                        ! Initial guess
-NONE                         !mld_sol.mtx                 ! Reference solution
+NONE                        !mld_sol.mtx                 ! Reference solution
 MM                          ! File format: MatrixMarket or Harwell-Boeing
 CSR                         ! Storage format: CSR COO JAD 
 GRAPH                       ! PART (partition method): BLOCK GRAPH
-BICGSTAB                   ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
+CG                          ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FCG GCR RGMRES
 2                           ! ISTOPC
 00500                       ! ITMAX
 1                           ! ITRACE
@@ -14,35 +14,35 @@ BICGSTAB                   ! Iterative method: BiCGSTAB BiCGSTABL BiCG CG CGS FC
 ML-VCYCLE-FBGS-R-UMF        ! Longer descriptive name for preconditioner (up to 20 chars) 
 ML                          ! Preconditioner type: NONE JACOBI BJAC AS ML
 1                           ! Number of sweeps for 1-level, outer sweeps for ML; lines below ignored for non-ML
-WCYCLE                      ! Type of multilevel CYCLE: VCYCLE WCYCLE KCYCLE MULT ADD
+VCYCLE                      ! Type of multilevel CYCLE: VCYCLE WCYCLE KCYCLE MULT ADD
 -3                          ! Max Number of levels in a multilevel preconditioner; if <0, lib default
--3                          ! Target coarse matrix size; if <0, lib default
+-500                         ! Target coarse matrix size; if <0, lib default
 SMOOTHED                    ! Type of aggregation: SMOOTHED UNSMOOTHED
-SYMDEC                         ! Parallel aggregation: DEC, SYMDEC
+DEC                         ! Parallel aggregation: DEC, SYMDEC
 NATURAL                     ! Ordering of aggregation NATURAL DEGREE
 NOFILTER                    ! Filtering of matrix ? FILTER NOFILTER
 -1.5                        ! Coarsening ratio, if < 0 use library default
-2                           ! Number of thresholds in vector, next line ignored if <= 0
+-2                          ! Number of thresholds in vector, next line ignored if <= 0
 0.05 0.025                  ! Thresholds 
-0.0100d0                    ! Smoothed aggregation threshold: >= 0.0
+0.0d0                       ! Smoothed aggregation threshold: >= 0.0
 BJAC                        ! Smoother type JACOBI GS BGS BJAC AS; ignored for non-ML
 1                           ! Number of sweeps for (pre) smoother
 0                           ! Number of overlap layers for AS preconditioner
 HALO                        ! AS restriction operator: NONE HALO 
 NONE                        ! AS prolongation operator: NONE SUM AVG
-ILU                      ! Subdomain solver: JACOBI BS BGS ILU ILUT MILU MUMPS SLU UMF
+GS                          ! Subdomain solver: JACOBI BS BGS ILU ILUT MILU MUMPS SLU UMF
 0                           ! Fill level P for ILU(P) and ILU(T,P)
 1.d-4                       ! Threshold T for ILU(T,P) 
-NONE                        ! Second (post) smoother, ignored if NONE
+BJAC                        ! Second (post) smoother, ignored if NONE
 1                           ! Number of sweeps for (post) smoother
 0                           ! Number of overlap layers for AS preconditioner
 HALO                        ! AS restriction operator: NONE HALO 
 NONE                        ! AS prolongation operator: NONE SUM AVG
-JACOBI                      ! Subdomain solver: JACOBI BS BGS ILU ILUT MILU MUMPS SLU UMF
+BWGS                        ! Subdomain solver: JACOBI BS BGS ILU ILUT MILU MUMPS SLU UMF
 0                           ! Fill level P for ILU(P) and ILU(T,P)
 1.d-4                       ! Threshold T for ILU(T,P) 
 DEFLT                       ! Coarsest-level solver: MUMPS UMF SLU SLUDIST JACOBI GS BJAC 
-ILU                        ! Coarsest-level subsolver: ILU UMF SLU MUMPS SLUDIST
+ILU                         ! Coarsest-level subsolver: ILU UMF SLU MUMPS SLUDIST
 DIST                        ! Coarsest-level matrix distribution: DIST  REPL           
 1                           ! Coarsest-level fillin P for ILU(P) and ILU(T,P)
 1.d-4                       ! Coarsest-level threshold T for ILU(T,P) 

tests/pdegen/mld_d_pde2d.f90

@@ -62,8 +62,314 @@
 ! data distribution.
 !
 module mld_d_pde2d_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_2d(x,y) result(val)
+      import :: psb_dpk_
+      real(psb_dpk_), intent(in) :: x,y
+      real(psb_dpk_) :: val
+    end function d_func_2d
+  end interface
+
+  interface mld_gen_pde2d 
+    module procedure mld_d_gen_pde2d
+  end interface mld_gen_pde2d
 contains
 
+  function d_null_func_2d(x,y) result(val)
+
+    real(psb_dpk_), intent(in) :: x,y
+    real(psb_dpk_) :: val
+
+    val = dzero
+
+  end function d_null_func_2d
+
+
+  !
+  !  subroutine to allocate and fill in the coefficient matrix and
+  !  the rhs. 
+  !
+  subroutine mld_d_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,&
+       & a1,a2,b1,b2,c,g,info,f,amold,vmold,imold,nrl)
+    use psb_base_mod
+    !
+    !   Discretizes the partial differential equation
+    ! 
+    !   a1 dd(u)  a2 dd(u)    b1 d(u)  b2 d(u) 
+    ! -   ------ -  ------ +  -----  +  ------  + c u = f
+    !      dxdx     dydy         dx       dy     
+    !
+    ! with Dirichlet boundary conditions
+    !   u = g 
+    !
+    !  on the unit square  0<=x,y<=1.
+    !
+    !
+    ! Note that if b1=b2=c=0., the PDE is the  Laplace equation.
+    !
+    implicit none
+    procedure(d_func_2d)  :: b1,b2,c,a1,a2,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_2d), 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 :: nrl
+
+    ! 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,glob_row,nlr,i,ii,ib,k
+    integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
+    integer(psb_ipk_) :: np, iam, nr, nt
+    integer(psb_ipk_) :: icoeff
+    integer(psb_ipk_), allocatable     :: irow(:),icol(:),myidx(:)
+    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=0.e0,one=1.e0,zero=0.e0
+    real(psb_dpk_)    :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
+    integer(psb_ipk_) :: err_act
+    procedure(d_func_2d), 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_2d
+    end if
+
+    deltah   = 1.e0/(idim+2)
+    sqdeltah = deltah*deltah
+    deltah2  = 2.e0* deltah
+
+    ! initialize array descriptor and sparse matrix storage. provide an
+    ! estimate of the number of non zeroes 
+
+    m   = idim*idim
+    n   = m
+    nnz = ((n*7)/(np))
+    if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
+
+    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
+    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(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
+
+
+    myidx = desc_a%get_global_indices()
+    nlr = size(myidx)
+
+    ! 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
+        if (mod(glob_row,(idim)) == 0) then
+          ix = glob_row/(idim)
+        else
+          ix = glob_row/(idim)+1
+        endif
+        iy = (glob_row-(ix-1)*idim)
+        ! x, y
+        x = (ix-1)*deltah
+        y = (iy-1)*deltah
+
+        zt(k) = f_(x,y)
+        ! internal point: build discretization
+        !   
+        !  term depending on   (x-1,y)
+        !
+        val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2
+        if (ix == 1) then 
+          zt(k) = g(dzero,y)*(-val(icoeff)) + zt(k)
+        else
+          icol(icoeff) = (ix-2)*idim+iy
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+        endif
+        !  term depending on     (x,y-1)
+        val(icoeff)  = -a2(x,y)/sqdeltah-b2(x,y)/deltah2
+        if (iy == 1) then 
+          zt(k) = g(x,dzero)*(-val(icoeff))   + zt(k)
+        else
+          icol(icoeff) = (ix-1)*idim+(iy-1)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+        endif
+
+        !  term depending on     (x,y)
+        val(icoeff)=2.e0*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y)
+        icol(icoeff) = (ix-1)*idim+iy
+        irow(icoeff) = glob_row
+        icoeff       = icoeff+1                  
+        !  term depending on     (x,y+1)
+        val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2
+        if (iy == idim) then 
+          zt(k) = g(x,done)*(-val(icoeff))   + zt(k)
+        else
+          icol(icoeff) = (ix-1)*idim+(iy+1)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+        endif
+        !  term depending on     (x+1,y)
+        val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2
+        if (ix==idim) then 
+          zt(k) = g(done,y)*(-val(icoeff))   + zt(k)
+        else
+          icol(icoeff) = (ix)*idim+(iy)
+          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(:)=0.e0
+      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 mld_d_gen_pde2d
+
+
   !
   ! functions parametrizing the differential equation 
   !  
@@ -249,13 +555,13 @@ program mld_d_pde2d
   !
   call psb_barrier(ictxt)
   t1 = psb_wtime()
-  call psb_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,&
+  call mld_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,&
        & a1,a2,b1,b2,c,g,info)  
   call psb_barrier(ictxt)
   t2 = psb_wtime() - t1
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
-    ch_err='psb_gen_pde2d'
+    ch_err='mld_gen_pde2d'
     call psb_errpush(info,name,a_err=ch_err)
     goto 9999
   end if

tests/pdegen/mld_d_pde3d.f90

@@ -63,7 +63,354 @@
 ! data distribution.
 !
 module mld_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 mld_gen_pde3d
+    module procedure  mld_d_gen_pde3d
+  end interface mld_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 mld_d_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,&
+       & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,nrl,iv)
+    use psb_base_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 :: 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,glob_row,nlr,i,ii,ib,k
+    integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
+    integer(psb_ipk_) :: np, iam, nr, nt
+    integer(psb_ipk_) :: icoeff
+    integer(psb_ipk_), allocatable     :: irow(:),icol(:),myidx(:)
+    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   = 1.d0/(idim+2)
+    sqdeltah = deltah*deltah
+    deltah2  = 2.d0* deltah
+
+    ! 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
+
+    if (.not.present(iv)) then 
+      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
+    else
+      if (size(iv) /= m) then
+        write(psb_err_unit,*) iam, 'Initialization error IV',size(iv),m
+        info = -1
+        call psb_barrier(ictxt)
+        call psb_abort(ictxt)
+        return    
+      end if
+
+    end if
+    call psb_barrier(ictxt)
+    t0 = psb_wtime()
+    if (present(iv)) then 
+      call psb_cdall(ictxt,desc_a,info,vg=iv)
+    else
+      call psb_cdall(ictxt,desc_a,info,nl=nr)
+    end if
+
+    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
+
+    myidx = desc_a%get_global_indices()
+    nlr = size(myidx)
+
+    ! 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
+        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-1)*deltah
+        y = (iy-1)*deltah
+        z = (iz-1)*deltah
+        zt(k) = f_(x,y,z)
+        ! 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.d0*(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(:)=0.d0
+      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 mld_d_gen_pde3d
+
   !
   ! functions parametrizing the differential equation 
   !  
@@ -262,13 +609,13 @@ program mld_d_pde3d
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
-  call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
+  call mld_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
        & a1,a2,a3,b1,b2,b3,c,g,info)  
   call psb_barrier(ictxt)
   t2 = psb_wtime() - t1
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
-    ch_err='psb_gen_pdeed'
+    ch_err='mld_gen_pde3d'
     call psb_errpush(info,name,a_err=ch_err)
     goto 9999
   end if

tests/pdegen/mld_s_pde2d.f90

@@ -62,8 +62,314 @@
 ! data distribution.
 !
 module mld_s_pde2d_mod
+  use psb_base_mod, only : psb_spk_, psb_ipk_, psb_desc_type,&
+       &  psb_sspmat_type, psb_s_vect_type, szero,&
+       &  psb_s_base_sparse_mat, psb_s_base_vect_type, psb_i_base_vect_type
+
+  interface 
+    function s_func_2d(x,y) result(val)
+      import :: psb_spk_
+      real(psb_spk_), intent(in) :: x,y
+      real(psb_spk_) :: val
+    end function s_func_2d
+  end interface
+
+  interface mld_gen_pde2d 
+    module procedure mld_s_gen_pde2d
+  end interface mld_gen_pde2d
 contains
 
+  function s_null_func_2d(x,y) result(val)
+
+    real(psb_spk_), intent(in) :: x,y
+    real(psb_spk_) :: val
+
+    val = szero
+
+  end function s_null_func_2d
+
+
+  !
+  !  subroutine to allocate and fill in the coefficient matrix and
+  !  the rhs. 
+  !
+  subroutine mld_s_gen_pde2d(ictxt,idim,a,bv,xv,desc_a,afmt,&
+       & a1,a2,b1,b2,c,g,info,f,amold,vmold,imold,nrl)
+    use psb_base_mod
+    !
+    !   Discretizes the partial differential equation
+    ! 
+    !   a1 dd(u)  a2 dd(u)    b1 d(u)  b2 d(u) 
+    ! -   ------ -  ------ +  -----  +  ------  + c u = f
+    !      dxdx     dydy         dx       dy     
+    !
+    ! with Dirichlet boundary conditions
+    !   u = g 
+    !
+    !  on the unit square  0<=x,y<=1.
+    !
+    !
+    ! Note that if b1=b2=c=0., the PDE is the  Laplace equation.
+    !
+    implicit none
+    procedure(s_func_2d)  :: b1,b2,c,a1,a2,g
+    integer(psb_ipk_)     :: idim
+    type(psb_sspmat_type) :: a
+    type(psb_s_vect_type) :: xv,bv
+    type(psb_desc_type)   :: desc_a
+    integer(psb_ipk_)     :: ictxt, info
+    character(len=*)      :: afmt
+    procedure(s_func_2d), optional :: f
+    class(psb_s_base_sparse_mat), optional :: amold
+    class(psb_s_base_vect_type), optional :: vmold
+    class(psb_i_base_vect_type), optional :: imold
+    integer(psb_ipk_), optional :: nrl
+
+    ! Local variables.
+
+    integer(psb_ipk_), parameter :: nb=20
+    type(psb_s_csc_sparse_mat)  :: acsc
+    type(psb_s_coo_sparse_mat)  :: acoo
+    type(psb_s_csr_sparse_mat)  :: acsr
+    real(psb_spk_)           :: zt(nb),x,y,z
+    integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k
+    integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
+    integer(psb_ipk_) :: np, iam, nr, nt
+    integer(psb_ipk_) :: icoeff
+    integer(psb_ipk_), allocatable     :: irow(:),icol(:),myidx(:)
+    real(psb_spk_), allocatable :: val(:)
+    ! deltah dimension of each grid cell
+    ! deltat discretization time
+    real(psb_spk_)            :: deltah, sqdeltah, deltah2
+    real(psb_spk_), parameter :: rhs=0.e0,one=1.e0,zero=0.e0
+    real(psb_dpk_)    :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
+    integer(psb_ipk_) :: err_act
+    procedure(s_func_2d), 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_ => s_null_func_2d
+    end if
+
+    deltah   = 1.e0/(idim+2)
+    sqdeltah = deltah*deltah
+    deltah2  = 2.e0* deltah
+
+    ! initialize array descriptor and sparse matrix storage. provide an
+    ! estimate of the number of non zeroes 
+
+    m   = idim*idim
+    n   = m
+    nnz = ((n*7)/(np))
+    if(iam == psb_root_) write(psb_out_unit,'("Generating Matrix (size=",i0,")...")')n
+
+    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
+    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(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
+
+
+    myidx = desc_a%get_global_indices()
+    nlr = size(myidx)
+
+    ! 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
+        if (mod(glob_row,(idim)) == 0) then
+          ix = glob_row/(idim)
+        else
+          ix = glob_row/(idim)+1
+        endif
+        iy = (glob_row-(ix-1)*idim)
+        ! x, y
+        x = (ix-1)*deltah
+        y = (iy-1)*deltah
+
+        zt(k) = f_(x,y)
+        ! internal point: build discretization
+        !   
+        !  term depending on   (x-1,y)
+        !
+        val(icoeff) = -a1(x,y)/sqdeltah-b1(x,y)/deltah2
+        if (ix == 1) then 
+          zt(k) = g(szero,y)*(-val(icoeff)) + zt(k)
+        else
+          icol(icoeff) = (ix-2)*idim+iy
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+        endif
+        !  term depending on     (x,y-1)
+        val(icoeff)  = -a2(x,y)/sqdeltah-b2(x,y)/deltah2
+        if (iy == 1) then 
+          zt(k) = g(x,szero)*(-val(icoeff))   + zt(k)
+        else
+          icol(icoeff) = (ix-1)*idim+(iy-1)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+        endif
+
+        !  term depending on     (x,y)
+        val(icoeff)=2.e0*(a1(x,y) + a2(x,y))/sqdeltah + c(x,y)
+        icol(icoeff) = (ix-1)*idim+iy
+        irow(icoeff) = glob_row
+        icoeff       = icoeff+1                  
+        !  term depending on     (x,y+1)
+        val(icoeff)=-a2(x,y)/sqdeltah+b2(x,y)/deltah2
+        if (iy == idim) then 
+          zt(k) = g(x,sone)*(-val(icoeff))   + zt(k)
+        else
+          icol(icoeff) = (ix-1)*idim+(iy+1)
+          irow(icoeff) = glob_row
+          icoeff       = icoeff+1
+        endif
+        !  term depending on     (x+1,y)
+        val(icoeff)=-a1(x,y)/sqdeltah+b1(x,y)/deltah2
+        if (ix==idim) then 
+          zt(k) = g(sone,y)*(-val(icoeff))   + zt(k)
+        else
+          icol(icoeff) = (ix)*idim+(iy)
+          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(:)=0.e0
+      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 mld_s_gen_pde2d
+
+
   !
   ! functions parametrizing the differential equation 
   !  
@@ -249,13 +555,13 @@ program mld_s_pde2d
   !
   call psb_barrier(ictxt)
   t1 = psb_wtime()
-  call psb_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,&
+  call mld_gen_pde2d(ictxt,idim,a,b,x,desc_a,afmt,&
        & a1,a2,b1,b2,c,g,info)  
   call psb_barrier(ictxt)
   t2 = psb_wtime() - t1
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
-    ch_err='psb_gen_pde2d'
+    ch_err='mld_gen_pde2d'
     call psb_errpush(info,name,a_err=ch_err)
     goto 9999
   end if

tests/pdegen/mld_s_pde3d.f90

@@ -63,7 +63,354 @@
 ! data distribution.
 !
 module mld_s_pde3d_mod
+  use psb_base_mod, only : psb_spk_, psb_ipk_, psb_desc_type,&
+       &  psb_sspmat_type, psb_s_vect_type, szero,&
+       &  psb_s_base_sparse_mat, psb_s_base_vect_type, psb_i_base_vect_type
+
+  interface 
+    function s_func_3d(x,y,z) result(val)
+      import :: psb_spk_
+      real(psb_spk_), intent(in) :: x,y,z
+      real(psb_spk_) :: val
+    end function s_func_3d
+  end interface 
+
+  interface mld_gen_pde3d
+    module procedure  mld_s_gen_pde3d
+  end interface mld_gen_pde3d
+
+  
 contains
+
+  function s_null_func_3d(x,y,z) result(val)
+
+    real(psb_spk_), intent(in) :: x,y,z
+    real(psb_spk_) :: val
+    
+    val = szero
+
+  end function s_null_func_3d
+
+  !
+  !  subroutine to allocate and fill in the coefficient matrix and
+  !  the rhs. 
+  !
+  subroutine mld_s_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,&
+       & a1,a2,a3,b1,b2,b3,c,g,info,f,amold,vmold,imold,nrl,iv)
+    use psb_base_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(s_func_3d)  :: b1,b2,b3,c,a1,a2,a3,g
+    integer(psb_ipk_)     :: idim
+    type(psb_sspmat_type) :: a
+    type(psb_s_vect_type) :: xv,bv
+    type(psb_desc_type)   :: desc_a
+    integer(psb_ipk_)     :: ictxt, info
+    character(len=*)      :: afmt
+    procedure(s_func_3d), optional :: f
+    class(psb_s_base_sparse_mat), optional :: amold
+    class(psb_s_base_vect_type), optional :: vmold 
+    class(psb_i_base_vect_type), optional :: imold
+    integer(psb_ipk_), optional :: nrl,iv(:)
+
+    ! Local variables.
+
+    integer(psb_ipk_), parameter :: nb=20
+    type(psb_s_csc_sparse_mat)  :: acsc
+    type(psb_s_coo_sparse_mat)  :: acoo
+    type(psb_s_csr_sparse_mat)  :: acsr
+    real(psb_spk_)           :: zt(nb),x,y,z
+    integer(psb_ipk_) :: m,n,nnz,glob_row,nlr,i,ii,ib,k
+    integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
+    integer(psb_ipk_) :: np, iam, nr, nt
+    integer(psb_ipk_) :: icoeff
+    integer(psb_ipk_), allocatable     :: irow(:),icol(:),myidx(:)
+    real(psb_spk_), allocatable :: val(:)
+    ! deltah dimension of each grid cell
+    ! deltat discretization time
+    real(psb_spk_)            :: deltah, sqdeltah, deltah2
+    real(psb_spk_), parameter :: rhs=szero,one=sone,zero=szero
+    real(psb_dpk_)    :: t0, t1, t2, t3, tasb, talc, ttot, tgen, tcdasb
+    integer(psb_ipk_) :: err_act
+    procedure(s_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_ => s_null_func_3d
+    end if
+
+    deltah   = 1.d0/(idim+2)
+    sqdeltah = deltah*deltah
+    deltah2  = 2.d0* deltah
+
+    ! 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
+
+    if (.not.present(iv)) then 
+      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
+    else
+      if (size(iv) /= m) then
+        write(psb_err_unit,*) iam, 'Initialization error IV',size(iv),m
+        info = -1
+        call psb_barrier(ictxt)
+        call psb_abort(ictxt)
+        return    
+      end if
+
+    end if
+    call psb_barrier(ictxt)
+    t0 = psb_wtime()
+    if (present(iv)) then 
+      call psb_cdall(ictxt,desc_a,info,vg=iv)
+    else
+      call psb_cdall(ictxt,desc_a,info,nl=nr)
+    end if
+
+    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
+
+    myidx = desc_a%get_global_indices()
+    nlr = size(myidx)
+
+    ! 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
+        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-1)*deltah
+        y = (iy-1)*deltah
+        z = (iz-1)*deltah
+        zt(k) = f_(x,y,z)
+        ! 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(szero,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,szero,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,szero)*(-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.d0*(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,sone)*(-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,sone,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(sone,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(:)=0.d0
+      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 mld_s_gen_pde3d
+
   !
   ! functions parametrizing the differential equation 
   !  
@@ -262,13 +609,13 @@ program mld_s_pde3d
 
   call psb_barrier(ictxt)
   t1 = psb_wtime()
-  call psb_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
+  call mld_gen_pde3d(ictxt,idim,a,b,x,desc_a,afmt,&
        & a1,a2,a3,b1,b2,b3,c,g,info)  
   call psb_barrier(ictxt)
   t2 = psb_wtime() - t1
   if(info /= psb_success_) then
     info=psb_err_from_subroutine_
-    ch_err='psb_gen_pdeed'
+    ch_err='mld_gen_pde3d'
     call psb_errpush(info,name,a_err=ch_err)
     goto 9999
   end if