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