Actual source code: ex258.c

  1: static char help[] = "Test MatProductReplaceMats() \n\
  2: Modified from the code contributed by Pierre Jolivet \n\n";

  4: #include <petscmat.h>

  6: int main(int argc, char **args)
  7: {
  8:   PetscInt  n = 2, convert;
  9:   Mat       A, B, Bdense, Conjugate;
 10:   PetscBool conjugate = PETSC_FALSE, equal, flg;

 12:   PetscFunctionBeginUser;
 13:   PetscCall(PetscInitialize(&argc, &args, NULL, help));

 15:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 16:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
 17:   PetscCall(MatSetType(A, MATDENSE));
 18:   PetscCall(MatSetFromOptions(A));
 19:   PetscCall(MatSeqDenseSetPreallocation(A, NULL));
 20:   PetscCall(MatMPIDenseSetPreallocation(A, NULL));
 21:   PetscCall(MatSetRandom(A, NULL));
 22:   PetscCall(MatViewFromOptions(A, NULL, "-A_view"));
 23:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-conjugate", &conjugate, NULL));

 25:   for (convert = 0; convert < 2; convert++) {
 26:     /* convert dense matrix A to aij format */
 27:     if (convert) PetscCall(MatConvert(A, MATAIJ, MAT_INPLACE_MATRIX, &A));

 29:     /* compute B = A^T * A or  B = A^H * A */
 30:     PetscCall(MatProductCreate(A, A, NULL, &B));

 32:     flg = PETSC_FALSE;
 33:     PetscCall(PetscOptionsGetBool(NULL, NULL, "-atb", &flg, NULL));
 34:     if (flg) {
 35:       PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
 36:     } else {
 37:       PetscCall(PetscOptionsGetBool(NULL, NULL, "-ptap", &flg, NULL));
 38:       if (flg) {
 39:         PetscCall(MatProductSetType(B, MATPRODUCT_PtAP));
 40:       } else {
 41:         PetscCall(PetscOptionsGetBool(NULL, NULL, "-abt", &flg, NULL));
 42:         if (flg) {
 43:           PetscCall(MatProductSetType(B, MATPRODUCT_ABt));
 44:         } else {
 45:           PetscCall(MatProductSetType(B, MATPRODUCT_AB));
 46:         }
 47:       }
 48:     }
 49:     PetscCall(MatProductSetFromOptions(B));
 50:     PetscCall(MatProductSymbolic(B));

 52:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &Conjugate));
 53:     if (conjugate) PetscCall(MatConjugate(Conjugate));

 55:     /* replace input A by Conjugate */
 56:     PetscCall(MatProductReplaceMats(Conjugate, NULL, NULL, B));

 58:     PetscCall(MatProductNumeric(B));
 59:     PetscCall(MatViewFromOptions(B, NULL, "-product_view"));

 61:     PetscCall(MatDestroy(&Conjugate));
 62:     if (!convert) {
 63:       Bdense = B;
 64:       B      = NULL;
 65:     }
 66:   }

 68:   /* Compare Bdense and B */
 69:   PetscCall(MatMultEqual(Bdense, B, 10, &equal));
 70:   PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Bdense != B");

 72:   PetscCall(MatDestroy(&Bdense));
 73:   PetscCall(MatDestroy(&B));
 74:   PetscCall(MatDestroy(&A));
 75:   PetscCall(PetscFinalize());
 76:   return 0;
 77: }

 79: /*TEST

 81:    test:
 82:       suffix: 1
 83:       args: -conjugate false -atb
 84:       output_file: output/ex258_1.out

 86:    test:
 87:       suffix: 2
 88:       args: -conjugate true -atb
 89:       output_file: output/ex258_1.out

 91:    test:
 92:       suffix: 3
 93:       args: -conjugate false
 94:       output_file: output/ex258_1.out

 96:    test:
 97:       suffix: 4
 98:       args: -ptap
 99:       output_file: output/ex258_1.out

101:    test:
102:       suffix: 5
103:       args: -abt
104:       output_file: output/ex258_1.out

106:    test:
107:       suffix: 6
108:       nsize: 2
109:       args: -conjugate false -atb
110:       output_file: output/ex258_1.out

112:    test:
113:       suffix: 7
114:       nsize: 2
115:       args: -conjugate true -atb
116:       output_file: output/ex258_1.out

118:    test:
119:       suffix: 8
120:       nsize: 2
121:       args: -conjugate false
122:       output_file: output/ex258_1.out

124:    test:
125:       suffix: 9
126:       nsize: 2
127:       args: -ptap
128:       output_file: output/ex258_1.out

130:    test:
131:       suffix: 10
132:       nsize: 2
133:       args: -abt
134:       output_file: output/ex258_1.out

136:    test:
137:       suffix: 11
138:       nsize: 2
139:       args: -conjugate true -atb -mat_product_algorithm backend
140:       TODO: bug: MatProductReplaceMats() does not change the product for this test

142: TEST*/