From fd820c363f65cf7ce2876764c85dc470f64a61a7 Mon Sep 17 00:00:00 2001
From: Salvatore Filippone <salvatore.filippone@cranfield.ac.uk>
Date: Tue, 1 May 2018 14:01:46 +0100
Subject: [PATCH] Fixed use of C_BASE in init_vl, and C test program.

---
 cbind/base/psb_base_tools_cbind_mod.F90 |  11 +-
 cbind/test/pargen/ppdec.c               | 171 ++++++++++++------------
 test/pargen/psb_d_pde2d.f90             |   4 +-
 test/pargen/psb_d_pde3d.f90             |   1 -
 test/pargen/psb_s_pde2d.f90             |   4 +-
 test/pargen/psb_s_pde3d.f90             |   1 -
 6 files changed, 100 insertions(+), 92 deletions(-)

diff --git a/cbind/base/psb_base_tools_cbind_mod.F90 b/cbind/base/psb_base_tools_cbind_mod.F90
index ea40197a..b2a4fd1d 100644
--- a/cbind/base/psb_base_tools_cbind_mod.F90
+++ b/cbind/base/psb_base_tools_cbind_mod.F90
@@ -2,6 +2,7 @@ module psb_base_tools_cbind_mod
   use iso_c_binding
   use psb_base_mod
   use psb_objhandle_mod
+  use psb_cpenv_mod
   use psb_base_string_cbind_mod
     
 contains
@@ -62,7 +63,7 @@ contains
     integer(psb_c_lpk)        :: vl(*)
     type(psb_c_object_type) :: cdh
     type(psb_desc_type), pointer :: descp
-    integer               :: info
+    integer               :: info, ixb 
 
     res = -1
     if (nl <=0) then 
@@ -79,8 +80,14 @@ contains
 
     allocate(descp,stat=info)
     if (info < 0) return 
+
+    ixb = psb_c_get_index_base()
       
-    call psb_cdall(ictxt,descp,info,vl=vl(1:nl))
+    if (ixb == 1) then 
+      call psb_cdall(ictxt,descp,info,vl=vl(1:nl))
+    else
+      call psb_cdall(ictxt,descp,info,vl=(vl(1:nl)+(1-ixb)))
+    end if
     cdh%item = c_loc(descp)
     res = info
 
diff --git a/cbind/test/pargen/ppdec.c b/cbind/test/pargen/ppdec.c
index f4ebd6ed..88e9d944 100644
--- a/cbind/test/pargen/ppdec.c
+++ b/cbind/test/pargen/ppdec.c
@@ -98,15 +98,15 @@ double  c(double x, double y, double  z)
 }
 double  b1(double x, double y, double  z)
 {
-  return(1.0/sqrt(3.0));
+  return(0.0/sqrt(3.0));
 }
 double b2(double x, double y, double  z)
 {
-  return(1.0/sqrt(3.0));
+  return(0.0/sqrt(3.0));
 }
 double b3(double x, double y, double  z)
 {
-  return(1.0/sqrt(3.0));
+  return(0.0/sqrt(3.0));
 }
 
 double g(double x, double y, double z)
@@ -120,13 +120,13 @@ double g(double x, double y, double z)
   }
 }
 
-psb_i_t matgen(psb_i_t ictxt, psb_l_t ng, psb_i_t idim, psb_i_t vg[],
+psb_i_t matgen(psb_i_t ictxt, psb_i_t nl, psb_i_t idim, psb_l_t vl[],
 	       psb_c_dspmat *ah,psb_c_descriptor *cdh,
 	       psb_c_dvector *xh, psb_c_dvector *bh, psb_c_dvector *rh)
 {
   psb_i_t iam, np;
   psb_l_t ix, iy, iz, el,glob_row;
-  psb_i_t i, info,ret;
+  psb_i_t i, k, info,ret;
   double x, y, z, deltah, sqdeltah, deltah2;
   double val[10*NBMAX], zt[NBMAX];
   psb_l_t irow[10*NBMAX], icol[10*NBMAX];
@@ -137,79 +137,77 @@ psb_i_t matgen(psb_i_t ictxt, psb_l_t ng, psb_i_t idim, psb_i_t vg[],
   sqdeltah = deltah*deltah;
   deltah2  = 2.0* deltah;
   psb_c_set_index_base(0);
-  for (glob_row=0; glob_row < ng; glob_row++) {
-
-    /* Check if I have to do something about this entry */
-    if (vg[glob_row] == iam) {
-      el=0;
-      ix = glob_row/(idim*idim);
-      iy = (glob_row-ix*idim*idim)/idim;
-      iz = glob_row-ix*idim*idim-iy*idim;
-      x=(ix+1)*deltah;
-      y=(iy+1)*deltah;
-      z=(iz+1)*deltah;
-      zt[0] = 0.0;
-      /*  internal point: build discretization */
-      /*  term depending on   (x-1,y,z)        */
-      val[el] = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2;	
-      if (ix==0) {
-	zt[0] += g(0.0,y,z)*(-val[el]);
-      } else {
-	icol[el]=(ix-1)*idim*idim+(iy)*idim+(iz);
-	el=el+1;
-      }
-      /*  term depending on     (x,y-1,z) */
-      val[el]  = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2;	      
-      if (iy==0) { 
-	zt[0] += g(x,0.0,z)*(-val[el]);
-      } else {
-	icol[el]=(ix)*idim*idim+(iy-1)*idim+(iz);
-	el=el+1;
-      }
-      /* term depending on     (x,y,z-1)*/
-      val[el]=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2;
-      if (iz==0) { 
-	zt[0] += g(x,y,0.0)*(-val[el]);
-      } else {
-	icol[el]=(ix)*idim*idim+(iy)*idim+(iz-1);
-	el=el+1;
-      }
-      /* term depending on     (x,y,z)*/
-      val[el]=2.0*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah + c(x,y,z);
-      icol[el]=(ix)*idim*idim+(iy)*idim+(iz);
+  for (i=0; i<nl;  i++) {
+    glob_row=vl[i];
+    //if ((i%100000 == 0)||(i<10)) fprintf(stderr,"%d: generation loop at %d %ld \n",iam,i,glob_row);
+    el=0;
+    ix = glob_row/(idim*idim);
+    iy = (glob_row-ix*idim*idim)/idim;
+    iz = glob_row-ix*idim*idim-iy*idim;
+    x=(ix+1)*deltah;
+    y=(iy+1)*deltah;
+    z=(iz+1)*deltah;
+    zt[0] = 0.0;
+    /*  internal point: build discretization */
+    /*  term depending on   (x-1,y,z)        */
+    val[el] = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2;	
+    if (ix==0) {
+      zt[0] += g(0.0,y,z)*(-val[el]);
+    } else {
+      icol[el]=(ix-1)*idim*idim+(iy)*idim+(iz);
+      el=el+1;
+    }
+    /*  term depending on     (x,y-1,z) */
+    val[el]  = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2;	      
+    if (iy==0) { 
+      zt[0] += g(x,0.0,z)*(-val[el]);
+    } else {
+      icol[el]=(ix)*idim*idim+(iy-1)*idim+(iz);
+      el=el+1;
+    }
+    /* term depending on     (x,y,z-1)*/
+    val[el]=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2;
+    if (iz==0) { 
+      zt[0] += g(x,y,0.0)*(-val[el]);
+    } else {
+      icol[el]=(ix)*idim*idim+(iy)*idim+(iz-1);
+      el=el+1;
+    }
+    /* term depending on     (x,y,z)*/
+    val[el]=2.0*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah + c(x,y,z);
+    icol[el]=(ix)*idim*idim+(iy)*idim+(iz);
+    el=el+1;
+    /*  term depending on     (x,y,z+1) */
+    val[el] = -a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2;
+    if (iz==idim-1) { 
+      zt[0] += g(x,y,1.0)*(-val[el]);
+    } else {
+      icol[el]=(ix)*idim*idim+(iy)*idim+(iz+1);
       el=el+1;
-      /*  term depending on     (x,y,z+1) */
-      val[el] = -a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2;
-      if (iz==idim-1) { 
-	zt[0] += g(x,y,1.0)*(-val[el]);
-      } else {
-	icol[el]=(ix)*idim*idim+(iy)*idim+(iz+1);
-	el=el+1;
-      }
-      /* term depending on     (x,y+1,z) */
-      val[el] = -a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2;
-      if (iy==idim-1) { 
-	zt[0] += g(x,1.0,z)*(-val[el]);
-      } else {
-	icol[el]=(ix)*idim*idim+(iy+1)*idim+(iz);
-	el=el+1;
-      }
-      /*  term depending on     (x+1,y,z) */
-      val[el] = -a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2;
-      if (ix==idim-1) { 
-	zt[0] += g(1.0,y,z)*(-val[el]);
-      } else {
-	icol[el]=(ix+1)*idim*idim+(iy)*idim+(iz);
-	el=el+1;
-      }
-      for (i=0; i<el; i++) irow[i]=glob_row;
-      if ((ret=psb_c_dspins(el,irow,icol,val,ah,cdh))!=0) 
-	fprintf(stderr,"From psb_c_dspins: %d\n",ret);
-      irow[0] = glob_row; 
-      psb_c_dgeins(1,irow,zt,bh,cdh);
-      zt[0]=0.0;
-      psb_c_dgeins(1,irow,zt,xh,cdh);
     }
+    /* term depending on     (x,y+1,z) */
+    val[el] = -a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2;
+    if (iy==idim-1) { 
+      zt[0] += g(x,1.0,z)*(-val[el]);
+    } else {
+      icol[el]=(ix)*idim*idim+(iy+1)*idim+(iz);
+      el=el+1;
+    }
+    /*  term depending on     (x+1,y,z) */
+    val[el] = -a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2;
+    if (ix==idim-1) { 
+      zt[0] += g(1.0,y,z)*(-val[el]);
+    } else {
+      icol[el]=(ix+1)*idim*idim+(iy)*idim+(iz);
+      el=el+1;
+    }
+    for (k=0; k<el; k++) irow[k]=glob_row;
+    if ((ret=psb_c_dspins(el,irow,icol,val,ah,cdh))!=0) 
+      fprintf(stderr,"From psb_c_dspins: %d\n",ret);
+    irow[0] = glob_row; 
+    psb_c_dgeins(1,irow,zt,bh,cdh);
+    zt[0]=0.0;
+    psb_c_dgeins(1,irow,zt,xh,cdh);   
   }
   
   if ((info=psb_c_cdasb(cdh))!=0)  return(info);
@@ -232,8 +230,8 @@ int main(int argc, char *argv[])
   psb_c_dprec *ph;
   psb_c_dspmat *ah;
   psb_c_dvector *bh, *xh, *rh; 
-  psb_i_t *vg, nb,nlr;
-  psb_l_t i,ng;
+  psb_i_t nb,nlr, nl;
+  psb_l_t i,ng, *vl, k;
   double t1,t2,eps,err;
   double *xv, *bv, *rv;
   double one=1.0, zero=0.0, res2;
@@ -284,18 +282,23 @@ int main(int argc, char *argv[])
   psb_c_barrier(ictxt);
 
   cdh=psb_c_new_descriptor();
+  psb_c_set_index_base(0);
 
   /* Simple minded BLOCK data distribution */ 
   ng = ((psb_l_t) idim)*idim*idim;  
   nb = (ng+np-1)/np;
-  if ((vg=malloc(ng*sizeof(psb_i_t)))==NULL) {
+  nl = nb;
+  if ( (ng -iam*nb) < nl) nl = ng -iam*nb; 
+  fprintf(stderr,"%d: Input data %d %ld %d %d\n",iam,idim,ng,nb, nl);
+  if ((vl=malloc(nb*sizeof(psb_l_t)))==NULL) {
     fprintf(stderr,"On %d: malloc failure\n",iam);
     psb_c_abort(ictxt);
   }
-  for (i=0; i<ng; i++) {
-    vg[i] = i/nb; 
-  }
-  if ((info=psb_c_cdall_vg(ng,vg,ictxt,cdh))!=0) {
+  i = ((psb_l_t)iam) * nb;
+  for (k=0; k<nl; k++)
+    vl[k] = i+k; 
+
+  if ((info=psb_c_cdall_vl(nl,vl,ictxt,cdh))!=0) {
     fprintf(stderr,"From cdall: %d\nBailing out\n",info);
     psb_c_abort(ictxt);
   }
@@ -316,7 +319,7 @@ int main(int argc, char *argv[])
 
   
   /* Matrix generation  */
-  if (matgen(ictxt,ng,idim,vg,ah,cdh,xh,bh,rh) != 0) {
+  if (matgen(ictxt,nl,idim,vl,ah,cdh,xh,bh,rh) != 0) {
     fprintf(stderr,"Error during matrix build loop\n");
     psb_c_abort(ictxt);
   }    
diff --git a/test/pargen/psb_d_pde2d.f90 b/test/pargen/psb_d_pde2d.f90
index bb7ccb08..00fd7061 100644
--- a/test/pargen/psb_d_pde2d.f90
+++ b/test/pargen/psb_d_pde2d.f90
@@ -184,7 +184,7 @@ contains
     integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
     integer(psb_lpk_) :: m,n,glob_row,nt
     integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
-    ! For 3D partition
+    ! For 2D partition
     ! Note: integer control variables going directly into an MPI call
     ! must be 4 bytes, i.e. psb_mpk_
     integer(psb_mpk_) :: npdims(2), npp, minfo
@@ -448,7 +448,7 @@ contains
       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
+      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
diff --git a/test/pargen/psb_d_pde3d.f90 b/test/pargen/psb_d_pde3d.f90
index 2651a436..a6fbad78 100644
--- a/test/pargen/psb_d_pde3d.f90
+++ b/test/pargen/psb_d_pde3d.f90
@@ -612,7 +612,6 @@ program psb_d_pde3d
   if(psb_get_errstatus() /= 0) goto 9999
   name='pde3d90'
   call psb_set_errverbosity(itwo)
-  call psb_cd_set_large_threshold(itwo)
   !
   ! Hello world
   !
diff --git a/test/pargen/psb_s_pde2d.f90 b/test/pargen/psb_s_pde2d.f90
index 40f47afb..be487e17 100644
--- a/test/pargen/psb_s_pde2d.f90
+++ b/test/pargen/psb_s_pde2d.f90
@@ -184,7 +184,7 @@ contains
     integer(psb_ipk_) :: nnz,nr,nlr,i,j,ii,ib,k, partition_
     integer(psb_lpk_) :: m,n,glob_row,nt
     integer(psb_ipk_) :: ix,iy,iz,ia,indx_owner
-    ! For 3D partition
+    ! For 2D partition
     ! Note: integer control variables going directly into an MPI call
     ! must be 4 bytes, i.e. psb_mpk_
     integer(psb_mpk_) :: npdims(2), npp, minfo
@@ -448,7 +448,7 @@ contains
       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
+      zt(:)=szero
       call psb_geins(ib,myidx(ii:ii+ib-1),zt(1:ib),xv,desc_a,info)
       if(info /= psb_success_) exit
     end do
diff --git a/test/pargen/psb_s_pde3d.f90 b/test/pargen/psb_s_pde3d.f90
index e3a5ed21..19d880ba 100644
--- a/test/pargen/psb_s_pde3d.f90
+++ b/test/pargen/psb_s_pde3d.f90
@@ -612,7 +612,6 @@ program psb_s_pde3d
   if(psb_get_errstatus() /= 0) goto 9999
   name='pde3d90'
   call psb_set_errverbosity(itwo)
-  call psb_cd_set_large_threshold(itwo)
   !
   ! Hello world
   !