1 #include <petscmat.h> 2 3 static char help[] = "Example of MatMat ops with MatDense in PETSc.\n"; 4 5 int main(int argc, char **args) 6 { 7 Mat A, P, PtAP, RARt, ABC, Pt; 8 PetscInt n = 4, m = 2; // Example dimensions 9 PetscMPIInt size; 10 11 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 12 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 13 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This example is for sequential runs only."); 14 15 // Create dense matrix P (n x m) 16 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n, m, NULL, &P)); 17 PetscScalar P_values[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0}; 18 PetscInt rows[] = {0, 1, 2, 3}; 19 PetscInt cols[] = {0, 1}; 20 PetscCall(MatSetValues(P, n, rows, m, cols, P_values, INSERT_VALUES)); 21 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY)); 22 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY)); 23 PetscCall(MatTranspose(P, MAT_INITIAL_MATRIX, &Pt)); 24 25 // Create dense matrix A (n x n) 26 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, n, n, NULL, &A)); 27 PetscCall(MatSetBlockSize(A, n)); // Set block size for A 28 PetscScalar A_values[] = {4.0, 1.0, 2.0, 0.0, 1.0, 3.0, 0.0, 1.0, 2.0, 0.0, 5.0, 2.0, 0.0, 1.0, 2.0, 6.0}; 29 PetscInt indices[] = {0, 1, 2, 3}; 30 PetscCall(MatSetValues(A, n, indices, n, indices, A_values, INSERT_VALUES)); 31 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 32 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 33 34 PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &PtAP)); 35 PetscCall(MatRARt(A, Pt, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &RARt)); 36 PetscCall(MatMatMatMult(Pt, A, P, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &ABC)); 37 38 // View matrices 39 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix P:\n")); 40 PetscCall(MatView(P, PETSC_VIEWER_STDOUT_SELF)); 41 42 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix A:\n")); 43 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF)); 44 45 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix PtAP:\n")); 46 PetscCall(MatView(PtAP, PETSC_VIEWER_STDOUT_SELF)); 47 48 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix RARt:\n")); 49 PetscCall(MatView(RARt, PETSC_VIEWER_STDOUT_SELF)); 50 51 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix ABC:\n")); 52 PetscCall(MatView(ABC, PETSC_VIEWER_STDOUT_SELF)); 53 54 // Clean up 55 PetscCall(MatDestroy(&P)); 56 PetscCall(MatDestroy(&Pt)); 57 PetscCall(MatDestroy(&A)); 58 PetscCall(MatDestroy(&PtAP)); 59 PetscCall(MatDestroy(&RARt)); 60 PetscCall(MatDestroy(&ABC)); 61 62 PetscCall(PetscFinalize()); 63 return 0; 64 } 65 66 /*TEST 67 68 test: 69 diff_args: -j 70 suffix: 1 71 72 TEST*/ 73