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*/