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 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, test = PETSC_FALSE, test_null = PETSC_FALSE; 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 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_null", &test_null, NULL)); 71 if (test_null) { 72 IS is[2]; 73 PetscInt st; 74 75 PetscCall(MatGetLocalSize(A1, &m, NULL)); 76 PetscCall(MatGetOwnershipRange(A1, &st, NULL)); 77 PetscCall(ISCreateStride(comm, m, st, 2, &is[0])); 78 PetscCall(ISCreateStride(comm, m, st + 1, 2, &is[1])); 79 PetscCall(MatNestSetSubMats(nest, 2, is, 2, is, NULL)); 80 PetscCall(ISDestroy(&is[0])); 81 PetscCall(ISDestroy(&is[1])); 82 } else { 83 Mat mata[] = {A1, A2, A3, A4}; 84 85 PetscCall(MatNestSetSubMats(nest, 2, NULL, 2, NULL, mata)); 86 } 87 PetscCall(MatSetUp(nest)); 88 89 /* test matrix product error messages */ 90 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_nest*nest", &test, NULL)); 91 if (test) PetscCall(MatMatMult(nest, nest, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C)); 92 93 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matproductsetfromoptions", &test, NULL)); 94 if (test) { 95 PetscCall(MatProductCreate(nest, nest, NULL, &C)); 96 PetscCall(MatProductSetType(C, MATPRODUCT_AB)); 97 PetscCall(MatProductSymbolic(C)); 98 } 99 100 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matproductsymbolic", &test, NULL)); 101 if (test) { 102 PetscCall(MatProductCreate(nest, nest, NULL, &C)); 103 PetscCall(MatProductSetType(C, MATPRODUCT_AB)); 104 PetscCall(MatProductSetFromOptions(C)); 105 PetscCall(MatProductSymbolic(C)); 106 } 107 108 PetscCall(MatConvert(nest, MATAIJ, MAT_INITIAL_MATRIX, &aij)); 109 PetscCall(MatView(aij, PETSC_VIEWER_STDOUT_WORLD)); 110 111 /* create a dense matrix */ 112 PetscCall(MatGetSize(nest, &M, NULL)); 113 PetscCall(MatGetLocalSize(nest, &m, NULL)); 114 PetscCall(MatCreateDense(comm, m, PETSC_DECIDE, M, K, NULL, &B)); 115 PetscCall(MatSetRandom(B, NULL)); 116 117 /* C = nest*B_dense */ 118 PetscCall(MatMatMult(nest, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C)); 119 PetscCall(MatMatMult(nest, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C)); 120 PetscCall(MatMatMultEqual(nest, B, C, 10, &equal)); 121 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in C != nest*B_dense"); 122 123 /* Test B = nest*C, reuse C and B with MatProductCreateWithMat() */ 124 /* C has been obtained from nest*B. Clear internal data structures related to factors to prevent circular references */ 125 PetscCall(MatProductClear(C)); 126 PetscCall(MatProductCreateWithMat(nest, C, NULL, B)); 127 PetscCall(MatProductSetType(B, MATPRODUCT_AB)); 128 PetscCall(MatProductSetFromOptions(B)); 129 PetscCall(MatProductSymbolic(B)); 130 PetscCall(MatProductNumeric(B)); 131 PetscCall(MatMatMultEqual(nest, C, B, 10, &equal)); 132 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in B != nest*C_dense"); 133 PetscCall(MatConvert(nest, MATAIJ, MAT_INPLACE_MATRIX, &nest)); 134 PetscCall(MatEqual(nest, aij, &equal)); 135 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in aij != nest"); 136 137 /* Test with virtual block */ 138 PetscCall(MatCreateTranspose(A1, &A5)); /* A1 is symmetric */ 139 PetscCall(MatNestSetSubMat(nest, 0, 0, A5)); 140 PetscCall(MatMatMult(nest, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C1)); 141 PetscCall(MatMatMult(nest, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C1)); 142 PetscCall(MatMatMultEqual(nest, B, C1, 10, &equal)); 143 PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in C1 != C"); 144 PetscCall(MatDestroy(&C1)); 145 PetscCall(MatDestroy(&A5)); 146 147 PetscCall(MatDestroy(&nest)); 148 PetscCall(MatDestroy(&C)); 149 PetscCall(MatDestroy(&B)); 150 PetscCall(MatDestroy(&aij)); 151 PetscCall(MatDestroy(&A1)); 152 PetscCall(MatDestroy(&A2)); 153 PetscCall(MatDestroy(&A3)); 154 PetscCall(MatDestroy(&A4)); 155 PetscCall(PetscFinalize()); 156 return 0; 157 } 158 159 /*TEST 160 161 test: 162 diff_args: -j 163 nsize: 2 164 suffix: 1 165 166 test: 167 diff_args: -j 168 nsize: 2 169 suffix: 1_null 170 args: -test_null 171 172 test: 173 diff_args: -j 174 suffix: 2 175 176 test: 177 diff_args: -j 178 suffix: 2_null 179 args: -test_null 180 181 TEST*/ 182