Actual source code: ex90.c

  1: static char help[] = "Tests MatPtAP() \n";

  3: #include <petscmat.h>

  5: /*
  6:  * This is an extremely simple example to test MatPtAP. It is very useful when developing and debugging the code.
  7:  *
  8:  * A =

 10:    1   2   0   4
 11:    0   1   2   0
 12:    2   0   4   0
 13:    0   1   2   1
 14:  *
 15:  *
 16:  *
 17:  * P =

 19:    1.00000   0.00000
 20:    0.30000   0.50000
 21:    0.00000   0.80000
 22:    0.90000   0.00000
 23:  *
 24:  *
 25:  *AP =
 26:  *
 27:  * 5.20000   1.00000
 28:    0.30000   2.10000
 29:    2.00000   3.20000
 30:    1.20000   2.10000
 31:  *
 32:  * PT =

 34:    1.00000   0.30000   0.00000   0.90000
 35:    0.00000   0.50000   0.80000   0.00000

 37:  *
 38:  * C =

 40:    6.3700   3.5200
 41:    1.7500   3.6100
 42:  *
 43:  * */

 45: int main(int argc, char **argv)
 46: {
 47:   Mat         A, P, PtAP;
 48:   PetscInt    i1[] = {0, 3, 5}, i2[] = {0, 2, 5};
 49:   PetscInt    j1[] = {0, 1, 3, 1, 2}, j2[] = {0, 2, 1, 2, 3};
 50:   PetscScalar a1[] = {1, 2, 4, 1, 2}, a2[] = {2, 4, 1, 2, 1};
 51:   PetscInt    pi1[] = {0, 1, 3}, pi2[] = {0, 1, 2};
 52:   PetscInt    pj1[] = {0, 0, 1}, pj2[] = {1, 0};
 53:   PetscScalar pa1[] = {1, 0.3, 0.5}, pa2[] = {0.8, 0.9};
 54:   MPI_Comm    comm;
 55:   PetscMPIInt rank, size;

 57:   PetscFunctionBeginUser;
 58:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 59:   comm = PETSC_COMM_WORLD;
 60:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
 61:   PetscCallMPI(MPI_Comm_size(comm, &size));
 62:   PetscCheck(size == 2, comm, PETSC_ERR_WRONG_MPI_SIZE, "You have to use two processor cores to run this example ");
 63:   PetscCall(MatCreateMPIAIJWithArrays(comm, 2, 2, PETSC_DETERMINE, PETSC_DETERMINE, rank ? i2 : i1, rank ? j2 : j1, rank ? a2 : a1, &A));
 64:   PetscCall(MatCreateMPIAIJWithArrays(comm, 2, 1, PETSC_DETERMINE, PETSC_DETERMINE, rank ? pi2 : pi1, rank ? pj2 : pj1, rank ? pa2 : pa1, &P));
 65:   PetscCall(MatPtAP(A, P, MAT_INITIAL_MATRIX, 1.1, &PtAP));
 66:   PetscCall(MatView(A, NULL));
 67:   PetscCall(MatView(P, NULL));
 68:   PetscCall(MatView(PtAP, NULL));
 69:   PetscCall(MatPtAP(A, P, MAT_REUSE_MATRIX, 1.1, &PtAP));
 70:   PetscCall(MatView(A, NULL));
 71:   PetscCall(MatView(P, NULL));
 72:   PetscCall(MatView(PtAP, NULL));
 73:   PetscCall(MatDestroy(&A));
 74:   PetscCall(MatDestroy(&P));
 75:   PetscCall(MatDestroy(&PtAP));
 76:   PetscCall(PetscFinalize());
 77:   return 0;
 78: }

 80: /*TEST
 81:    test:
 82:      nsize: 2
 83:      args: -matptap_via allatonce
 84:      output_file: output/ex90_1.out

 86:    test:
 87:      nsize: 2
 88:      suffix: merged
 89:      args: -matptap_via allatonce_merged
 90:      output_file: output/ex90_1.out

 92: TEST*/