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; 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 PetscCall(MatConvert(nest, MATAIJ, MAT_INITIAL_MATRIX, &aij)); 75 PetscCall(MatView(aij, PETSC_VIEWER_STDOUT_WORLD)); 76 77 /* create a dense matrix */ 78 PetscCall(MatGetSize(nest, &M, NULL)); 79 PetscCall(MatGetLocalSize(nest, &m, NULL)); 80 PetscCall(MatCreateDense(comm, m, PETSC_DECIDE, M, K, NULL, &B)); 81 PetscCall(MatSetRandom(B, NULL)); 82 83 /* C = nest*B_dense */ 84 PetscCall(MatMatMult(nest, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C)); 85 PetscCall(MatMatMult(nest, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C)); 86 PetscCall(MatMatMultEqual(nest, B, C, 10, &equal)); 87 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in C != nest*B_dense"); 88 89 /* Test B = nest*C, reuse C and B with MatProductCreateWithMat() */ 90 /* C has been obtained from nest*B. Clear internal data structures related to factors to prevent circular references */ 91 PetscCall(MatProductClear(C)); 92 PetscCall(MatProductCreateWithMat(nest, C, NULL, B)); 93 PetscCall(MatProductSetType(B, MATPRODUCT_AB)); 94 PetscCall(MatProductSetFromOptions(B)); 95 PetscCall(MatProductSymbolic(B)); 96 PetscCall(MatProductNumeric(B)); 97 PetscCall(MatMatMultEqual(nest, C, B, 10, &equal)); 98 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in B != nest*C_dense"); 99 PetscCall(MatConvert(nest, MATAIJ, MAT_INPLACE_MATRIX, &nest)); 100 PetscCall(MatEqual(nest, aij, &equal)); 101 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in aij != nest"); 102 PetscCall(MatDestroy(&nest)); 103 104 if (size > 1) { /* Do not know why this test fails for size = 1 */ 105 PetscCall(MatCreateTranspose(A1, &A5)); /* A1 is symmetric */ 106 mata[0] = A5; 107 PetscCall(MatCreate(comm, &nest)); 108 PetscCall(MatSetType(nest, MATNEST)); 109 PetscCall(MatNestSetSubMats(nest, 2, NULL, 2, NULL, mata)); 110 PetscCall(MatSetUp(nest)); 111 PetscCall(MatMatMult(nest, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C1)); 112 PetscCall(MatMatMult(nest, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C1)); 113 114 PetscCall(MatMatMultEqual(nest, B, C1, 10, &equal)); 115 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in C1 != C"); 116 PetscCall(MatDestroy(&C1)); 117 PetscCall(MatDestroy(&A5)); 118 PetscCall(MatDestroy(&nest)); 119 } 120 121 PetscCall(MatDestroy(&C)); 122 PetscCall(MatDestroy(&B)); 123 PetscCall(MatDestroy(&aij)); 124 PetscCall(MatDestroy(&A1)); 125 PetscCall(MatDestroy(&A2)); 126 PetscCall(MatDestroy(&A3)); 127 PetscCall(MatDestroy(&A4)); 128 PetscCall(PetscFinalize()); 129 return 0; 130 } 131 132 /*TEST 133 134 test: 135 nsize: 2 136 137 test: 138 suffix: 2 139 140 TEST*/ 141