diff --git a/src/interface/multilinear.cxx b/src/interface/multilinear.cxx index 0ccfd58d..680e394a 100644 --- a/src/interface/multilinear.cxx +++ b/src/interface/multilinear.cxx @@ -339,7 +339,7 @@ namespace CTF { } Matrix A(nrow_U, ncol_VT, SP*dA.is_sparse, *dA.wrld, *dA.sr); if (need_transpose_A){ - Tensor T(dA.order, dA.is_sparse, dA.lens, *dA.wrld, *dA.sr); + Tensor T(dA.order, dA.is_sparse, unf_lens_A, *dA.wrld, *dA.sr); T[unf_idx_A] += dA.operator[](idx_A); A.reshape(T); } else { diff --git a/src/interface/tensor.cxx b/src/interface/tensor.cxx index 73863f88..1ddf6022 100644 --- a/src/interface/tensor.cxx +++ b/src/interface/tensor.cxx @@ -885,7 +885,16 @@ NORM1_INST(double) for (int i=0; i([](dtype a){ return (double)std::norm(a); })(A[inds])); + //nrm = std::sqrt((double)Function([](dtype a){ return (double)std::norm(a); })(A[inds])); + Tensor cA(A.order, A.is_sparse, A.lens, *A.wrld, *A.sr); + cA[inds] += A[inds]; + Transform([](dtype & a){ a = std::abs(a)*std::abs(a); })(cA[inds]); + Tensor sc(0, NULL, *A.wrld); + sc[""] = cA[inds]; + dtype val = ((dtype*)sc.data)[0]; + MPI_Bcast((char *)&val, sizeof(dtype), MPI_CHAR, 0, A.wrld->comm); + nrm = std::abs(val); + } diff --git a/src/summation/summation.cxx b/src/summation/summation.cxx index 97d57f1d..bd1ccf2d 100644 --- a/src/summation/summation.cxx +++ b/src/summation/summation.cxx @@ -125,8 +125,10 @@ namespace CTF_int { #endif //update_all_models(A->wrld->cdt.cm); int stat = home_sum_tsr(run_diag); - if (stat != SUCCESS) + if (stat != SUCCESS){ printf("CTF ERROR: Failed to perform summation\n"); + ASSERT(0); + } } double summation::estimate_time(){ diff --git a/src_python/ctf/core.pyx b/src_python/ctf/core.pyx index 7b5165cd..2b1fef7d 100644 --- a/src_python/ctf/core.pyx +++ b/src_python/ctf/core.pyx @@ -68,8 +68,10 @@ cdef extern from "ctf.hpp" namespace "CTF_int": cdef cppclass ctensor "CTF_int::tensor": World * wrld algstrct * sr + int * lens bool is_sparse int64_t nnz_tot + int order ctensor() ctensor(ctensor * other, bool copy, bool alloc_data) ctensor(ctensor * other, int * new_sym) @@ -190,7 +192,7 @@ cdef extern from "../ctf_ext.h" namespace "CTF_int": cdef void matrix_svd_rand(ctensor * A, ctensor * U, ctensor * S, ctensor * VT, int rank, int iter, int oversmap, ctensor * U_init); cdef void matrix_svd_rand_cmplx(ctensor * A, ctensor * U, ctensor * S, ctensor * VT, int rank, int iter, int oversmap, ctensor * U_init); cdef void tensor_svd(ctensor * dA, char * idx_A, char * idx_U, char * idx_VT, int rank, double threshold, bool use_svd_rand, int num_iter, int oversamp, ctensor ** USVT) - #cdef void tensor_svd_cmplx(Tensor * dA, char * idx_A, char * idx_U, char * idx_VT, int rank, double threshold, bool use_svd_rand, int num_iter, int oversamp, ctensor ** USVT) + cdef void tensor_svd_cmplx(ctensor * dA, char * idx_A, char * idx_U, char * idx_VT, int rank, double threshold, bool use_svd_rand, int num_iter, int oversamp, ctensor ** USVT) cdef void conv_type(int type_idx1, int type_idx2, ctensor * A, ctensor * B) cdef void delete_arr(ctensor * A, char * arr) cdef void delete_pairs(ctensor * A, char * pairs) @@ -505,19 +507,18 @@ cdef class itensor(term): self.it.multeq(s) return self - - def svd(self, row_string, col_string, rank=None, threshold=None, use_svd_rand=False, num_iter=1, oversamp=5): + def svd(self, U_string, VT_string, rank=None, threshold=None, use_svd_rand=False, num_iter=1, oversamp=5): """ - svd(self, row_string, col_string, rank=None, threshold=None, use_svd_rand=False, num_iter=1, oversamp=5) + svd(self, U_string, VT_string, rank=None, threshold=None, use_svd_rand=False, num_iter=1, oversamp=5) Compute Single Value Decomposition of a given transposition/matricization of tensor A. Parameters ---------- - row_string: char string - Indices indexing row of matricization, should be subset of the string of this itensor + U_string: char string + Indices indexing left singular vectors, should be subset of the string of this itensor plus an auxiliary index - col_string: char string - Indices indexing columns of matricization, should be subset of the string of this itensor + VT_string: char string + Indices indexing right singular vectors, should be subset of the string of this itensor, plus same auxiliary index as in U threshold: real double precision or None, optional threshold for truncation of singular values. Either rank or threshold must be set to None. @@ -531,13 +532,13 @@ cdef class itensor(term): Returns ------- U: tensor - A unitary CTF tensor with dimensions len(row_string)+1. + A unitary CTF tensor with dimensions len(U_string). S: tensor A 1-D tensor with singular values. VT: tensor - A unitary CTF tensor with dimensions len(col_string)+1. + A unitary CTF tensor with dimensions len(VT_string)+1. """ t_svd = timer("pyTSVD") t_svd.start() @@ -546,23 +547,37 @@ cdef class itensor(term): if threshold == None: threshold = 0. cdef ctensor ** ctsrs = malloc(sizeof(ctensor*)*3) - if self.tsr.order == 'F': - row_string = _rev_array(row_string) - col_string = _rev_array(col_string) + if _ord_comp(self.tsr.order, 'F'): + U_string = _rev_array(U_string) + VT_string = _rev_array(VT_string) if self.tsr.dtype == np.float64 or self.tsr.dtype == np.float32: - tensor_svd(self.tsr.dt, self.string.encode(), col_string.encode(), row_string.encode(), rank, threshold, use_svd_rand, num_iter, oversamp, ctsrs) -# elif A.dtype == np.complex128 or A.dtype == np.complex64: + tensor_svd(self.tsr.dt, self.string.encode(), VT_string.encode(), U_string.encode(), rank, threshold, use_svd_rand, num_iter, oversamp, ctsrs) + elif self.tsr.dtype == np.complex128 or self.tsr.dtype == np.complex64: + tensor_svd_cmplx(self.tsr.dt, self.string.encode(), VT_string.encode(), U_string.encode(), rank, threshold, use_svd_rand, num_iter, oversamp, ctsrs) else: raise ValueError('CTF PYTHON ERROR: SVD must be called on real or complex single/double precision tensor') - U = tensor([],dtype=self.tsr.dtype,order=self.tsr.order) - S = tensor([],dtype=self.tsr.dtype,order=self.tsr.order) - VT = tensor([],dtype=self.tsr.dtype,order=self.tsr.order) + cdef cnp.ndarray lens_U = cnp.ndarray(ctsrs[2].order,dtype=np.int) + cdef cnp.ndarray lens_S = cnp.ndarray(ctsrs[1].order,dtype=np.int) + cdef cnp.ndarray lens_VT = cnp.ndarray(ctsrs[0].order,dtype=np.int) + for i in range(ctsrs[0].order): + lens_VT[i] = ctsrs[0].lens[i] + if _ord_comp(self.tsr.order, 'F'): + lens_VT = _rev_array(lens_VT) + for i in range(ctsrs[1].order): + lens_S[i] = ctsrs[1].lens[i] + for i in range(ctsrs[2].order): + lens_U[i] = ctsrs[2].lens[i] + if _ord_comp(self.tsr.order, 'F'): + lens_U = _rev_array(lens_U) + U = tensor(lens_U,dtype=self.tsr.dtype,order=self.tsr.order) + S = tensor(lens_S,dtype=self.tsr.dtype,order=self.tsr.order) + VT = tensor(lens_VT,dtype=self.tsr.dtype,order=self.tsr.order) del U.dt del S.dt del VT.dt - U.dt = ctsrs[0] + U.dt = ctsrs[2] S.dt = ctsrs[1] - VT.dt = ctsrs[2] + VT.dt = ctsrs[0] free(ctsrs) t_svd.stop() return [U, S, VT] @@ -1380,6 +1395,10 @@ cdef class tensor: (self.dt).fill_random(mn,mx) elif self.dtype == np.float64: (self.dt).fill_random(mn,mx) + elif self.dtype == np.complex64: + (self.dt).fill_random(mn,mx) + elif self.dtype == np.complex128: + (self.dt).fill_random(mn,mx) else: raise ValueError('CTF PYTHON ERROR: bad dtype') @@ -2269,7 +2288,7 @@ cdef class tensor: delete_arr(self.dt, cdata) return inds, vals - def _tot_size(self, unpack=True): + def tot_size(self, unpack=True): return self.dt.get_tot_size(not unpack) def read_all(self, arr=None, unpack=True): @@ -2783,10 +2802,10 @@ cdef class tensor: Returns ------- - output: scalar + output: scalar np.float64 2-norm of the tensor. - Examples + xamples -------- >>> import ctf >>> a = ctf.ones([3,4], dtype=np.float64) @@ -2794,19 +2813,21 @@ cdef class tensor: 3.4641016151377544 """ if self.dtype == np.float64: - return (self.dt).norm2() + return np.float64((self.dt).norm2()) elif self.dtype == np.float32: - return (self.dt).norm2() + return np.float64((self.dt).norm2()) elif self.dtype == np.int64: - return (self.dt).norm2() + return np.float64((self.dt).norm2()) elif self.dtype == np.int32: - return (self.dt).norm2() + return np.float64((self.dt).norm2()) elif self.dtype == np.int16: - return (self.dt).norm2() + return np.float64((self.dt).norm2()) elif self.dtype == np.int8: - return (self.dt).norm2() -# elif self.dtype == np.complex128: -# return (self.dt).norm2() + return np.float64((self.dt).norm2()) + elif self.dtype == np.complex64: + return np.float64(np.abs(np.complex64((self.dt).norm2()))) + elif self.dtype == np.complex128: + return np.float64(np.abs(np.complex128((self.dt).norm2()))) else: raise ValueError('CTF PYTHON ERROR: norm not present for this dtype') @@ -2874,7 +2895,7 @@ cdef class tensor: [1., 1., 1., 1.], [1., 1., 1., 1.]]) """ - vals = np.zeros(self._tot_size(), dtype=self.dtype) + vals = np.zeros(self.tot_size(), dtype=self.dtype) self.read_all(vals) #return np.asarray(np.ascontiguousarray(np.reshape(vals, self.shape, order='F')),order='C') #return np.reshape(vals, _rev_array(self.shape)).transpose() @@ -2909,9 +2930,9 @@ cdef class tensor: if self.dt.wrld.np == 1: self.write_all(arr) elif self.dt.wrld.rank == 0: - #self.write(np.arange(0,self._tot_size(),dtype=np.int64),np.asfortranarray(arr).flatten()) - #self.write(np.arange(0,self._tot_size(),dtype=np.int64),np.asfortranarray(arr).flatten()) - self.write(np.arange(0,self._tot_size(),dtype=np.int64),arr.ravel()) + #self.write(np.arange(0,self.tot_size(),dtype=np.int64),np.asfortranarray(arr).flatten()) + #self.write(np.arange(0,self.tot_size(),dtype=np.int64),np.asfortranarray(arr).flatten()) + self.write(np.arange(0,self.tot_size(),dtype=np.int64),arr.ravel()) else: self.write([], []) @@ -5101,7 +5122,7 @@ def _comp_all(tensor A, axis=None, out=None, keepdims=None): raise ValueError("'keepdims' not supported for all yet") if axis is None: x = A._bool_sum() - return x == A._tot_size() + return x == A.tot_size() def transpose(init_A, axes=None): """ diff --git a/src_python/ctf_ext.cxx b/src_python/ctf_ext.cxx index df5a262a..3199ecfb 100644 --- a/src_python/ctf_ext.cxx +++ b/src_python/ctf_ext.cxx @@ -507,6 +507,48 @@ namespace CTF_int{ } + void tensor_svd_cmplx(tensor * dA, char * idx_A, char * idx_U, char * idx_VT, int rank, double threshold, bool use_svd_rand, int num_iter, int oversamp, tensor ** USVT){ + char idx_S[2]; + idx_S[1] = '\0'; + for (int i=0; i<(int)strlen(idx_U); i++){ + for (int j=0; j<(int)strlen(idx_VT); j++){ + if (idx_U[i] == idx_VT[j]) idx_S[0] = idx_U[i]; + } + } + switch (dA->sr->el_size){ + case 8: + { + CTF::Tensor> * U = new CTF::Tensor>(); + CTF::Tensor> * S = new CTF::Tensor>(); + CTF::Tensor> * VT = new CTF::Tensor>(); + ((CTF::Tensor>*)dA)->operator[](idx_A).svd(U->operator[](idx_U), S->operator[](idx_S), VT->operator[](idx_VT), rank, threshold, use_svd_rand, num_iter, oversamp); + USVT[0] = U; + USVT[1] = S; + USVT[2] = VT; + } + break; + + + case 16: + { + CTF::Tensor> * U = new CTF::Tensor>(); + CTF::Tensor> * S = new CTF::Tensor>(); + CTF::Tensor> * VT = new CTF::Tensor>(); + ((CTF::Tensor>*)dA)->operator[](idx_A).svd(U->operator[](idx_U), S->operator[](idx_S), VT->operator[](idx_VT), rank, threshold, use_svd_rand, num_iter, oversamp); + USVT[0] = U; + USVT[1] = S; + USVT[2] = VT; + //printf("A dims %d %d, U dims %d %d, S dim %d, mVT dms %d %d)\n",mA.nrow, mA.ncol, mU.nrow, mU.ncol, vS.len, mVT.nrow, mVT.ncol); + } + break; + + default: + printf("CTF ERROR: SVD called on invalid tensor element type\n"); + assert(0); + break; + } + } + /* template <> void tensor::conv_type>(tensor * B){ diff --git a/src_python/ctf_ext.h b/src_python/ctf_ext.h index 246e07ef..4bb02694 100644 --- a/src_python/ctf_ext.h +++ b/src_python/ctf_ext.h @@ -96,6 +96,7 @@ namespace CTF_int{ void matrix_svd_rand(tensor * A, tensor * U, tensor * S, tensor * VT, int rank, int iter, int oversamp, tensor * U_init); void matrix_svd_rand_cmplx(tensor * A, tensor * U, tensor * S, tensor * VT, int rank, int iter, int oversamp, tensor * U_init); void tensor_svd(tensor * dA, char * idx_A, char * idx_U, char * idx_VT, int rank, double threshold, bool use_svd_rand, int num_iter, int oversamp, tensor ** USVT); + void tensor_svd_cmplx(tensor * dA, char * idx_A, char * idx_U, char * idx_VT, int rank, double threshold, bool use_svd_rand, int num_iter, int oversamp, tensor ** USVT); /** * \brief convert tensor from one type to another diff --git a/test/python/test_la.py b/test/python/test_la.py index 4db9ecb7..ada82530 100644 --- a/test/python/test_la.py +++ b/test/python/test_la.py @@ -75,7 +75,7 @@ def test_svd(self): m = 9 n = 5 k = 5 - for dt in [numpy.float32, numpy.float64]: + for dt in [numpy.float32, numpy.float64, numpy.complex64, numpy.complex128]: A = ctf.random.random((m,n)) A = ctf.astensor(A,dtype=dt) [U,S,VT]=ctf.svd(A,k) @@ -84,38 +84,11 @@ def test_svd(self): self.assertTrue(allclose(ctf.eye(k), ctf.dot(U.T(), U))) self.assertTrue(allclose(ctf.eye(k), ctf.dot(VT, VT.T()))) - A = ctf.tensor((m,n),dtype=numpy.complex64) - rA = ctf.tensor((m,n),dtype=numpy.float32) - rA.fill_random() - A.real(rA) - iA = ctf.tensor((m,n),dtype=numpy.float32) - iA.fill_random() - A.imag(iA) - - [U,S,VT]=ctf.svd(A,k) - self.assertTrue(allclose(A, ctf.dot(U,ctf.dot(ctf.diag(S),VT)))) - - self.assertTrue(allclose(ctf.eye(k,dtype=numpy.complex64), ctf.dot(ctf.conj(U.T()), U))) - self.assertTrue(allclose(ctf.eye(k,dtype=numpy.complex64), ctf.dot(VT, ctf.conj(VT.T())))) - - A = ctf.tensor((m,n),dtype=numpy.complex128) - rA = ctf.tensor((m,n),dtype=numpy.float64) - rA.fill_random() - A.real(rA) - iA = ctf.tensor((m,n),dtype=numpy.float64) - iA.fill_random() - A.imag(iA) - - [U,S,VT]=ctf.svd(A,k) - self.assertTrue(allclose(A, ctf.dot(U,ctf.dot(ctf.diag(S),VT)))) - self.assertTrue(allclose(ctf.eye(k,dtype=numpy.complex128), ctf.dot(ctf.conj(U.T()), U))) - self.assertTrue(allclose(ctf.eye(k,dtype=numpy.complex128), ctf.dot(VT, ctf.conj(VT.T())))) - def test_svd_rand(self): m = 19 n = 15 k = 13 - for dt in [numpy.float32, numpy.float64]: + for dt in [numpy.float32, numpy.float64, numpy.complex64, numpy.complex128]: A = ctf.random.random((m,n)) A = ctf.astensor(A,dtype=dt) [U,S,VT]=ctf.svd_rand(A,k,5,1) @@ -124,12 +97,73 @@ def test_svd_rand(self): [U2,S2,VT2]=ctf.svd(A,k) rs1 = ctf.vecnorm(A - ctf.dot(U*S,VT)) rs2 = ctf.vecnorm(A - ctf.dot(U2*S2,VT2)) - self.assertTrue(abs(rs1-rs2)<1.e-2) + rA = ctf.vecnorm(A) + self.assertTrue(rs1 < rA) + self.assertTrue(rs2 < rs1) + self.assertTrue(numpy.abs(rs1 - rs2)<3.e-1) + + def test_tsvd(self): + lens = [4,5,6,3] + for dt in [numpy.float32, numpy.float64, numpy.complex64, numpy.complex128]: + A = ctf.tensor(lens,dtype=dt) + A.fill_random() + [U,S,VT]=A.i("ijkl").svd("ija","akl") + A.i("ijkl") << -1.0*U.i("ija")*S.i("a")*VT.i("akl") + self.assertTrue(ctf.vecnorm(A)/A.tot_size()<1.e-6) + + A = ctf.tensor(lens,dtype=dt) + A.fill_random() + [U,S,VT]=A.i("ijkl").svd("ika","ajl") + A.i("ijkl") << -1.0*U.i("ika")*S.i("a")*VT.i("ajl") + self.assertTrue(ctf.vecnorm(A)/A.tot_size()<1.e-6) + + A = ctf.tensor(lens,dtype=dt) + [U,S1,VT]=A.i("ijkl").svd("ika","ajl",4) + [U,S2,VT]=A.i("ijkl").svd("ika","ajl",numpy.abs(S[4])) + self.assertTrue(allclose(S1,S2)) + + [U,S,VT]=A.i("ijkl").svd("iakj","la") + A.i("ijkl") << -1.0*U.i("iakj")*S.i("a")*VT.i("la") + self.assertTrue(ctf.vecnorm(A)/A.tot_size()<1.e-6) + + A.fill_random() + [U,S,VT]=A.i("ijkl").svd("alk","jai") + A.i("ijkl") << -1.0*U.i("alk")*S.i("a")*VT.i("jai") + self.assertTrue(ctf.vecnorm(A)/A.tot_size()<1.e-6) + K = ctf.tensor((U.shape[0],U.shape[0]),dtype=dt) + K.i("ab") << U.i("alk") * U.i("blk") + self.assertTrue(allclose(K,ctf.eye(U.shape[0]))) + 0.*K.i("ab") << VT.i("jai") * VT.i("jbi") + self.assertTrue(allclose(K,ctf.eye(U.shape[0]))) + + A.fill_random() + [U,S,VT]=A.i("ijkl").svd("alk","jai",4,0,True) + nrm1 = ctf.vecnorm(A) + A.i("ijkl") << -1.0*U.i("alk")*S.i("a")*VT.i("jai") + self.assertTrue(ctf.vecnorm(A)