Actual source code: ex186.c
1: #include <petscmat.h>
3: static char help[] = "Example of MatMat ops with MatDense in PETSc.\n";
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;
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.");
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));
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));
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));
38: // View matrices
39: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix P:\n"));
40: PetscCall(MatView(P, PETSC_VIEWER_STDOUT_SELF));
42: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix A:\n"));
43: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF));
45: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix PtAP:\n"));
46: PetscCall(MatView(PtAP, PETSC_VIEWER_STDOUT_SELF));
48: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix RARt:\n"));
49: PetscCall(MatView(RARt, PETSC_VIEWER_STDOUT_SELF));
51: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Matrix ABC:\n"));
52: PetscCall(MatView(ABC, PETSC_VIEWER_STDOUT_SELF));
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));
62: PetscCall(PetscFinalize());
63: return 0;
64: }
66: /*TEST
68: test:
69: diff_args: -j
70: suffix: 1
72: TEST*/