Sparse Matrix Multiplikation wie (Maxmin) in C++ mit Oktavbibliotheken

Question

Ich implementiere eine Maxmin-Funktion, es funktioniert wie Matrix-Multiplikation, aber anstelle der Summierung von Produkten erhält es maximal von min zwischen zwei Zahlen punktweise. Ein Beispiel für eine naive Implementierung ist

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; 
    } 
} 

ich es als eine Oktave oct-Datei bin Codierung so i oct.h Datenstruktur verwenden. Das Problem ist, dass ich eine spärliche Version implementieren möchten, aber in der Regel müssen Sie einen Verweis auf das nächste nicht Nullelement in einer Zeile oder in einer Spalte wie in diesem Beispiel (4.3 Algorithmus e): http://www.eecs.harvard.edu/~ellard/Q-97/HTML/root/node20.html

Es row_p tun -> next gab das nächste von Null verschiedene Element der Zeile (das gleiche für die Spalte). Gibt es eine Möglichkeit, dasselbe mit der Octave SparseMatrix-Klasse zu tun? Oder gibt es eine andere Möglichkeit, die Multiplikation der Sparse-Matrix zu implementieren, die ich für meine Maxmin-Funktion übernehmen kann?

Answer 1

Ich weiß nicht, ob anyoe jemals daran interessiert wäre, aber ich habe es geschafft, eine Lösung zu finden. Der Code der Lösung ist Teil von fl-core1.0, einem Fuzzy-Logic-Kernpaket für Octave, das unter der LGPL-Lizenz veröffentlicht wurde. (Der Code beruht auf einigen Oktavenfunktionen)

// 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(); 
}