<div dir="ltr"><div>Hi, Steven</div><div> Thank you for the bug report and test example. You were right. The MatCopy(A, B,..) implementation was wrong when B was on the device.</div><div> I have a fix at <a href="https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/merge_requests/8288__;!!G_uCfscf7eWS!e15LlybplgvD-3BTmMELbBy9EiU1biKGxn_qlVmPVnTPUjMYT9nZ7kf3ZtjPPhbQL7zFDI8QGuKk_V-tW9IYCL3xFVri$">https://gitlab.com/petsc/petsc/-/merge_requests/8288</a>, and will add your test to the MR tomorrow.</div><div><br></div><div> Thanks!</div><div><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr">--Junchao Zhang</div></div></div><br></div><br><div class="gmail_quote gmail_quote_container"><div dir="ltr" class="gmail_attr">On Mon, Apr 7, 2025 at 4:57 PM Steven Dargaville <<a href="mailto:dargaville.steven@gmail.com">dargaville.steven@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Hi</div><div><br></div><div>I have some code that computes a MatMatMult, then modifies one of the matrices and calls the MatMatMult again with MAT_REUSE_MATRIX. This gives the correct results with mat_types aij, aijkokkos and either run on a CPU or a GPU.</div><div><br></div><div>In some cases one of the matrices is diagonal so I have been modifying my code to call MatDiagonalScale instead and have seen some peculiar behaviour. If the mat type is aij, everything works correctly on the CPU and GPU. If the mat type is aijkokkos, everything works correctly on the CPU but when run on a machine with a GPU the results differ. <br></div><div><br></div><div>I have some example code below that shows this, it prints out four matrices; the first two should match and the last two should match. To see the failure run with "-mat_type aijkokkos" on a machine with a GPU (again this gives the correct results with aijkokkos run on a CPU). I get the output:</div><div><br></div><div>Mat Object: 1 MPI process<br> type: seqaijkokkos<br>row 0: (0, 4.) (1, 6.) <br>row 1: (0, 10.) (1, 14.) <br>Mat Object: 1 MPI process<br> type: seqaijkokkos<br>row 0: (0, 4.) (1, 6.) <br>row 1: (0, 10.) (1, 14.) <br>Mat Object: 1 MPI process<br> type: seqaijkokkos<br>row 0: (0, 6.) (1, 9.) <br>row 1: (0, 15.) (1, 21.) <br>Mat Object: 1 MPI process<br> type: seqaijkokkos<br>row 0: (0, 12.) (1, 18.) <br>row 1: (0, 30.) (1, 42.)</div><div><br></div><div>I have narrowed down this failure to a MatCopy call, where the values of result_diag should be overwritten with A before calling MatDiagonalScale. The results with aijkokos on a GPU suggest the values of result_diag are not being changed. If instead of using the MatCopy I destroy result_diag and call MatDuplicate, the results are correct. You can trigger the correct behaviour with "-mat_type aijkokkos -correct". <br></div><div><br></div><div>To me this looks like the values are out of sync between the device/host, but I would have thought a call to MatCopy would be safe. Any thoughts on what I might be doing wrong?</div><div><br></div><div>Thanks for your help!</div><div>Steven</div><div><br></div><div>~~~~~~~~~~~~~~~~~~~~~~~~~~~~~<br><br>static char help[] = "Test matmat products with matdiagonal on gpus \n\n";<br><br>#include <iostream><br>#include <petscmat.h><br><br><br>int main(int argc, char **args)<br>{<br> const PetscInt inds[] = {0, 1};<br> PetscScalar avals[] = {2, 3, 5, 7};<br> Mat A, B_diag, B_aij_diag, result, result_diag;<br> Vec diag;<br> <br> PetscCall(PetscInitialize(&argc, &args, NULL, help));<br><br> // Create matrix to start<br> PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 2, 2, 2, 2, &A));<br> PetscCall(MatSetUp(A));<br> PetscCall(MatSetValues(A, 2, inds, 2, inds, avals, INSERT_VALUES));<br> PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));<br> PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); <br><br> // Create a matdiagonal matrix<br> // Will be the matching vec type as A<br> PetscCall(MatCreateVecs(A, &diag, NULL));<br> PetscCall(VecSet(diag, 2.0));<br> PetscCall(MatCreateDiagonal(diag, &B_diag));<br><br> // Create the same matrix as the matdiagonal but in aij format<br> PetscCall(MatCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, 2, 2, 2, 2, &B_aij_diag));<br> PetscCall(MatSetUp(B_aij_diag));<br> PetscCall(MatDiagonalSet(B_aij_diag, diag, INSERT_VALUES));<br> PetscCall(MatAssemblyBegin(B_aij_diag, MAT_FINAL_ASSEMBLY));<br> PetscCall(MatAssemblyEnd(B_aij_diag, MAT_FINAL_ASSEMBLY)); <br> PetscCall(VecDestroy(&diag));<br><br> // ~~~~~~~~~~~~~<br> // Do an initial matmatmult<br> // A * B_aij_diag<br> // and then<br> // A * B_diag but just using MatDiagonalScale<br> // ~~~~~~~~~~~~~<br><br> // aij * aij<br> PetscCall(MatMatMult(A, B_aij_diag, MAT_INITIAL_MATRIX, 1.5, &result));<br> PetscCall(MatView(result, PETSC_VIEWER_STDOUT_WORLD));<br><br> // aij * diagonal <br> PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &result_diag));<br> PetscCall(MatDiagonalGetDiagonal(B_diag, &diag));<br> PetscCall(MatDiagonalScale(result_diag, NULL, diag));<br> PetscCall(MatDiagonalRestoreDiagonal(B_diag, &diag));<br> PetscCall(MatView(result_diag, PETSC_VIEWER_STDOUT_WORLD));<br><br> // ~~~~~~~~~~~~~<br> // Now let's modify the diagonal and do it again with "reuse"<br> // ~~~~~~~~~~~~~ <br> PetscCall(MatDiagonalGetDiagonal(B_diag, &diag));<br> PetscCall(VecSet(diag, 3.0));<br> PetscCall(MatDiagonalSet(B_aij_diag, diag, INSERT_VALUES));<br> PetscCall(MatDiagonalRestoreDiagonal(B_diag, &diag));<br><br> // aij * aij<br> PetscCall(MatMatMult(A, B_aij_diag, MAT_REUSE_MATRIX, 1.5, &result));<br> PetscCall(MatView(result, PETSC_VIEWER_STDOUT_WORLD));<br><br> PetscBool correct = PETSC_FALSE;<br> PetscCall(PetscOptionsGetBool(NULL, NULL, "-correct", &correct, NULL));<br> <br> if (!correct)<br> {<br> // This gives the wrong results below when run on gpu<br> // Results suggest this isn't copied <br> PetscCall(MatCopy(A, result_diag, SAME_NONZERO_PATTERN));<br> }<br> else<br> {<br> // This gives the correct results below<br> PetscCall(MatDestroy(&result_diag));<br> PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &result_diag));<br> }<br> <br> // aij * diagonal <br> PetscCall(MatDiagonalGetDiagonal(B_diag, &diag));<br> PetscCall(MatDiagonalScale(result_diag, NULL, diag));<br> PetscCall(MatDiagonalRestoreDiagonal(B_diag, &diag)); <br> PetscCall(MatView(result_diag, PETSC_VIEWER_STDOUT_WORLD));<br><br> PetscCall(PetscFinalize());<br> return 0;<br>}</div></div>
</blockquote></div>