// Evgenii B. Rudnyi, http://Evgenii.Rudnyi.Ru #include #include #include extern "C" { #include } using namespace::std; //The class keeps the matrix by rows class SparseMatrix : public vector< map > { void init() { m = 0; } public: size_t m; SparseMatrix() {init();} explicit SparseMatrix(size_t n_) : vector< map >(n_) {init();} void clear() { vector< map >::clear(); init(); } size_t NRows() const {return size();} size_t NCols() const {return m;} void swap(SparseMatrix &y) { vector< map >::swap(y); std::swap(m, y.m); } void clearMemory() { SparseMatrix empty; swap(empty); } }; // C = A A' // A is assumed to be unsymmetric void multiply(SparseMatrix &A, SparseMatrix &C) { // C_ij = Sum_k A_ik*B_kj // B = A' b_kj = a_jk // C_ij = Sum_k A_ik*A_jk // Let us compute only the upper half of C C.clear(); int dim = A.NRows(); C.resize(dim); C.m = dim; for (int i = 0; i < dim; ++i) for (int j = i; j < dim; ++j) { double sum = 0.; if (i == j) { for (map::const_iterator k = A[i].begin(); k != A[i].end(); ++k) sum += (*k).second*(*k).second; } else { for (map::const_iterator k = A[i].begin(); k != A[i].end(); ++k) { map::const_iterator kk = A[j].find((*k).first); if (kk != A[j].end()) sum += (*k).second*(*kk).second; } } if (sum != 0) C[i][j] = sum; } } // Symmetric matrix C above is actually equivalent to CCS required by TAUCS // Attention! This is working only for symmetric matrices taucs_ccs_matrix* convert2taucs(SparseMatrix &C) { taucs_ccs_matrix *A; int dim = C.size(); int nnz = 0; for (int i = 0; i < dim; ++i) { int num = C[i].size(); nnz += num; } A = taucs_dccs_create(dim, dim, nnz); A->flags |= TAUCS_SYMMETRIC | TAUCS_LOWER; int pos = 0; for (int i = 0; i < dim; ++i) { A->colptr[i] = pos; for (map::const_iterator j = C[i].begin(); j != C[i].end(); ++j) { A->rowind[pos] = (*j).first; A->values.d[pos] = (*j).second; ++pos; } } A->colptr[dim] = nnz; return A; }