Implemented a bunch of useless shit because I am scared of writing my qual proposal:

added matrix and vector types
added 3D FFT test

examples/Makefile

 include ../config.mk
 include ../rules.mk
 
-examples: fft gemm gemm_4D
+examples: fft fft_3D gemm gemm_4D
 fft: ${bindir}/test/fft
 ${bindir}/test/fft: fft.o
 
+fft_3D: ${bindir}/test/fft_3D
+${bindir}/test/fft_3D: fft_3D.o
+
 gemm_4D: ${bindir}/test/gemm_4D
 ${bindir}/test/gemm_4D: gemm_4D.o
 

examples/fft.cxx

@@ -23,12 +23,10 @@ int main(int argc, char ** argv){
     logn = 5;
   }
   n = 1<<logn;
-  int edge_lens[2] = {n,n};
-  int sym[2]       = {SY,NS};
 
   tCTF_World< std::complex<double> > * wrld = new tCTF_World< std::complex<double> >();
-  tCTF_Tensor< std::complex<double> > DFT(2, edge_lens, sym, wrld);
-  tCTF_Tensor< std::complex<double> > IDFT(2, edge_lens, sym, wrld);
+  tCTF_Matrix< std::complex<double> > DFT(n, n, SY, wrld);
+  tCTF_Matrix< std::complex<double> > IDFT(n, n, SY, wrld);
 
   DFT.get_local_data(&np, &idx, &data);
 

examples/fft_3D.cxx

+#include <ctf.hpp>
+#include <assert.h>
+#include <stdlib.h>
+
+
+
+
+/**
+ * \brief Forms N-by-N DFT matrix A and inverse-dft iA and checks A*iA=I
+ */
+int main(int argc, char ** argv){
+  int myRank, numPes, logn, i, j;
+  int64_t  n, np;
+  int64_t * idx;
+  std::complex<double> * data;
+  std::complex<double> imag(0,1);
+  
+  MPI_Init(&argc, &argv);
+  MPI_Comm_size(MPI_COMM_WORLD, &numPes);
+  MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
+
+  if (argc > 1){
+    logn = atoi(argv[1]);
+    if (logn<0) logn = 5;
+  } else {
+    logn = 5;
+  }
+  n = 1<<logn;
+
+  tCTF_World< std::complex<double> > * wrld = new tCTF_World< std::complex<double> >();
+  tCTF_Matrix< std::complex<double> > DFT(n, n, SY, wrld);
+  tCTF_Matrix< std::complex<double> > IDFT(n, n, SY, wrld);
+  tCTF_Tensor< std::complex<double> > MESH(3, (int[]){n, n, n}, (int[]){NS, NS, NS}, wrld);
+
+  DFT.get_local_data(&np, &idx, &data);
+
+  for (i=0; i<np; i++){
+    data[i] = (1./n)*exp(-2.*(idx[i]/n)*(idx[i]%n)*(M_PI/n)*imag);
+  }
+  DFT.write_remote_data(np, idx, data);
+  //DFT.print(stdout);
+  free(idx);
+  free(data); 
+  
+  IDFT.get_local_data(&np, &idx, &data);
+
+  for (i=0; i<np; i++){
+    data[i] = (1./n)*exp(2.*(idx[i]/n)*(idx[i]%n)*(M_PI/n)*imag);
+  }
+  IDFT.write_remote_data(np, idx, data);
+  //IDFT.print(stdout);
+  free(idx);
+  free(data); 
+
+  MESH.get_local_data(&np, &idx, &data);
+  for (i=0; i<np; i++){
+    for (j=0; j<n; j++){
+      data[i] += exp(imag*((-2.*M_PI*(j/(double)(n)))
+                      *((idx[i]%n) + ((idx[i]/n)%n) +(idx[i]/(n*n)))));
+    }
+  }
+  MESH.write_remote_data(np, idx, data);
+  //MESH.print(stdout);
+  free(idx);
+  free(data); 
+  
+  MESH["ijk"] = MESH["pqr"]*DFT["ip"]*DFT["jq"]*DFT["kr"];
+ 
+  MESH.get_local_data(&np, &idx, &data);
+  //MESH.print(stdout);
+  for (i=0; i<np; i++){
+    if (idx[i]%n == (idx[i]/n)%n && idx[i]%n == idx[i]/(n*n))
+      assert(fabs(data[i].real() - 1.)<=1.E-9);
+    else 
+      assert(fabs(data[i].real())<=1.E-9);
+  }
+  
+  if (myRank == 0)
+    printf("{ 3D_IDFT(3D_DFT(I))) = I } confirmed\n");
+
+  MPI_Barrier(MPI_COMM_WORLD);
+
+  free(idx);
+  free(data);
+  
+}

examples/gemm.cxx

@@ -29,16 +29,10 @@ void gemm(int const  m,
     printf("m = %d, n = %d, k = %d, p = %d, niter = %d\n", 
             m,n,k,num_pes,niter);
   
-  int shape_NS[] = {NS,NS};
-  int shape[] = {sym,NS};
-  int size_A[] = {m,k};
-  int size_B[] = {k,n};
-  int size_C[] = {m,n};
-
   //* Creates distributed tensors initialized with zeros
-  CTF_Tensor A = CTF_Tensor(2, size_A, shape, dw);
-  CTF_Tensor B = CTF_Tensor(2, size_B, shape, dw);
-  CTF_Tensor C = CTF_Tensor(2, size_C, shape_NS, dw);
+  CTF_Matrix A(m, k, sym, dw);
+  CTF_Matrix B(k, n, sym, dw);
+  CTF_Matrix C(m, n, NS,  dw);
 
   if (rank == 0)
     printf("tensor creation succeed\n");
@@ -76,8 +70,8 @@ void gemm(int const  m,
     if (rank == 0)
       printf("Verifying associativity\n");
     /* verify D=(A*B)*C = A*(B*C) */
-    CTF_Tensor D = CTF_Tensor(2, size_C, shape_NS, dw);
-    CTF_Tensor E = CTF_Tensor(2, size_C, shape_NS, dw);
+    CTF_Matrix D(m, n, NS, dw);
+    CTF_Matrix E(m, n, NS, dw);
     
     D["ij"] = A["ik"]*B["kj"];
     D["ij"] = D["ik"]*C["kj"];

examples/gemm_4D.cxx

@@ -31,9 +31,9 @@ void gemm_4D(int const  n,
   int sizeN4[] = {n,n,n,n};
 
   //* Creates distributed tensors initialized with zeros
-  CTF_Tensor A = CTF_Tensor(4, sizeN4, shapeN4, dw);
-  CTF_Tensor B = CTF_Tensor(4, sizeN4, shapeN4, dw);
-  CTF_Tensor C = CTF_Tensor(4, sizeN4, shapeN4, dw);
+  CTF_Tensor A(4, sizeN4, shapeN4, dw);
+  CTF_Tensor B(4, sizeN4, shapeN4, dw);
+  CTF_Tensor C(4, sizeN4, shapeN4, dw);
 
   if (rank == 0)
     printf("tensor creation succeed\n");
@@ -73,7 +73,7 @@ void gemm_4D(int const  n,
   }
   
   /* verify D=(A*B)*C = A*(B*C) */
-  CTF_Tensor D = CTF_Tensor(4, sizeN4, shapeN4, dw);
+  CTF_Tensor D(4, sizeN4, shapeN4, dw);
   
   D["ijkl"] = A["ijmn"]*B["mnkl"];
   D["ijkl"] = D["ijmn"]*C["mnkl"];
@@ -129,8 +129,7 @@ int main(int argc, char ** argv){
   } else niter = 3;
 
 
-  CTF_World * dw;
-  dw = new CTF_World();
+  CTF_World * dw = new CTF_World();
 
   if (rank == 0){
     printf("Computing C_ijkl = A_ijmn*B_klmn\n");
@@ -146,7 +145,6 @@ int main(int argc, char ** argv){
   }
   gemm_4D(n, AS, niter, dw, dir);
 
-  delete dw;
 
   MPI_Finalize();
   return 0;

include/ctf.hpp

@@ -52,6 +52,11 @@ class tCTF_Tensor {
     tCTF_World<dtype> * world;
 
   public:
+    /**
+     * \breif default constructor sets nothing 
+     */
+    tCTF_Tensor(){};
+
     /**
      * \brief copies a tensor (setting data to zero or copying A)
      * \param[in] A tensor to copy
@@ -196,6 +201,45 @@ class tCTF_Tensor {
     ~tCTF_Tensor();
 };
 
+/**
+ * \brief Matrix class which encapsulates a 2D tensor 
+ */
+template<typename dtype> 
+class tCTF_Matrix : public tCTF_Tensor<dtype> {
+  public:
+    int nrow, ncol, sym;
+
+    /**
+     * \brief constructor for a matrix
+     * \param[in] nrow number of matrix rows
+     * \param[in] ncol number of matrix columns
+     * \param[in] sym symmetry of matrix
+     * \param[in] world CTF world where the tensor will live
+     */ 
+    tCTF_Matrix(int const           nrow_, 
+                int const           ncol_, 
+                int const           sym_,
+                tCTF_World<dtype> * world);
+
+};
+
+/**
+ * \brief Vector class which encapsulates a 1D tensor 
+ */
+template<typename dtype> 
+class tCTF_Vector : public tCTF_Tensor<dtype> {
+  public:
+    int len;
+
+    /**
+     * \brief constructor for a vector
+     * \param[in] len_ dimension of vector
+     * \param[in] world CTF world where the tensor will live
+     */ 
+    tCTF_Vector(int const          len_,
+               tCTF_World<dtype> * world);
+};
+
 template<typename dtype> static
 tCTF_Idx_Tensor<dtype>& operator*(double d, tCTF_Idx_Tensor<dtype>& tsr){
   return tsr*d;
@@ -285,5 +329,7 @@ class tCTF_Idx_Tensor {
 
 typedef tCTF<double> CTF;
 typedef tCTF_Tensor<double> CTF_Tensor;
+typedef tCTF_Matrix<double> CTF_Matrix;
+typedef tCTF_Vector<double> CTF_Vector;
 typedef tCTF_World<double> CTF_World;
 #endif

rules.mk

 all: $(DEFAULT_COMPONENTS)
 
-ALL_COMPONENTS = ctf test_model pgemm_test nonsq_pgemm_test bench bench_model nonsq_pgemm_bench examples fft gemm gemm_4D
+ALL_COMPONENTS = src ctf test_model pgemm_test nonsq_pgemm_test bench bench_model nonsq_pgemm_bench examples fft fft_3D gemm gemm_4D
 
-examples: fft gemm
-test_model bench_model pgemm_test nonsq_pgemm_test nonsq_pgemm_bench fft gemm gemm_4D: ctf 
+examples: fft fft_3D gemm
+test_model bench_model pgemm_test nonsq_pgemm_test nonsq_pgemm_bench fft fft_3D gemm gemm_4D: ctf 
 
 bindir = ${top_dir}/bin
 libdir = ${top_dir}/lib

src/Makefile

@@ -21,6 +21,8 @@ ctf_SUBDIRS = ctr_seq ctr_comm shared dist_tensor interface
 ctf: ${libdir}/libctf.a
 ${libdir}/libctf.a: interface/ctf_world.o \
                     interface/ctf_tensor.o \
+                    interface/ctf_matrix.o \
+                    interface/ctf_vector.o \
                     interface/ctf_idx_tensor.o \
                     shared/comm.o \
                     shared/util.o \
@@ -30,5 +32,5 @@ ${libdir}/libctf.a: interface/ctf_world.o \
                     dist_tensor/dist_tensor.o \
                     unit_test/unit_test.o 
 
-INCLUDES += -I${top_dir}/src/ctr_comm -I${top_dir}/src/ctr_seq -I${top_dir}/dist_tensor
+INCLUDES += -I${top_dir}/src/ctr_comm -I${top_dir}/src/ctr_seq -I${top_dir}/src/dist_tensor -I${top_dir}/src/util -I${top_dir}/src/interface
 

src/interface/ctf_idx_tensor.cxx

@@ -39,6 +39,7 @@ tCTF_Idx_Tensor<dtype> * get_intermediate(tCTF_Idx_Tensor<dtype>* A,
   idx_C = (char*)malloc(sizeof(char)*ndim_C);
   sym_C = (int*)malloc(sizeof(int)*ndim_C);
   len_C = (int*)malloc(sizeof(int)*ndim_C);
+  idx = 0;
   for (i=0; i<A->parent->ndim; i++){
     for (j=0; j<B->parent->ndim; j++){
       if (A->idx_map[i] == B->idx_map[j]){
@@ -81,8 +82,8 @@ tCTF_Idx_Tensor<dtype> * get_intermediate(tCTF_Idx_Tensor<dtype>* A,
 
 template<typename dtype>
 tCTF_Idx_Tensor<dtype>::tCTF_Idx_Tensor(tCTF_Tensor<dtype> * parent_, const char * idx_map_){
-  idx_map = (char*)malloc((strlen(idx_map_)+1)*sizeof(char));
-  strcpy(idx_map, idx_map_);
+  idx_map = (char*)malloc(parent_->ndim*sizeof(char));
+  memcpy(idx_map, idx_map_,parent_->ndim*sizeof(char));
   parent        = parent_;
   has_contract  = 0;
   has_scale     = 0;
@@ -94,7 +95,6 @@ tCTF_Idx_Tensor<dtype>::tCTF_Idx_Tensor(tCTF_Tensor<dtype> * parent_, const char
 template<typename dtype>
 tCTF_Idx_Tensor<dtype>::~tCTF_Idx_Tensor(){
   free(idx_map);
-  if (is_intm) free(parent);
 }
 
 template<typename dtype>
@@ -116,8 +116,10 @@ void tCTF_Idx_Tensor<dtype>::operator*=(tCTF_Idx_Tensor<dtype>& tsr){
 
 template<typename dtype>
 tCTF_Idx_Tensor<dtype>& tCTF_Idx_Tensor<dtype>::operator* (tCTF_Idx_Tensor<dtype>& tsr){
-  if (has_contract || has_sum)
-    return (*NBR)*tsr;
+  if (has_contract || has_sum){
+    (*NBR)*tsr;
+    return *this;
+  }
   NBR = &tsr;
   has_contract = 1;
   return *this;

src/interface/ctf_matrix.cxx

+#include <ctf.hpp>
+
+template<typename dtype>
+tCTF_Matrix<dtype>::tCTF_Matrix(int const           nrow_,
+                                int const           ncol_,
+                                int const           sym_,
+                                tCTF_World<dtype> * world) :
+  tCTF_Tensor<dtype>(2, (int[]){nrow_, ncol_}, (int[]){sym_, NS}, world) {
+  nrow = nrow_;
+  ncol = ncol_;
+  sym = sym_;
+  
+}
+
+template class tCTF_Matrix<double>;
+template class tCTF_Matrix< std::complex<double> >;

src/interface/ctf_vector.cxx

+#include <ctf.hpp>
+
+template<typename dtype>
+tCTF_Vector<dtype>::tCTF_Vector(int const           len_,
+                                tCTF_World<dtype> * world) :
+  tCTF_Tensor<dtype>(2, &len, (int[]){NS}, world) {
+  len = len_;
+  
+}
+
+template class tCTF_Vector<double>;
+template class tCTF_Vector< std::complex<double> >;