在C++中使用Octave庫稀疏矩陣乘法（maxmin）

Question

我正在實現一個maxmin函數，它的工作原理與矩陣乘法相似，但它不是求和產品，而是逐點地在兩個數字之間獲取最大值。天真實施的一個例子是

double mx = 0; 
double mn = 0; 
for (i = 0; i < rowsC;i++) 
{ 
    for(j = 0; j < colsC;j++) 
    { 
     mx = 0; 
     for(k = 0; k < colsA; k++) 
     { 
      if (a(i, k) < b(k, j)) 
       mn = a(i,k); 
      else 
       mn = b(k,j); 

      if (mn > mx) 
       mx = mn; 
     } 
     c(i, j) = mx; 
    } 
} 

我把它編碼爲八度oct-file，所以我必須使用oct.h數據結構。問題是，我想實現一個稀疏的版本，但通常你需要到下一個非零元素的引用行或像在這個例子中一列（見4.3算法）： http://www.eecs.harvard.edu/~ellard/Q-97/HTML/root/node20.html

有做row_p - >下一行給出了該行的下一個非零元素（對於列是相同的）。有沒有辦法對八度SparseMatrix類做同樣的事情？還是有另一種實現我可以採用我的maxmin函數的稀疏矩陣乘法的方法？

Answer 1

我不知道anyoe是否會有興趣，但我設法找到一個解決方案。 解決方案的代碼是用於Octave的fl-core1.0模糊邏輯核心軟件包的一部分，並且在LGPL許可下發布。 （該代碼依賴於一些倍頻程功能）

// Calculate the S-Norm/T-Norm composition of sparse matrices (single thread) 
void sparse_compose(octave_value_list args) 
{ 
    // Create constant versions of the input matrices to prevent them to be filled by zeros on reading. 
    // a is the const reference to the transpose of a because octave sparse matrices are column compressed 
    // (to cycle on the rows, we cycle on the columns of the transpose). 
    SparseMatrix atmp = args(0).sparse_matrix_value(); 
    const SparseMatrix a = atmp.transpose(); 
    const SparseMatrix b = args(1).sparse_matrix_value(); 

    // Declare variables for the T-Norm and S-Norm values 
    float snorm_val;  
    float tnorm_val;  

    // Initialize the result sparse matrix 
    sparseC = SparseMatrix((int)colsB, (int)rowsA, (int)(colsB*rowsA)); 

    // Initialize the number of nonzero elements in the sparse matrix c 
    int nel = 0; 
    sparseC.xcidx(0) = 0; 

    // Calculate the composition for each element 
    for (int i = 0; i < rowsC; i++) 
    { 
     for(int j = 0; j < colsC; j++) 
     { 

      // Get the index of the first element of the i-th column of a transpose (i-th row of a) 
      // and the index of the first element of the j-th column of b 
      int ka = a.cidx(i); 
      int kb = b.cidx(j); 
      snorm_val = 0; 

      // Check if the values of the matrix are really not 0 (it happens if the column of a or b hasn't any value) 
      // because otherwise the cidx(i) or cidx(j) returns the first nonzero element of the previous column 
      if(a(a.ridx(ka),i)!=0 && b(b.ridx(kb),j)!=0) 
      { 
       // Cicle on the i-th column of a transpose (i-th row of a) and j-th column of b 
       // From a.cidx(i) to a.cidx(i+1)-1 there are all the nonzero elements of the column i of a transpose (i-th row of a) 
       // From b.cidx(j) to b.cidx(j+1)-1 there are all the nonzero elements of the column j of b 
       while ((ka <= (a.cidx(i+1)-1)) && (kb <= (b.cidx(j+1)-1))) 
       { 

        // If a.ridx(ka) == b.ridx(kb) is true, then there's a nonzero value on the same row 
        // so there's a k for that a'(k, i) (equals to a(i, k)) and b(k, j) are both nonzero 
        if (a.ridx(ka) == b.ridx(kb)) 
        { 
         tnorm_val = calc_tnorm(a.data(ka), b.data(kb)); 
         snorm_val = calc_snorm(snorm_val, tnorm_val); 
         ka++; 
         kb++; 
        } 

        // If a.ridx(ka) == b.ridx(kb) ka should become the index of the next nonzero element on the i column of a 
        // transpose (i row of a) 
        else if (a.ridx(ka) < b.ridx(kb))   
         ka++; 
        // If a.ridx(ka) > b.ridx(kb) kb should become the index of the next nonzero element on the j column of b 
        else 
         kb++; 
       } 
      } 

      if (snorm_val != 0) 
      { 
       // Equivalent to sparseC(i, j) = snorm_val; 
       sparseC.xridx(nel) = j; 
       sparseC.xdata(nel++) = snorm_val; 
      } 
     } 
     sparseC.xcidx(i+1) = nel; 
    } 

    // Compress the result sparse matrix because it is initialized with a number of nonzero element probably greater than the real one 
    sparseC.maybe_compress(); 

    // Transpose the result 
    sparseC = sparseC.transpose(); 
}