[klausb]

I need additional, temporary matrices to perform Linear Algebra operations including matrix-matrix multiplications, subtractions and additions. The content of the matrices differs from the coefficient matrix and could be e.g. a subset like the (strictly) lower part. There are multiple matrix classes like lduMatrix or RectangularMatrix.

Which one should be used and how would I construct a new, empty one and e.g. add the lower part or the coefficient matrix?

[dlahaye]

Please elaborate. Are your matrices square or rectangular, sparse or full, complex or real valued, for scalar fields or vector fields?

[klausb]

See my code below, right now, I can't figure out how to do the actual matrix operations *, +, - like G=E*F or G=E+F or G=E-F

    // Set the size of the matrices
    const label N = 4; // matrix size
    
    // Declare and initialize lduMatrices E, F, and G
    RectangularMatrix <scalar> E(N,N);
    RectangularMatrix <scalar> F(N,N);
    RectangularMatrix <scalar> G(N,N);

    // Set coefficients for matrix E diagonal
    scalar diagonalValuesE[] = {1.0, 2.0, 3.0, 4.0};
    for (label i = 0; i < N; ++i)
    {
        E[i][i] = diagonalValuesE[i];
    }

    // Set coefficients for matrix F diagonal
    scalar diagonalValuesF[] = {5.0, 6.0, 7.0, 8.0};
    for (label i = 0; i < N; ++i)
    {
        F[i][i] = diagonalValuesF[i];
    }

    //Compute G = E * F ----> things don't work from here onwards
    G=E*F; 
    
    // Print the diagonal element G[0][0]
    Info << "Diagonal element G[0][0]: " << G[0][0] << endl;

// https://www.openfoam.com/documentati...trices_8H.html

[dlahaye]

1/ does printing of E and F work fine?

2/ should this be multiply(G, E, F) instead?

[klausb]

1/ does printing of E and F work fine? YES

2/ should this be multiply(G, E, F) instead? NO, that's not supported by the RectangularMatrix class


The Matrix class seems to support multiply(G, E, F) but the parameters are different and I couldn't figure out how to define a Matrix instead of a RectangularMatrix.

[dlahaye]

I see an implementation of multiply() starting at L238 of scalarMatricesTemplaces.C at https://www.openfoam.com/documentati...8C_source.html

I therefore imagine that it might be a matter of given scalarMatrix a suitable form and type.

Does this help?

[klausb]

I am still stuck with an error message "multiply, there's no valid function", I tried various matrix types but I can't get a matrix x matrix multiplication to work.

// Set the size of the matrices
    const label N = 4; // matrix size
    
    // Declare and initialize lduMatrices E, F, and G RectangularMatrix
    RectangularMatrix <scalar> E(N,N); 
    RectangularMatrix <scalar> F(N,N);
    RectangularMatrix <scalar> G(N,N);
    
    Matrix <scalarRectangularMatrix, scalar> J(N,N); 
    Matrix <scalarRectangularMatrix, scalar> K(N,N);
    Matrix <scalarRectangularMatrix, scalar> L(N,N);
    

    // Set coefficients for matrix E diagonal
    scalar diagonalValuesE[] = {1.0, 2.0, 3.0, 4.0};
    for (label i = 0; i < N; ++i)
    {
        E[i][i] = diagonalValuesE[i];
    }
    
    Info << "element from E[1][1]: " << E[1][1] << endl;
    
    // Set coefficients for matrix F diagonal
    scalar diagonalValuesF[] = {5.0, 6.0, 7.0, 8.0};
    for (label i = 0; i < N; ++i)
    {
        F[i][i] = diagonalValuesF[i];
    }
    
    Info << "element from F[1][1]: " << F[1][1] << endl;
    
    //-----------------------------------------------------
    
    // Set coefficients for matrix J diagonal
    scalar diagonalValuesJ[] = {5.0, 6.0, 7.0, 8.0};
    for (label i = 0; i < N; ++i)
    {
        J[i][i] = diagonalValuesJ[i];
    }
    
    Info << "element from J[1][1]: " << J[1][1] << endl;
    
        
    // Set coefficients for matrix K diagonal
    scalar diagonalValuesK[] = {5.0, 6.0, 7.0, 8.0};
    for (label i = 0; i < N; ++i)
    {
        K[i][i] = diagonalValuesK[i];
    }
    
    Info << "element from K[1][1]: " << K[1][1] << endl;
    
    
    // J = J * K ; note: the middle position would be for a matrix defined as Diagonal
    multiply(J,Zero,K);

[dlahaye]

Maybe scalarRectangularMatrix?

See L387 to L392 of

https://www.openfoam.com/documentati...8C_source.html

[klausb]

multiply(L,J,K); is still the problem

It seems, the matrix in the middle is supposed to be a DiagonalMatrix so I changed it but it's still not working. There are settings that compile but the application (I test it adding the code to a preconditioner) will crash when the computation reaches multiply(L,J,K);

[dlahaye]

1/ My always limited understanding is that OpenFoam does *not* provide sparse-matrix-times-sparse-matrix as the multiplication of sparse matrices is tricky and typically not required.

Assume A and B to be sparse matrices and v to be a vector.

Typically one requires (A * B) * v.

However, (A * B) * v = A * (B * v) and only the matrix-vector multiplication with A and B is required.

2/ Diagonal preconditioning is an application of above. Indeed

Assume A x = b to be a linear system to be solved for x. Assume D = diagonal(A). Then diagonal preconditioning can be implemented by adding the operation
z = inv(D) * r

at each step of the outer Krylov (CG, GMRES or other). The computation of inv(D) * A is avoided.

3/ I do agreed that inv(D) * A is attractive to build explicitly in a testing phase. This is referred to as diagonal scaling. Not sure whether OpenFoam has is. Possibly looking into the implementation of the diagonal preconditioner gives clues.