Added Jacobi example and made corrections to expression compiler and diagonal iteration.

Renamed Accumulator->Transform, since Accumulator does not make so much sense for Endomorphism

Makefile

@@ -15,7 +15,7 @@ all: $(BDIR)/lib/libctf.a
 
 EXAMPLES = dft dft_3D gemm gemm_4D scalar trace weigh_4D subworld_gemm \
            permute_multiworld strassen slice_gemm ccsd sparse_permuted_slice qinformatics endomorphism endomorphism_cust endomorphism_cust_sp \
-           univar_function univar_accumulator_cust univar_accumulator_cust_sp spmv spmm
+           univar_function univar_accumulator_cust univar_accumulator_cust_sp spmv spmm jacobi
 
 TESTS = test_suite pgemm_test nonsq_pgemm_test diag_sym sym3 readwrite_test \
         ccsdt_t3_to_t2 ccsdt_map_test multi_tsr_sym diag_ctr readall_test  sy_times_ns repack \

examples/Makefile

@@ -2,7 +2,7 @@ EXAMPLES = dft dft_3D gemm gemm_4D scalar trace weigh_4D subworld_gemm \
            permute_multiworld strassen slice_gemm ccsd sparse_permuted_slice \
            qinformatics endomorphism endomorphism_cust endomorphism_cust_sp \
            univar_function univar_accumulator_cust univar_accumulator_cust_sp \
-           spmv spmm
+           spmv spmm jacobi
 
 .PHONY:
 $(EXAMPLES): %: $(BDIR)/bin/%

examples/endomorphism.cxx

@@ -31,7 +31,7 @@ int endomorphism(int     n,
 
   double scale = 1.0;
 
-  CTF::Accumulator<double> endo([=](double & d){ d=scale*d*d*d; });
+  CTF::Transform<double> endo([=](double & d){ d=scale*d*d*d; });
   // below is equivalent to A.scale(1.0, "ijkl", endo);
   endo(A["ijkl"]);
 

examples/endomorphism_cust.cxx

@@ -59,7 +59,7 @@ int endomorphism_cust(int     n,
  
   A.write(nvals, inds, vals);
 
-  CTF::Accumulator<cust_type> endo(
+  CTF::Transform<cust_type> endo(
     [](cust_type & a){
       a.len_name = strlen(a.name);
     });

examples/endomorphism_cust_sp.cxx

@@ -42,7 +42,7 @@ int endomorphism_cust_sp(int     n,
   } else
     A.write(0, NULL, NULL);
 
-  CTF::Accumulator<cust_sp_type> endo(comp_len);
+  CTF::Transform<cust_sp_type> endo(comp_len);
   // below is equivalent to A.scale(NULL, "ijkl", endo);
   endo(A["ijkl"]);
 

examples/jacobi.cxx

+/** \addtogroup examples 
+  * @{ 
+  * \defgroup jacobi jacobi
+  * @{ 
+  * \brief Multiplication of a random square sparse matrix by a vector
+  */
+
+#include <ctf.hpp>
+using namespace CTF;
+
+void jacobi_iter(Matrix<> & R, Vector<> & b, Vector<> & d, Vector<> &x){
+  x["i"] = -R["ij"]*x["j"];
+  x["i"] += b["i"];
+  x["i"] *= d["i"];
+}
+
+int jacobi(int     n,
+           World & dw){
+
+  Matrix<> spA(true,  n, n, NS, dw, "spA");
+  Matrix<> dnA(false, n, n, NS, dw, "dnA");
+  Vector<> b(n, dw);
+  Vector<> c1(n, dw);
+  Vector<> c2(n, dw);
+  Vector<> res(n, dw);
+
+  srand48(dw.rank);
+  b.fill_random(0.0,1.0);
+  c1.fill_random(0.0,1.0);
+  c2["i"] = c1["i"];
+
+  //make diagonally dominant matrix
+  dnA.fill_random(0.0,1.0);
+  spA["ij"] += dnA["ij"];
+  //sparsify
+  spA.sparsify(.5);
+  spA["ii"] += 2.*n;
+  dnA["ij"] = spA["ij"];
+
+  Vector<> d(n, dw);
+  d["i"] = spA["ii"];
+  Transform<> inv([](double & d){ d=1./d; });
+  inv(d["i"]);
+  
+  Matrix<> spR(true,  n, n, NS, dw, "spR");
+  Matrix<> dnR(false, n, n, NS, dw, "dnR");
+  spR["ij"] = spA["ij"];
+  dnR["ij"] = dnA["ij"];
+  spR["ii"] = 0;
+  dnR["ii"] = 0;
+
+/*  spR.print(); 
+  dnR.print(); */
+ 
+  //do up to 100 iterations
+  double res_norm;
+  int iter;
+  for (iter=0; iter<100; iter++){
+    jacobi_iter(dnR, b, d, c1);
+
+    res["i"]  = b["i"];
+    res["i"] -= dnA["ij"]*c1["j"];
+
+    res_norm = res.norm2();
+    if (res_norm < 1.E-4) break;
+  }
+#ifndef TEST_SUITE
+  if (dw.rank == 0)
+    printf("Completed %d iterations of Jacobi with dense matrix, residual F-norm is %E\n", iter, res_norm);
+#endif
+
+  for (iter=0; iter<100; iter++){
+    jacobi_iter(spR, b, d, c2);
+
+    res["i"]  = b["i"];
+    res["i"] -= spA["ij"]*c2["j"];
+
+    res_norm = res.norm2();
+    if (res_norm < 1.E-4) break;
+  }
+#ifndef TEST_SUITE
+  if (dw.rank == 0)
+    printf("Completed %d iterations of Jacobi with sparse matrix, residual F-norm is %E\n", iter, res_norm);
+#endif
+
+  c2["i"] -= c1["i"];
+
+  bool pass = c2.norm2() <= 1.E-6;
+
+  if (dw.rank == 0){
+    if (pass) 
+      printf("{ Jacobi x[\"i\"] = (1./A[\"ii\"])*(b[\"j\"] - (A[\"ij\"]-A[\"ii\"])*x[\"j\"]) with sparse A } passed \n");
+    else
+      printf("{ Jacobi x[\"i\"] = (1./A[\"ii\"])*(b[\"j\"] - (A[\"ij\"]-A[\"ii\"])*x[\"j\"]) with sparse A } failed \n");
+  }
+  return pass;
+} 
+
+
+#ifndef TEST_SUITE
+char* getCmdOption(char ** begin,
+                   char ** end,
+                   const   std::string & option){
+  char ** itr = std::find(begin, end, option);
+  if (itr != end && ++itr != end){
+    return *itr;
+  }
+  return 0;
+}
+
+
+int main(int argc, char ** argv){
+  int rank, np, n, pass;
+  int const in_num = argc;
+  char ** input_str = argv;
+
+  MPI_Init(&argc, &argv);
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  MPI_Comm_size(MPI_COMM_WORLD, &np);
+
+  if (getCmdOption(input_str, input_str+in_num, "-n")){
+    n = atoi(getCmdOption(input_str, input_str+in_num, "-n"));
+    if (n < 0) n = 7;
+  } else n = 7;
+
+
+  {
+    World dw(argc, argv);
+
+    if (rank == 0){
+      printf("Running Jacobi method on random %d-by-%d sparse matrix\n",n,n);
+    }
+    pass = jacobi(n, dw);
+    assert(pass);
+  }
+
+  MPI_Finalize();
+  return 0;
+}
+/**
+ * @} 
+ * @}
+ */
+
+#endif

examples/univar_accumulator_cust.cxx

@@ -48,7 +48,7 @@ int univar_accumulator_cust(int     n,
   free(loc_frcs);
   free(finds);
 
-  CTF::Accumulator<force,particle> uacc([](force f, particle & p){
+  CTF::Transform<force,particle> uacc([](force f, particle & p){
       p.dx += f.fx*p.coeff;
       p.dy += f.fy*p.coeff;
     });

examples/univar_accumulator_cust_sp.cxx

@@ -61,7 +61,7 @@ int univar_accumulator_cust_sp(int     n,
     
   F.write(my_forces.size(), &my_forces[0]);
 
-  CTF::Accumulator<force,particle> uacc(&acc_force);
+  CTF::Transform<force,particle> uacc(&acc_force);
 
   //FIXME = does not work because it sets beta to addid :/
   F2["ij"] += F["ij"];

src/interface/functions.h

@@ -99,7 +99,7 @@ namespace CTF {
    *          e.g. B["ij"] = f(A["ij"],B["ij"])
    */
   template<typename dtype_A=double, typename dtype_B=dtype_A>
-  class Univar_Accumulator : public CTF_int::univar_function {
+  class Univar_Transform : public CTF_int::univar_function {
     public:
       /**
        * \brief function signature for element-wise multiplication, compute b=f(a)
@@ -111,7 +111,7 @@ namespace CTF {
        * \brief constructor takes function pointers to compute B=f(A));
        * \param[in] f_ linear function (type_A)->(type_B)
        */
-      Univar_Accumulator(std::function<void(dtype_A, dtype_B &)> f_){ f = f_; }
+      Univar_Transform(std::function<void(dtype_A, dtype_B &)> f_){ f = f_; }
 
       /** 
        * \brief evaluate B=f(A) 
@@ -224,23 +224,23 @@ namespace CTF {
 
   
   template<typename dtype_A=double, typename dtype_B=dtype_A, typename dtype_C=dtype_A>
-  class Accumulator {
+  class Transform {
     public:
       bool is_endo;
       Endomorphism<dtype_A> * endo;
       bool is_univar;
-      Univar_Accumulator<dtype_A, dtype_B> * univar;
+      Univar_Transform<dtype_A, dtype_B> * univar;
 
-      Accumulator(std::function<void(dtype_A&)> f_){
+      Transform(std::function<void(dtype_A&)> f_){
         is_endo = true;
         is_univar = false;
         endo = new Endomorphism<dtype_A>(f_);
       }
       
-      Accumulator(std::function<void(dtype_A, dtype_B&)> f_){
+      Transform(std::function<void(dtype_A, dtype_B&)> f_){
         is_endo = false;
         is_univar = true;
-        univar = new Univar_Accumulator<dtype_A, dtype_B>(f_);
+        univar = new Univar_Transform<dtype_A, dtype_B>(f_);
       }
 
       void operator()(CTF_int::Term const & A) const {
@@ -253,7 +253,7 @@ namespace CTF {
         univar->operator()(A,B);
       }
       
-      operator Univar_Accumulator<dtype_A, dtype_B>(){
+      operator Univar_Transform<dtype_A, dtype_B>(){
         assert(is_univar);
         return *univar;
       }

src/interface/idx_tensor.cxx

@@ -206,6 +206,21 @@ namespace CTF {
       sr->safecopy(scale,sr->mulid());
     }
   }
+  
+  void Idx_Tensor::operator=(double scl){ execute() = Idx_Tensor(sr,scl); }
+  void Idx_Tensor::operator+=(double scl){ execute() += Idx_Tensor(sr,scl); }
+  void Idx_Tensor::operator-=(double scl){ execute() -= Idx_Tensor(sr,scl); }
+  void Idx_Tensor::operator*=(double scl){ execute() *= Idx_Tensor(sr,scl); }
+
+  void Idx_Tensor::operator=(int64_t scl){ execute() = Idx_Tensor(sr,scl); }
+  void Idx_Tensor::operator+=(int64_t scl){ execute() += Idx_Tensor(sr,scl); }
+  void Idx_Tensor::operator-=(int64_t scl){ execute() -= Idx_Tensor(sr,scl); }
+  void Idx_Tensor::operator*=(int64_t scl){ execute() *= Idx_Tensor(sr,scl); }
+
+  void Idx_Tensor::operator=(int scl){ execute() = Idx_Tensor(sr,(int64_t)scl); }
+  void Idx_Tensor::operator+=(int scl){ execute() += Idx_Tensor(sr,(int64_t)scl); }
+  void Idx_Tensor::operator-=(int scl){ execute() -= Idx_Tensor(sr,(int64_t)scl); }
+  void Idx_Tensor::operator*=(int scl){ execute() *= Idx_Tensor(sr,(int64_t)scl); }
 
   /*Idx_Tensor Idx_Tensor::operator-() const {
 
@@ -241,8 +256,6 @@ namespace CTF {
   }
 
 
-  void Idx_Tensor::operator=(double scl){ execute() = Idx_Tensor(sr,scl); }
-  void Idx_Tensor::operator=(int64_t scl){ execute() = Idx_Tensor(sr,scl); }
 
 
   void Idx_Tensor::execute(Idx_Tensor output) const {

src/interface/idx_tensor.h

@@ -96,7 +96,17 @@ namespace CTF {
 
       //same as in parent (Term) but not inherited in C++
       void operator=(double scl);
+      void operator+=(double scl);
+      void operator-=(double scl);
+      void operator*=(double scl);
       void operator=(int64_t scl);
+      void operator+=(int64_t scl);
+      void operator-=(int64_t scl);
+      void operator*=(int64_t scl);
+      void operator=(int scl);
+      void operator+=(int scl);
+      void operator-=(int scl);
+      void operator*=(int scl);
 
       /**
        * \brief A += B, compute any operations on operand B and add

src/interface/tensor.cxx

@@ -395,6 +395,12 @@ namespace CTF {
     int ret = CTF_int::tensor::sparsify((char*)&threshold, take_abs);
     assert(ret == CTF_int::SUCCESS);
   }
+  
+  template<typename dtype>
+  void Tensor<dtype>::sparsify(std::function<bool(dtype)> filter){
+    int ret = CTF_int::tensor::sparsify([&](char const * c){ filter(((dtype*)c)[0]); });
+    assert(ret == CTF_int::SUCCESS);
+  }
 
 
   template<typename dtype>

src/interface/tensor.h

@@ -651,7 +651,13 @@ namespace CTF {
       void sparsify(dtype threshold,
                     bool  take_abs=true);
 
-      
+     
+      /**
+       * \brief sparsifies tensor keeping only values v such that filter(v) = true
+       * \param[in] filter boolean function to apply to values to determine whether to keep them
+       */ 
+      void sparsify(std::function<bool(dtype)> filter);
+
      /**
        * \brief accumulates this tensor to a tensor object defined on a different world
        * \param[in] tsr a tensor object of the same characteristic as this tensor,

src/interface/term.cxx

@@ -127,7 +127,10 @@ namespace CTF_int {
   
   Term::~Term(){
     delete sr;
-    if (scale != NULL) free(scale);
+    if (scale != NULL){
+      free(scale);
+      scale = NULL;
+    }
   }
 
   Contract_Term Term::operator*(Term const & A) const {
@@ -166,6 +169,17 @@ namespace CTF_int {
   void Term::operator+=(int64_t scl){ execute() += Idx_Tensor(sr,scl); }
   void Term::operator-=(int64_t scl){ execute() -= Idx_Tensor(sr,scl); }
   void Term::operator*=(int64_t scl){ execute() *= Idx_Tensor(sr,scl); }
+  void Term::operator=(int scl){ execute() = Idx_Tensor(sr,(int64_t)scl); }
+  void Term::operator+=(int scl){ execute() += Idx_Tensor(sr,(int64_t)scl); }
+  void Term::operator-=(int scl){ execute() -= Idx_Tensor(sr,(int64_t)scl); }
+  void Term::operator*=(int scl){ execute() *= Idx_Tensor(sr,(int64_t)scl); }
+
+
+
+  Term & Term::operator-(){
+    sr->safeaddinv(scale,scale);
+    return *this;
+  }
 /*
   Contract_Term Contract_Term::operator-() const {
     Contract_Term trm(*this);
@@ -485,6 +499,7 @@ namespace CTF_int {
         c.execute();
       }
       if (tscale != NULL) free(tscale);
+      tscale = NULL;
       delete pop_A;
       delete pop_B;
     } 
@@ -524,7 +539,10 @@ namespace CTF_int {
       }
       delete pop_A;
       delete pop_B;
-      if (tscale != NULL) free(tscale);
+      if (tscale != NULL){
+        free(tscale);
+        tscale = NULL;
+      }
     } 
     Idx_Tensor rtsr = tmp_ops[0]->execute();
     delete tmp_ops[0];

src/interface/term.h

@@ -115,6 +115,8 @@ namespace CTF_int {
       Sum_Term operator-(double scl) const;
       Sum_Term operator-(int64_t scl) const;
       
+      Term & operator-();
+      
       /**
        * \brief A = B, compute any operations on operand B and set
        * \param[in] B tensor on the right hand side
@@ -134,6 +136,11 @@ namespace CTF_int {
       void operator+=(int64_t scl);
       void operator-=(int64_t scl);
       void operator*=(int64_t scl);
+
+      void operator=(int scl);
+      void operator+=(int scl);
+      void operator-=(int scl);
+      void operator*=(int scl);
       /**
        * \brief figures out what world this term lives on
        */

src/summation/spr_seq_sum.cxx

@@ -187,6 +187,7 @@ namespace CTF_int{
                char *&                 pprs_new,
                univar_function const * func,
                int64_t                 map_pfx){
+
     // determine how many unique keys there are in prs_tsr and prs_Write
     nnew = nB;
     for (int64_t t=0,ww=0; ww<nA*map_pfx; ww++){
@@ -218,6 +219,8 @@ namespace CTF_int{
       int64_t mw = ww%map_pfx;
       if (t<nB && (w==nA || prs_B[t].k() < prs_A[w].k()*map_pfx+mw)){
         memcpy(prs_new[n].ptr, prs_B[t].ptr, sr_B->pair_size());
+        if (beta != NULL)
+          sr_B->mul(prs_B[t].d(), beta, prs_new[n].d());
         t++;
       } else {
         if (t>=nB || prs_B[t].k() > prs_A[w].k()*map_pfx+mw){
@@ -292,13 +295,18 @@ namespace CTF_int{
                        int64_t                 size_A,
                        algstrct const *        sr_A,
                        char const *            beta,
-                       char const *            B,
+                       char *                  B,
                        int64_t                 size_B,
                        char *&                 new_B,
                        int64_t &               new_size_B,
                        algstrct const *        sr_B,
                        univar_function const * func,
                        int64_t                 map_pfx){
+
+/*      if (!sr_B->isequal(beta, sr_B->mulid())){
+        printf("scaling B by 0\n");
+        sr_B->scal(size_B, beta, B, 1);
+      }*/
       spspsum(sr_A, size_A, ConstPairIterator(sr_A, A), beta,
               sr_B, size_B, ConstPairIterator(sr_B, B),alpha,
               new_size_B, new_B, func, map_pfx);

src/summation/spr_seq_sum.h

@@ -87,7 +87,7 @@ namespace CTF_int {
                        int64_t                 size_A,
                        algstrct const *        sr_A,
                        char const *            beta,
-                       char const *            B,
+                       char *                  B,
                        int64_t                 size_B,
                        char *&                 new_B,
                        int64_t &               new_size_B,

src/tensor/untyped_tensor.cxx

@@ -1110,17 +1110,44 @@ namespace CTF_int {
 
   int tensor::sparsify(char const * threshold,
                        bool         take_abs){
+    if ((threshold == NULL && sr->addid() == NULL) ||
+        (threshold != NULL && !sr->is_ordered())){
+      return SUCCESS;
+    }
+    if (threshold == NULL)
+      return sparsify([&](char const* c){ return !sr->isequal(c, sr->addid()); });
+    else if (!take_abs)
+      return sparsify([&](char const* c){ 
+        char tmp[sr->el_size];
+        sr->max(c,threshold,tmp);
+        return !sr->isequal(tmp, threshold);
+      });
+    else
+      return sparsify([&](char const* c){ 
+        char tmp[sr->el_size];
+        sr->abs(c,tmp);
+        sr->max(tmp,threshold,tmp);
+        return !sr->isequal(tmp, threshold);
+      });
+  }
+
+  int tensor::sparsify(std::function<bool(char*)> f){
     if (is_sparse){
-      if ((threshold == NULL && sr->addid() == NULL) ||
-          (threshold != NULL && !sr->is_ordered())){
-        return SUCCESS;
-      }
       int64_t nnz_loc_new = 0;
       PairIterator pi(sr, data);
       int64_t nnz_blk_old[calc_nvirt()];
       memcpy(nnz_blk_old, nnz_blk, calc_nvirt()*sizeof(int64_t));
       memset(nnz_blk, 0, calc_nvirt()*sizeof(int64_t));
-      if (threshold == NULL){
+      int64_t i=0; 
+      for (int v=0; v<calc_nvirt(); v++){
+        for (int64_t j=0; j<nnz_blk_old[v]; j++,i++){
+          if (f(pi[i].d())){
+            nnz_loc_new++;
+            nnz_blk[v]++;
+          }
+        }
+      }
+/*      if (threshold == NULL){
         int64_t i=0; 
         for (int v=0; v<calc_nvirt(); v++){
           for (int64_t j=0; j<nnz_blk_old[v]; j++,i++){
@@ -1146,13 +1173,19 @@ namespace CTF_int {
             }
           }
         }
-      }
+      }*/
       // if we don't have any actual zeros don't do anything
       if (nnz_loc_new == nnz_loc) return SUCCESS;
       alloc_ptr(nnz_loc_new*sr->pair_size(), (void**)&data);
       PairIterator pi_new(sr, data);
       nnz_loc_new = 0;
-      if (threshold == NULL){
+      for (int64_t i=0; i<nnz_loc; i++){
+        if (f(pi[i].d())){
+          memcpy(pi_new[nnz_loc_new].ptr, pi[i].ptr, sr->pair_size());
+          nnz_loc_new++;
+        }
+      }
+      /*if (threshold == NULL){
         for (int64_t i=0; i<nnz_loc; i++){
           if (!sr->isequal(pi[i].d(), sr->addid())){
             memcpy(pi_new[nnz_loc_new].ptr, pi[i].ptr, sr->pair_size());
@@ -1172,7 +1205,7 @@ namespace CTF_int {
             nnz_loc_new++;
           }
         }
-      }
+      }*/
       nnz_loc = nnz_loc_new;
       //FIXME compute max nnz_loc?
     } else {
@@ -1196,7 +1229,7 @@ namespace CTF_int {
       write(num_pairs, sr->mulid(), sr->addid(), all_pairs);
       cdealloc(all_pairs);
       //sparsify sparse->sparse
-      sparsify();
+      sparsify(f);
     }
     return SUCCESS;
   }
@@ -1999,11 +2032,11 @@ namespace CTF_int {
           }
           if (is_sparse){
             int64_t lda_i=1, lda_j=1;
-            for (int ii=1; ii<i; ii++){
-              lda_i *= lens[ii-1];
+            for (int ii=0; ii<i; ii++){
+              lda_i *= lens[ii];
             }
-            for (int jj=1; jj<j; jj++){
-              lda_j *= lens[jj-1];
+            for (int jj=0; jj<j; jj++){
+              lda_j *= lens[jj];
             }
             if (rw){
               PairIterator pi(sr, data);
@@ -2045,6 +2078,14 @@ namespace CTF_int {
                 ((int64_t*)(wdata[p].ptr))[0] = k_new;
 //                printf("k = %ld, k_new = %ld lda_i = %ld lda_j = %ld lens[0] = %d lens[1] = %d\n", k,k_new,lda_i,lda_j,lens[0],lens[1]);
               }
+              PairIterator pi(sr, this->data);
+              for (int p=0; p<nnz_loc; p++){
+                int64_t k = pi[p].k();
+                if ((k/lda_i)%lens[i] == (k/lda_j)%lens[j]){
+                  pi[p].write_val(sr->addid());
+                }
+              }
+
               this->write(nw, sr->mulid(), sr->addid(), pwdata);
               cdealloc(pwdata);
             }

src/tensor/untyped_tensor.h

@@ -415,6 +415,12 @@ namespace CTF_int {
        */
       int sparsify(char const * threshold=NULL,
                    bool         take_abs=true);
+  
+      /**
+       * \brief sparsifies tensor keeping only values v such that filter(v) = true
+       * \param[in] f boolean function to apply to values to determine whether to keep them
+       */ 
+      int sparsify(std::function<bool(char*)> f);
 
       /**
        * \brief read tensor data pairs local to processor including those with zero values

test/test_suite.cxx

@@ -34,6 +34,7 @@
 #include "../examples/univar_accumulator_cust_sp.cxx"
 #include "../examples/spmv.cxx"
 #include "../examples/spmm.cxx"
+#include "../examples/jacobi.cxx"
 
 using namespace CTF;
 
@@ -242,6 +243,10 @@ int main(int argc, char ** argv){
     if (rank == 0)
       printf("Testing sparse-matrix times matrix with n=%d k=%d:\n",n*n,n);
     pass.push_back(spmm(n*n,n,dw));
+    
+    if (rank == 0)
+      printf("Testing Jacobi iteration with n=%d:\n",n);
+    pass.push_back(jacobi(n,dw));
   }
   int num_pass = std::accumulate(pass.begin(), pass.end(), 0);
   if (rank == 0)