2
votes

Matlab code here

 N = 4096;
 x       =   linspace(-1,1,N);     
 A = sparse(100000,100000); 
 index = 1:32;

 A(index,index)    =   kernel(x(index),x(index));
  //kernel here outputs a 32 x 32 matrix 

I have to convert the MATLAB code to C++,but I am stuck at sparse function, I tried using:

 N=4096;
 Eigen::VectorXd x = Eigen::VectorXd::LinSpaced(N,-1,1); 
 Eigen::SparseMatrix<double> A(Asize,Asize);
 A.block(1,1,index.size(), index.size()) = Kernel();

But SparseMatrix has block as read only function and therefore can't be useful for updating the matrix.

Another Point:

I have gone through the Eigen documentation and checked a different form of SparseMatrix declaration :

 typedef Eigen::Triplet<double> T;
 std::vector<T>tripleList;
 tripleList.reserve(nnz);

 for(...)
 {
  // ...
  tripletList.push_back(T(i,j,v_ij));   //?? what are these values?
  }
 A.setFromTriplets(tripleList.begin(), tripleList.end());

But I don't understand what should be value of nnz, should it be the value that I get from Matlab code and what value should I push through the for loop? Will they be random, how do I select the "pushed values" given that matrix size is so large.

Also, the last question remains, how will the sparse matrix after declaration gets updated block-wise?

1
nnz means "number of non zeroes". I assume that is for memory allocation in Eigen. Additionally, ` A(index,index) = kernel(x(index),x(index));` this can be done with a simple double nested loop in C++, you dont need a one-linerAnder Biguri
According to the documentaion for performance reasons, writing to a sub-sparse-matrix is limited. and only you can update contiguous sets of columns. with the size of submatrix should be known at compile timerahnema1
@AnderBiguri How can it be done in a double nested loop? are u suggesting using the function "insert"? but then if the same line is repeated multiple times I have to first store the result from Kernel() every time and then insert it in sparse matrix.user7440094

1 Answers

3
votes

In general, Eigen sparse matrix blocks are const. The exception to this rule is e.g. m.col(i) in a column major matrix or .row(i) in a row major one.

As for your second point, nnz is the number of non-zeros. By specifying how many non-zeros the matrix needs to hold at the beginning, you minimize the number of reallocations/copies needed when constructing the sparse matrix. See the documentation. Also, note that other overload is useful if you know how many elements are needed in each column. Regarding the for loop, assume you had three arrays: columnIndices, rowIndices and values, all of length nnz. Your for loop would look like:

for(int i = 0; i < nnz; i++)
{
    tripletList.push_back(T(columnIndices[i], rowIndices[i], values[i]));
}

Alternatively, you could compute the values for each triplet within the for loop.