|
[Sponsors] |
How to add additional, temporary matrices to perform Linear Algebra operations? |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
December 11, 2023, 12:20 |
How to add additional, temporary matrices to perform Linear Algebra operations?
|
#1 |
Senior Member
Klaus
Join Date: Mar 2009
Posts: 260
Rep Power: 22 |
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? |
|
December 11, 2023, 15:11 |
|
#3 |
Senior Member
Klaus
Join Date: Mar 2009
Posts: 260
Rep Power: 22 |
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
Code:
// 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; |
|
December 11, 2023, 17:57 |
|
#5 |
Senior Member
Klaus
Join Date: Mar 2009
Posts: 260
Rep Power: 22 |
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. |
|
December 12, 2023, 03:47 |
|
#6 |
Senior Member
|
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? |
|
December 12, 2023, 04:27 |
|
#7 |
Senior Member
Klaus
Join Date: Mar 2009
Posts: 260
Rep Power: 22 |
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.
Code:
// 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); |
|
December 12, 2023, 07:42 |
|
#8 |
Senior Member
|
Maybe scalarRectangularMatrix?
See L387 to L392 of https://www.openfoam.com/documentati...8C_source.html |
|
December 12, 2023, 12:36 |
|
#9 |
Senior Member
Klaus
Join Date: Mar 2009
Posts: 260
Rep Power: 22 |
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); |
|
December 13, 2023, 03:02 |
|
#10 |
Senior Member
|
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. |
|
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Where and how to create additional, temporary scalarFields? | klausb | OpenFOAM Programming & Development | 0 | December 6, 2023 16:10 |
Interfoam UEqn.H, where to add additional terms? | Luke99 | OpenFOAM Programming & Development | 3 | August 29, 2023 03:26 |
[PyFoam] and paraview | eelcovv | OpenFOAM Community Contributions | 28 | May 30, 2016 09:23 |
emag beta feature: charge density | charlotte | CFX | 4 | March 22, 2011 09:14 |
Replace periodic by inlet-outlet pair | lego | CFX | 3 | November 5, 2002 20:09 |