1 /* 2 * ex195.c 3 * 4 * Created on: Aug 24, 2015 5 * Author: Fande Kong <fdkong.jd@gmail.com> 6 */ 7 8 static char help[] = " Demonstrate the use of MatConvert_Nest_AIJ\n"; 9 10 #include <petscmat.h> 11 12 int main(int argc, char **args) 13 { 14 Mat A1, A2, A3, A4, A5, B, C, C1, nest; 15 Mat mata[4]; 16 Mat aij; 17 MPI_Comm comm; 18 PetscInt m, M, n, istart, iend, ii, i, J, j, K = 10; 19 PetscScalar v; 20 PetscMPIInt size; 21 PetscBool equal, test = PETSC_FALSE; 22 23 PetscFunctionBeginUser; 24 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 25 comm = PETSC_COMM_WORLD; 26 PetscCallMPI(MPI_Comm_size(comm, &size)); 27 28 /* 29 Assemble the matrix for the five point stencil, YET AGAIN 30 */ 31 PetscCall(MatCreate(comm, &A1)); 32 m = 2, n = 2; 33 PetscCall(MatSetSizes(A1, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n)); 34 PetscCall(MatSetFromOptions(A1)); 35 PetscCall(MatSetUp(A1)); 36 PetscCall(MatGetOwnershipRange(A1, &istart, &iend)); 37 for (ii = istart; ii < iend; ii++) { 38 v = -1.0; 39 i = ii / n; 40 j = ii - i * n; 41 if (i > 0) { 42 J = ii - n; 43 PetscCall(MatSetValues(A1, 1, &ii, 1, &J, &v, INSERT_VALUES)); 44 } 45 if (i < m - 1) { 46 J = ii + n; 47 PetscCall(MatSetValues(A1, 1, &ii, 1, &J, &v, INSERT_VALUES)); 48 } 49 if (j > 0) { 50 J = ii - 1; 51 PetscCall(MatSetValues(A1, 1, &ii, 1, &J, &v, INSERT_VALUES)); 52 } 53 if (j < n - 1) { 54 J = ii + 1; 55 PetscCall(MatSetValues(A1, 1, &ii, 1, &J, &v, INSERT_VALUES)); 56 } 57 v = 4.0; 58 PetscCall(MatSetValues(A1, 1, &ii, 1, &ii, &v, INSERT_VALUES)); 59 } 60 PetscCall(MatAssemblyBegin(A1, MAT_FINAL_ASSEMBLY)); 61 PetscCall(MatAssemblyEnd(A1, MAT_FINAL_ASSEMBLY)); 62 PetscCall(MatView(A1, PETSC_VIEWER_STDOUT_WORLD)); 63 64 PetscCall(MatDuplicate(A1, MAT_COPY_VALUES, &A2)); 65 PetscCall(MatDuplicate(A1, MAT_COPY_VALUES, &A3)); 66 PetscCall(MatDuplicate(A1, MAT_COPY_VALUES, &A4)); 67 68 /*create a nest matrix */ 69 PetscCall(MatCreate(comm, &nest)); 70 PetscCall(MatSetType(nest, MATNEST)); 71 mata[0] = A1, mata[1] = A2, mata[2] = A3, mata[3] = A4; 72 PetscCall(MatNestSetSubMats(nest, 2, NULL, 2, NULL, mata)); 73 PetscCall(MatSetUp(nest)); 74 75 /* test matrix product error messages */ 76 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nest*nest", &test, NULL)); 77 if (test) PetscCall(MatMatMult(nest, nest, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C)); 78 79 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matproductsetfromoptions", &test, NULL)); 80 if (test) { 81 PetscCall(MatProductCreate(nest, nest, NULL, &C)); 82 PetscCall(MatProductSetType(C, MATPRODUCT_AB)); 83 PetscCall(MatProductSymbolic(C)); 84 } 85 86 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matproductsymbolic", &test, NULL)); 87 if (test) { 88 PetscCall(MatProductCreate(nest, nest, NULL, &C)); 89 PetscCall(MatProductSetType(C, MATPRODUCT_AB)); 90 PetscCall(MatProductSetFromOptions(C)); 91 PetscCall(MatProductSymbolic(C)); 92 } 93 94 PetscCall(MatConvert(nest, MATAIJ, MAT_INITIAL_MATRIX, &aij)); 95 PetscCall(MatView(aij, PETSC_VIEWER_STDOUT_WORLD)); 96 97 /* create a dense matrix */ 98 PetscCall(MatGetSize(nest, &M, NULL)); 99 PetscCall(MatGetLocalSize(nest, &m, NULL)); 100 PetscCall(MatCreateDense(comm, m, PETSC_DECIDE, M, K, NULL, &B)); 101 PetscCall(MatSetRandom(B, NULL)); 102 103 /* C = nest*B_dense */ 104 PetscCall(MatMatMult(nest, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C)); 105 PetscCall(MatMatMult(nest, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C)); 106 PetscCall(MatMatMultEqual(nest, B, C, 10, &equal)); 107 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in C != nest*B_dense"); 108 109 /* Test B = nest*C, reuse C and B with MatProductCreateWithMat() */ 110 /* C has been obtained from nest*B. Clear internal data structures related to factors to prevent circular references */ 111 PetscCall(MatProductClear(C)); 112 PetscCall(MatProductCreateWithMat(nest, C, NULL, B)); 113 PetscCall(MatProductSetType(B, MATPRODUCT_AB)); 114 PetscCall(MatProductSetFromOptions(B)); 115 PetscCall(MatProductSymbolic(B)); 116 PetscCall(MatProductNumeric(B)); 117 PetscCall(MatMatMultEqual(nest, C, B, 10, &equal)); 118 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in B != nest*C_dense"); 119 PetscCall(MatConvert(nest, MATAIJ, MAT_INPLACE_MATRIX, &nest)); 120 PetscCall(MatEqual(nest, aij, &equal)); 121 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in aij != nest"); 122 PetscCall(MatDestroy(&nest)); 123 124 if (size > 1) { /* Do not know why this test fails for size = 1 */ 125 PetscCall(MatCreateTranspose(A1, &A5)); /* A1 is symmetric */ 126 mata[0] = A5; 127 PetscCall(MatCreate(comm, &nest)); 128 PetscCall(MatSetType(nest, MATNEST)); 129 PetscCall(MatNestSetSubMats(nest, 2, NULL, 2, NULL, mata)); 130 PetscCall(MatSetUp(nest)); 131 PetscCall(MatMatMult(nest, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C1)); 132 PetscCall(MatMatMult(nest, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C1)); 133 134 PetscCall(MatMatMultEqual(nest, B, C1, 10, &equal)); 135 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in C1 != C"); 136 PetscCall(MatDestroy(&C1)); 137 PetscCall(MatDestroy(&A5)); 138 PetscCall(MatDestroy(&nest)); 139 } 140 141 PetscCall(MatDestroy(&C)); 142 PetscCall(MatDestroy(&B)); 143 PetscCall(MatDestroy(&aij)); 144 PetscCall(MatDestroy(&A1)); 145 PetscCall(MatDestroy(&A2)); 146 PetscCall(MatDestroy(&A3)); 147 PetscCall(MatDestroy(&A4)); 148 PetscCall(PetscFinalize()); 149 return 0; 150 } 151 152 /*TEST 153 154 test: 155 nsize: 2 156 157 test: 158 suffix: 2 159 160 TEST*/ 161